根据定义式，可写出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复原出来的波形与原始波形看起来没有太大区别。

习题：

编写matlab脚本程序，实现DFT和IDFT算法；

利用所编写的程序对信号

进行DFT分析，并绘制频谱序列图，注意点数

和采样频率

的选择，避免出现频率混叠和栅栏效应；

利用所编写的程序对第2问的频率序列进行IDFT重构，并绘制重构之后的时间序列图。