根据定义式，可写出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重构，并绘制重构之后的时间序列图。
