数字信号入门笔记4—离散傅里叶变换

Posted 2022-11-29 vonchenchen1

目录

数字信号入门笔记4—离散傅里叶变换

4.1离散傅里叶变换DFT定义

4.2 从典型复正弦信号看DFT的运作方式

4.3实正弦信号DFT

4.4采样不是正好一个周期复正弦信号 观察其DFT

4.4.1主瓣泄露

4.4.2旁瓣泄露

4.5 从等效滤波器的角度看DFT的运作方式

---

数字信号入门笔记4—离散傅里叶变换

本节知识点提要：

1. 从向量的角度看复正弦信号的DFT，直观解释DFT挑选频率的过程
2. 由复正弦信号DFT如何推广到正弦信号的DFT
3. 从不是一个恰好完整周期的信号DFT看频谱泄露与加窗
4. DFT结果实部、虚部、模长的意义和使用场景
5. 现实中不存在虚信号，真正要处理的是是处理完毕后信号的实部
6. 实正弦信号增加虚部后频谱不会出现对称，而是集中到真正的频率上
7. DFT的滤波器解释
8. DFT周期性与循环卷积
9. 脉冲信号的DFT与sinc函数
10. 复正弦信号与正弦信号的关系
11. 信号补零
12. 信号的分辨率

4.1离散傅里叶变换DFT定义

对于长度为N的离散时间信号x(n)，0≤n≤N-1，离散傅里叶变换DFT定义为

参考之前离散时间傅里叶变换的公式

离散傅里叶变换相当于对角频率w进行采样，原本连续的w被采样为2kPi/N，N个采样点映射到2Pi上。离散傅里叶变换也可以表示为幅度相位的形式

 

4.2 从典型复正弦信号看DFT的运作方式

将一个频率w=π/4，初始化相位为0，幅度为1的复正弦信号，采样8个点，对其做离散傅里叶变换。

对上述复正弦信号进行采样得到x(n)，第一幅图是matlab画出的x(n)，其实这里只是cos(angle)的幅值。 按照DFT的定义，对采样得到的x(n)序列做DFT运算，得到Xk(n)序列，第二幅图就是matlab绘制的Xk(n)序列。

N=8;

xn = (zeros(1,N));%sig exp(jwn)

w = 2*pi/N;

fi = 0;

for index=1:N

    angle = w*(index-1)+fi;

    xn(index)=cos(angle)+sin(angle)*1i;

end

subplot(2,1,1);stem(xn);grid;xlabel('n');ylabel('y');

Xk = (zeros(1,N));

for k=1:N

     opera = (zeros(1,N));   % exp(-j(2pi*k/N)n). %当前k对应的傅立叶算子序列 可以等效为当前k对应的滤波器，中心频率为(2pi*k/N)

    for n=1:N

        angle = -((2*pi*(k-1)/N)*(n-1)); %中心频率为2*pi*(k-1)/N的滤波器

        opera(n)=cos(angle)+sin(angle)*1i;        

        Xk(k) = Xk(k)+xn(n)*opera(n);

    end

end

subplot(2,1,2);stem(Xk);grid;xlabel('n');ylabel('y');

 

下面我们手推体会一下，下面时x(n) 序列的值

x(0) = e^(j*π/4)*0

x(1) = e^(j*π/4)*1

x(2) = e^(j*π/4)*2

x(3) = e^(j*π/4)*3

x(4) = e^(j*π/4)*4

x(5) = e^(j*π/4)*5

x(6) = e^(j*π/4)*6

x(7) = e^(j*π/4)*7

下面看一下k等于0时X(0)的推算

X(0) = (e^(j*π/4)*0)*1 + (e^(j*π/4)*1)*1 + (e^(j*π/4)*2)*1 + (e^(j*π/4)*3)*1 + (e^(j*π/4)*4)*1 + (e^(j*π/4)*5)*1 + (e^(j*π/4)*6)*1 + (e^(j*π/4)*7)*1

        =     向量0度      +    向量45度     +     向量90度        +     向量135度    +   向量180度     +    向量225度    +    向量270度     +     向量315度

        =     向量和发散 幅度为0

