1周学FFT——第2天 DFT和IDFT的MATLAB实现

根据定义式，可写出DFT的MATLAB代码如下[从玉良，2009，p72]：

function [f, Xk] = mydft(xn, fs, N)
    % DFT
    n = [0:1:N-1];
    k = n;
    WN= exp(-j*2*pi/N);
    nk = n' * k;            % N^2 times multiply
    Xk = xn(1:N) * WN.^nk;  % N^3 times multiply
    f = 0:fs/N:fs/N*(N-1);
% end of function

该函数的使用方法如下：

clear;
fs = 200;   % sample frequency
N  = 40;    % the point of DFT, try to change it.

% compose the input signal
d = 0:(1/fs):1;
xn = sin(2*pi*5*d) + 0.3*sin(2*pi*10*d);
subplot(2,1,1);
stem(d, xn);
title('离散时间序列x(n)，(1, 5Hz)和(0.3, 10Hz)叠加，采样频率200Hz');

[f, Xk] = mydft(xn, fs, N);
Xk_amp = abs(Xk)/N*2;   % amplitude of Xk, should be divided by N/2
subplot(2,1,2);
stem(f(1:N/2), Xk_amp(1:N/2));  % display only half of N
title('离散频率序列X(k)，N=40，\Delta f = fs/N=5Hz，避免栅栏效应')

运行结果为：

DFT结果

如果预先知道信号的基频及谐波频率的话，可以有针对性的设计fs和N，以避免频率混叠和栅栏效应。

比如上面的例子中信号基频为5Hz，含有2次谐波（10Hz），则可令fs=200，远大于奈奎斯特频率，可避免频率混叠。同时，N=40，则细分频率，正好可整除基波及谐波，可避免出现栅栏效应。

如果fs与N选择不合适的话，效果如下：

频率混叠和栅栏效应

IDFT的MATLAB代码如下：

function [nn, xn] = myidft(Xk, fs, N)
    % DFT
    n = [0:1:N-1];
    k = n;
    WN= exp(-j*2*pi/N);
    nk = n' * k;            % N^2 times multiply
    xn = (Xk(1:N) * WN.^(-nk))/N;  % N^3 times multiply
    nn = 0:1/fs:1/fs*(N-1);
% end of function

该函数的使用方法如下：

%% IDFT
clear;
fs = 200;   % sample frequency
N  = 40;    % the point of DFT, try to change it.

% compose signal
d = 0:(1/fs):1;
xn = sin(2*pi*5*d) + 0.3*sin(2*pi*10*d);
subplot(2,1,1);
stem(d, xn);
title('原始信号序列');

[f, Xk] = mydft(xn, fs, N);
Xk_amp = abs(Xk)/N*2;   % amplitude of Xk, should be divided by N/2

[nn, xnn] = myidft(Xk, fs, N);
subplot(2,1,2);
stem([nn, 1/5+nn, 2/5+nn, 3/5+nn, 4/5+nn],[xnn, xnn, xnn, xnn, xnn], 'r');
title('IDFT复原序列（复制5次）');

运行效果如下：

IDFT运行效果

如果此处刻意的制造栅栏效应（比如令N=100），根据栅栏后的频谱IDFT复原出来的波形与原始波形看起来没有太大区别。

习题：

1. 编写matlab脚本程序，实现DFT和IDFT算法；
2. 利用所编写的程序对信号进行DFT分析，并绘制频谱序列图，注意点数和采样频率的选择，避免出现频率混叠和栅栏效应；
3. 利用所编写的程序对第2问的频率序列进行IDFT重构，并绘制重构之后的时间序列图。