下面看一下k等于1时X(1)的推算

X(1) = (e^(j*π/4)*0)*(e^-(j*2π/8)*1*0) + (e^(j*π/4)*1)*(e^-(j*2π/8)*1*1) +
         向量0度     逆时针转0度                向量45度     逆时针转45度
       (e^(j*π/4)*2)*(e^-(j*2π/8)*1*2) + (e^(j*π/4)*3)*(e^-(j*2π/8)*1*3)
         向量90度     逆时针转90度         向量135度     逆时针转135度

       + (e^(j*π/4)*4)*(e^-(j*2π/8)*1*4) + (e^(j*π/4)*5)*(e^-(j*2π/8)*1*5) +
            向量180度   逆时针转180度           向量225度   逆时针转225度
         (e^(j*π/4)*6)*(e^-(j*2π/8)*1*6) + (e^(j*π/4)*7)**(e^-(j*2π/8)*1*7)
            向量270度   逆时针转270度        向量315度   逆时针转315度

    =    所有原始信号的向量都被旋转到0度，向量长度叠加幅度为8

下图蓝色箭头旋转代表w=π/4的信号旋转，红色箭头代表(e^-(j*2π/8)*k*n)，以w旋转的的向量，都恰好会被红色箭头旋转到0度。

 

这里可以看到，(e^-(j*2π/8)*k*n) 这一项相当于一个选择器，对于(e^(j*π/4)*n)在复平面旋转的信号（也就是角速度为π/4的信号），都能恰好转到0，可以理解为这样就挑出了信号x(n)中频率为π/4的信号。那么沿着这个逻辑，当k=2时，其实就是要挑出信号中是否有2π*2/8 也就是要挑出x(n)中是否包含角速度为π/2的信号，下面接着走一遍公式，看看结果。

X(2) = (e^(j*π/4)*0)*(e^-(j*2π/8)*2*0) + (e^(j*π/4)*1)*(e^-(j*2π/8)*2*1) +

           向量0度逆时针转0度=0度        向量45度逆时针转90度=-45度
       (e^(j*π/4)*2)*(e^-(j*2π/8)*2*2) + (e^(j*π/4)*3)*(e^-(j*2π/8)*2*3) +

           向量90度逆时针转180度=-90度   向量135度逆时针转270度=-135度
       +(e^(j*π/4)*4)*(e^-(j*2π/8)*2*4) +     (e^(j*π/4)*5)*(e^-(j*2π/8)*2*5) +      
          向量180度逆时针转360度=-180度   向量225度逆时针转450度=-225度   
       (e^(j*π/4)*6)*(e^-(j*2π/8)*2*6) +   (e^(j*π/4)*7)*(e^-(j*2π/8)*2*7) +

          向量270度逆时针转540度=-270度     向量315度逆时针转630度=-315度
    
    =    向量和发散 幅度为0

当k=2时，DFT挑选信号x(n)中2π*2/8频率的分量为0，说明x(n)中没有这个频率的信号。

这样就最特殊的复正弦信号和采样频率，看到了DFT是如何把信号的频率选择出来的。

4.3实正弦信号DFT

上文使用了最理想的复正弦信号，展现了一些DFT的工作原理，自然界中信号都是实信号，下面看一下实正弦信号的的DFT。这里没有找到特别形象的选择频率的过程，但是依然可以从向量的角度解释DFT的结果。

修改上文matlab代码中x(n)

%xn(index)=cos(angle)+sin(angle)*1i;

xn(index)=cos(angle);

DFT后结果如上图，在第2(k==1)点处仍有一条谱线，但是高度变为了4，而第8点(k==7)出也有一条高度为4的谱线。下面依然用向量的观点看一下k==1时向量的情况。由于每项变成了 实正弦信号x(e^-(j*2π/8)*k*n)的形式，那么就不能看做是把向量旋转，而是改变(e^-(j*2π/8)*k*n)向量的模长或者方向了。

x(0) = 1

x(1) = sqrt(2)/2

x(2) = 0

x(3) = -sqrt(2)/2

x(4) = -1

x(5) = -sqrt(2)/2

x(6) = 0

x(7) = sqrt(2)/2

当k==1时X(1)的推算如下

X(1) = 1*(e^-(j*2π/8)*1*0) +     (sqrt(2)/2)*(e^-(j*2π/8)*1*1) +

          0度 模长1                  -45度 模长(sqrt(2)/2)

       (0)*(e^-(j*2π/8)*1*2) +    (-sqrt(2)/2)*(e^-(j*2π/8)*1*3) +

          -90度 模长0                 -(-135)度=45度   模长(sqrt(2)/2)
       + (-1)*(e^-(j*2π/8)*1*4) +  (-sqrt(2)/2)*(e^-(j*2π/8)*1*5) +       

         -(-180)度=0度 模长1           -(-225)度=-45度 模长(sqrt(2)/2)
       (0)*(e^-(j*2π/8)*1*6) +      (sqrt(2)/2)*(e^-(j*2π/8)*1*7) +

         -270度 模长0                 -315度=45度  模长(sqrt(2)/2)

     =   向量和为4

这里有个问题，为什么同样的频率，复正弦信号结果恰好可以定位到对应的频率，而实正弦信号的结果却是两根共轭对称的谱线？这个问题用向量的观点可以推出结果，但是似乎无法直接说明。下面看一下简单的数学推导，这里可以得出一个结论，就是谱线共轭对称和傅里叶变换的周期性相关。

DFT  X(k)中第k个点需要乘一个 e^-j(2πnk/N)，而第N-k个点需要乘e^-j(2πn(N-k)/N)，而e^-j(2πn(N-k)/N) -> e^-jn(2π-2πk/N) 这里根据e^jw的周期性(2π-2πk/N)==(2πk/N)。可以看出第k个点对应乘e^-j(2πnk/N)，第N-k对应乘e^-jn(2π-2πk/N)，所以这个值是关于π对称的，我猜想DFT的周期性也就是因为e^-j(2πnk/N)的周期性。而x(n)为实数，求每个X(k)时都是实数序列x(n)与e^-j(2πnk/N)做卷积，所以X(k)与X(N-k)的关系是由e^-j(2πnk/N)决定，为关于π对称。

还有一种比较民科的看法，由于实信号在极坐标上实轴，可以试想复正弦信号绕一圈，实轴上的轨迹，先从1到0，从0到-1，再从-1到0，从0到1，整个过程是对称的。

4.4采样不是正好一个周期复正弦信号 观察其DFT

将上文代码中w改为w = pi/3，N=8;

这里看其DFT中包含了各个频率的信息。为什么原本输入频率是π/3，而频谱中却包含多种频率呢。上文知道DFT是有隐含的周期性。参考下图，如果对这个时域信号进行周期性延展，n=9的时候值是与n=1的值是一样的，也就是整个波形并非一个标准的正弦曲线，而且还会多出一个拐角，这样就会引入其他频率的信号，这样看截断后得到波形其实并非真实信号。

 

如下图，a恰好截断到，b截断后的波形经过延拓后则明显不是原先波形，而且产生一个突变，导致其他频率大量引入。

这样，真实频率的大小在频谱上会变小(主瓣泄露)，同时还多出其他频率(旁瓣泄露)，产生频谱泄露。上面是从信号DFT的角度来看的，下面从时域截断的角度来看一下泄露的产生。

4.4.1主瓣泄露

这里需要了解一下脉冲信号，信号的截断相当于脉冲信号与输入信号时域上每个点相乘得到的序列，由卷积定理可知时域相乘相当于频域卷积，那么首先需要知道脉冲信号的傅里叶变换结果。可以回顾第一章理想低通滤波器，这里需要记住矩形框和sinc互为傅里叶变换对。

单频复正弦信号在频谱上为一根直线，现在用sinc与单频复正弦信号的傅里叶变换做卷积，如下图

 

这里图a是复正弦信号的傅里叶变换，b是矩形信号的傅里叶变换，两个信号在频域卷积就得到c图，这样就从频域的角度解释了为什么单频复正弦信号的傅里叶变换为什么频谱中有多个频率。

xn = (zeros(1,N));

xn(3) = 1;

xn(4) = 1;

Xn = fft(xn);

subplot(2,1,1);stem(xn, 'filled');grid;xlabel('n');ylabel('y');

subplot(2,1,2);stem(Xn, 'filled');grid;xlabel('n');ylabel('y');

 

上面两图代表不同宽度的矩形脉冲和其对应的频谱，可以看到矩形越宽，旁瓣越小，能量越集中在主瓣。对应时域采样就是做DFT时要多采点，相当于增加脉冲宽度，抑制主瓣泄露。

4.4.2旁瓣泄露

产生旁瓣泄露是由于截断时截断处产生了突变，导致引入了高频信号。从时域上看，如果要避免这些高频引入，必须使截断后信号边缘变缓，减少突变。在处理时常采用时域加窗的方式。下面是几种常见的窗函数

窗函数时域波形

窗函数频域波形

 

这里可以看到，主瓣变窄和旁瓣变低是不能同时满足的，需要根据实际信号的特性选择不同的窗函数。

下面我们尝试加窗，比对加窗前后变化

N=8;

xn = (zeros(1,N));%sig exp(jwn)

w = pi/3;

fi = 0;

for index=1:N

    angle = w*(index-1)+fi;

    xn(index)=cos(angle)+sin(angle)*1i;

end

subplot(2,2,1);stem(xn);grid;xlabel('n');ylabel('y');

Xk = (zeros(1,N));

for k=1:N

    opera = (zeros(1,N));

    for n=1:N

        angle = -((2*pi*(k-1)/N)*(n-1));

        opera(n)=cos(angle)+sin(angle)*1i;

        Xk(k) = Xk(k)+xn(n)*opera(n);

    end

end

subplot(2,2,2);stem(Xk);grid;xlabel('k');ylabel('y');

w_ham = hamming(N);w_ham = w_ham.';

xn = xn.*w_ham;

subplot(2,2,3);stem(xn);grid;xlabel('n');ylabel('y');

Xk = (zeros(1,N));

for k=1:N

    opera = (zeros(1,N));

    for n=1:N  

        angle = -((2*pi*(k-1)/N)*(n-1));

        opera(n)=cos(angle)+sin(angle)*1i;

        Xk(k) = Xk(k)+xn(n)*opera(n);

    end

end

subplot(2,2,4);stem(Xk);grid;xlabel('k');ylabel('y');

可见左边两幅图加窗后，左下角加窗后边缘值变小，波形明显畸变。DFT后频谱k=2处能量集中，高频分量相比原始信号的频谱减小很多。

4.5 从等效滤波器的角度看DFT的运作方式

前面计算傅里叶变换，对于不同的k，输入都是x(n)序列，与e^-j(2πnk/N)序列卷积，得到X(k)。这样其实就相当于序列通过滤波器e^-j(2πnk/N)。

 

现在把 e^-j(2πnk/N) 当做一个系统，计算频率响应，看一下这个滤波器的特性。

 

为了得到稳定系统，用如下公式近似上述滤波器

下面计算这个系统的频率响应

这里取不同的k值，幅频响应形状不变，为sinc，中心频率随k变化。这样DFT就可以等效于有k个不同的滤波器，每次x(n)与之卷积，得到X(k)，做k次卷积得到X(k)序列。

频域相乘

如果采样点在两个滤波器中心频率之间，则会被分散给两个主瓣。

这里再转到向量，采样点在两个滤波器中心频率之间，再单位圆上就是相当于在两个频率点之间。4.1.1中恰好在中心频率，所以经过旋转后都在中心频率上，向量叠加后结果也比较整齐，要么都在实轴，要么发散为0。如果在中心频率之间，每次旋转完毕后都会有个小偏转角度。这样向量和就会每次都旋转一些，如下图

这里就从向量的角度解释了为什么不在滤波器中心频率，最终得到的幅值会变小。