一、实验目的
1.掌握离散傅里叶变换（DFT）及快速傅里叶变换（FFT）的计算机实现方法。
2.检验序列DFT的性质。
3.掌握利用DFT（FFT）计算序列线性卷积的方法。
4.学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。
5.了解采样频率对谱分析的影响。
6.了解利用FFT进行语音信号分析的方法。
二、实验设备
1.计算机
2.Matlab软件7.0以上版本。
三、实验内容
1.对不同序列进行离散傅里叶变换并进行分析；DFT共轭对称性质的应用（通过1次N点FFT计算2个N点实序列的DFT）。
2.线性卷积及循环卷积的关系，以及利用DFT（FFT）进行线性卷积的方法。
3.比较计算序列的DFT和FFT的运算时间。
4.利用FFT实现带噪信号检测。
5.利用FFT计算信号频谱及功率谱。
6.扩展部分主要是关于离散系统采样频率、时域持续时间、谱分辨率等参数之间的关系，频谱的内插恢复，对语音信号进行简单分析。
四、实验原理
1.序列的离散傅里叶变换及性质
离散傅里叶变换的定义：
离散傅里叶变换的性质：
（1）DFT的共轭对称性。若，，则：， 。
（2）实序列DFT的性质。若为实序列，则其离散傅里叶变换为共轭对称，即。
（3）实偶序列DFT的性质。若为实偶序列，则其离散傅里叶变换为实偶对称，即。
（4）实奇序列DFT的性质。若为实奇序列，则其离散傅里叶变换为纯虚奇对称，即。
2.利用DFT计算线性卷积
时域循环卷积定理：设有限长序列为和，它们的N点DFT分别为和，如果，则其IDFT为两序列的循环卷积。
循环卷积和线性卷积的关系：循环卷积是线性卷积周期延拓的主值序列。当循环卷积的点数大于等于线性卷积的长度，这时循环卷积等于线性卷积，就可以利用DFT计算线性卷积，方法如图4-1所示。
L点
DFT
L点IDFT
x1(n)
y(n)
补L-N
个零点
补L-M
个零点
L点
DFT
x2n)
图4-1利用DFT计算线性卷积流程图
3、利用DFT对信号进行谱分析
（1）用DFT进行谱分析的参数选择
用DFT进行连续信号谱分析时参数的选择要遵循一定的原则。
， ，
（2）用DFT进行谱分析的误差问题
由于DFT计算频谱仅限于离散点上的频谱，而采样点之间的频谱函数值是不知道的。因此会引起栅栏效应，解决方法是：对有限长序列在原数据末尾补0，对无限长序列增大截取长度和DFT变换区间长度。由于实际中信号可能是无限长的，而DFT处理的是有限长序列，因此实际中经常要把观测的时域信号截断，相当于时域加窗，会引起截断效应，表现为频谱泄露和谱间干扰，解决的方法是增加窗长、改变窗函数形状。栅栏效应和谱分辨率是不同的概念，对信号补0，可以减小栅栏效应，但是由于截断已经使频谱变模糊了，所以补零后虽采样间隔变小，但得到的频谱包络仍是变模糊的频谱，因此频率分辨率没有提高，要提高频率分辨率，必须使时域截取长度增加。
4.DFT和FFT的运算量
离散傅里叶变换在实际应用中非常重要，利用它可以计算信号的频谱、功率谱和线性卷积等。但其计算量太大(与N的平方成正比）, 很难实时地处理问题, 因此引出了快速傅里叶变换(FFT) 。
设序列的长度为N，则计算DFT和计算FFT的运算量如表4.1所示。
表4.1 DFT和FFT运算量比较
| 计算DFT | 计算FFT | |
| 复数乘法次数 | ||
| 复数加法次数 |
5.周期图法计算信号功率谱
周期图法是一种估计信号功率谱密度的方法。因为序列的DFT隐含有周期性，所以这种功率谱也有周期性，称为周期图。该方法的优点是能应用离散傅里叶变换的快速算法来进行估值。定义如式（4-1）所示。
（4-1）
6.信噪比
设纯净信号为，噪声信号为，带噪信号为，则信噪比的定义如式（4-2）所示，单位为dB。
（4-2）
五、实验步骤
1.序列的离散傅里叶变换及分析
分别对复序列，实序列，实偶序列，实奇序列，虚奇序列进行离散傅里叶变换，得到实验结果并对其特点进行分析。实验所需序列自选。
x1n=[1 1 1 1];
x2n=[1 2 3 4];
xn=x1n+1i*x2n;
Xk=fft(xn,4);
k=0:3;wk=2*k/4;
subplot(3,2,1);
h=stem(wk,abs(Xk),'o','fill');
set(h,'LineWidth',3)
title('复序列的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk1=fft(x1n,4);
k=0:3;wk=2*k/4;
subplot(3,2,2);
h=stem(wk,abs(Xk1),'o','fill');
set(h,'LineWidth',3)
title('实序列的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk2=fft(x2n,4);
k=0:3;wk=2*k/4;
subplot(3,2,1);
h=stem(wk,abs(Xk2),'o','fill');
set(h,'LineWidth',3)
title('复序列的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk3=(0.5)*(Xk+conj(Xk(mod(-k,4)+1)));
k=0:3;wk=2*k/4;
subplot(3,2,3);
h=stem(wk,abs(Xk3),'o','fill');
set(h,'LineWidth',3)
title('实偶序列的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk4=(0.5)*(Xk-conj(Xk(mod(-k,4)+1)));
k=0:3;wk=2*k/4;
subplot(3,2,4);
h=stem(wk,abs(Xk4),'o','fill');
set(h,'LineWidth',3)
title('实奇序列的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk5=-1i*(0.5)*(Xk-conj(Xk(mod(-k,4)+1)));
k=0:3;wk=2*k/4;
subplot(3,2,5);
h=stem(wk,abs(Xk5),'o','fill');
set(h,'LineWidth',3)
title('虚奇序列的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
2.利用共轭对称性，设计高效算法计算2个N点实序列的DFT。
用一个N点FFT计算两个长度为N的实序列N点离散傅里叶变换，并将结果和直接使用两个N点DFT得到的结果进行比较。
x1n=[1 1 1 1];
x2n=[1 2 3 4];
xn=x1n+1i*x2n;
Xk=fft(xn,4);
k=0:3;wk=2*k/4;
subplot(3,2,1);
h=stem(wk,abs(Xk),'o','fill');
set(h,'LineWidth',3)
title('xn的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk3=(0.5)*(Xk+conj(Xk(mod(-k,4)+1)));
k=0:3;wk=2*k/4;
subplot(3,2,2);
h=stem(wk,abs(Xk3),'o','fill');
set(h,'LineWidth',3)
title('经计算的X1k的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk4=-1i*(0.5)*(Xk-conj(Xk(mod(-k,4)+1)));
k=0:3;wk=2*k/4;
subplot(3,2,3);
h=stem(wk,abs(Xk4),'o','fill');
set(h,'LineWidth',3)
title('经计算的X2k的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk1=fft(x1n,4);
k=0:3;wk=2*k/4;
subplot(3,2,4);
h=stem(wk,abs(Xk1),'o','fill');
set(h,'LineWidth',3)
title('x1n的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
Xk2=fft(x2n,4);
k=0:3;wk=2*k/4;
subplot(3,2,5);
h=stem(wk,abs(Xk2),'o','fill');
set(h,'LineWidth',3)
title('x2n的4点DFT的幅频特性图');
xlabel('\\omega/\\pi');
ylabel('幅度');
3.线性卷积及循环卷积的实现及二者关系分析
计算两序列的线性卷积及循环卷积，循环卷积采用2种计算方法（时域、频域方法）。设序列x1长度为M，序列x2长度为N，循环卷积长度为L，分别计算L大于、等于、小于（M+N-1）时的循环卷积。序列x1、x2、L自选，得到实验结果并对线性卷积及循环卷积的关系进行分析。
x1=[1,2,3,4,5,6];
x2=[1,2,3,4,5];
N=12;
y=circonv(x1,x2,N);
subplot(3,1,1)
stem(x1)
subplot(3,1,2)
stem(x2)
subplot(3,1,3)
stem(y)
k1=1:6;
f1=[1,2,3,4,5,6];
k2=1:5;
f2=[1,2,3,4,5];
[f,k]= lsjuanji (f1,f2,k1,k2);
subplot(3,1,1)
stem(f1)
subplot(3,1,2)
stem(f2)
subplot(3,1,3)
stem(f)
线性卷积
4.比较DFT和FFT的运算时间。
（1）自行选择进行计算的序列（或产生随机序列）。
（2）利用计时函数tic、toc，计算点数N=、128、256、……
tic
N=256
n=0:N-1;
x=2*cos(0.25*pi*n)+cos(0.65*pi*n);
Y=dft(x,N);
toc
N =
时间已过 0.100923 秒。
>> Untitled
N =
128
时间已过 0.016722 秒。
>> Untitled
N =
256
时间已过 0.197153 秒。
tic
N=256;
n=0:N-1;
x=2*cos(0.25*pi*n)+cos(0.65*pi*n);
Y=fft(x,N);
toc
时间已过 0.920454 秒。
>> Untitled
时间已过 0.000888 秒。
>> Untitled
时间已过 0.002787 秒。
5.利用FFT计算信号功率谱。
自行设计信号，利用周期图法实现信号功率谱并画图。
N=0:1:63;
x1=[1 2 3 4 5];
x2=[5 4 3 2 1];
xk1=x1+1i*x2;
xk2=x1-1i*x2;
Xk1=fft(xk1,);
Xk2=fft(xk2,);
P=1/.*Xk1.*Xk2;
subplot(2,1,1)
stem(N,abs(P))
subplot(2,1,2)
stem(N,angle(P))
6.利用FFT求信号频谱及分析采样频率对频谱的影响。
设模拟信号，以 进行取样，求N点DFT的幅值谱（N分别取50，256），利用FFT实现。分析信号频谱所对应频率轴的数字频率和频率之间的关系。在截取长度不变的条件下改变采样频率，观察信号频谱的变化，分析产生变化的原因。
N=256;
n=0:N-1;
t=0.01*n;
x=3*cos(8*pi*t)+6*cos(20*pi*t);
X=fft(x,N);
stem(n,x)
7.利用FFT分析信号频率成分。
对信号加入白噪声（0.2*randn(1,N)），比较有无噪声时的信号谱，加大噪声到2*randn(1,N)和10*randn(1,N)，计算信噪比，画出并比较不同噪声下时域波形和频谱，讨论噪声对信号分析的影响。信号自行选择。
程序：
w=randn(1,80);
n=1:80;
x=cos(0.04*pi*n);
subplot(2,2,1);
stem(n,x);
title('余弦信号cos(0.04*pi*n)');
xlabel('2*pi\\\\omega/');
ylabel('幅度')
y=(cos(0.04*pi*n)+(0.2)*w(1,n))*1;
subplot(2,2,2);
stem(n,y);
title('带噪信号cos(0.04*pi*n)+(0.2)*w(1,n)');
xlabel('2*pi\\\\omega/');
ylabel('幅度')
y1=(cos(0.04*pi*n)+(2)*w(1,n))*1;
subplot(2,2,3);
stem(n,y1);
title('带噪信号cos(0.04*pi*n)+(2)*w(1,n)');
xlabel('2*pi\\\\omega/');
ylabel('幅度')
y2=(cos(0.04*pi*n)+(10)*w(1,n))*1;
subplot(2,2,4);
stem(n,y2);
title('带噪信号cos(0.04*pi*n)+(10)*w(1,n)');
xlabel('2*pi\\\\omega/');
ylabel('幅度')
8.创新训练拓展内容
（1）信号持续时间、频谱分析范围、采样点数和谱分辨率的关系。
1）已知信号最高频率，谱分辨率，通过程序确定采样频率，采样点数，时域最小记录时间（信号持续时间），时域最大采样间隔。
设F（谱分辨率）≤10Hz，信号最高频率fc=2.5kHz.
Tp（最小记录时间）≥1/F=1/10=0.1s
Tmax（最大采样间隔）=1/Fsmin=1/(2fc)=1/(2*2500)=0.2*10-3
Nmin(最小采样点数)=2fc/F=2*2500/10=500
Fft进行频谱分析，需要的点数为2的n次幂，距20最近的2的整数次幂为2^5=32
2）改变时域记录时间（小于，大于最小记录时间），分析相应信号谱分辨率不同。
F≥1/Tp
3）若要提高谱分辨率或谱分辨率不满足要求，而又要保持谱分析范围不变，采用适当的方法实现上述要求。
提高谱分辨率，最小记录时间变大（Tp≥1/F）。若fc（最高频率）不变，则最少采样点数增加。
（2）频谱的内插函数恢复方法。
利用频域采样及内插公式实现离散谱到连续谱的转换。
clear;
N=;
Fs=6000;
xt=zeros(1,401);
for m=-N/2:N/2-1
xn(m+N/2+1)=exp(-1000*abs(m/Fs));
for k=0:400
t=(k-200)/30000;
if pi*(t-m/Fs)*Fs==0
xt(k+1)=xt(k+1)+exp(-1000*abs(m/Fs));
else
xt(k+1)=xt(k+1)+exp(-1000*abs(m/Fs))*sin(pi*(t-m/Fs)*Fs)./(pi*(t-m/Fs)*Fs);
end;
end;
end
figure(1)
subplot(2,1,1)
stem((0:N-1)-N/2,xn);
title('离散')
subplot(2,1,2)
plot(xt)
title('连续')
（3）对语音信号进行简单分析。
对采集（或给定）的语音信号进行谱分析，获得采样频率，画出语音信号的波形及频谱，并分析语音信号的频率分布特点。
nChannels=1;
Fs=4000;
nBits=8;
recObj = audiorecorder(Fs,nBits,nChannels);
disp('Start speaking.')
recordblocking(recObj, 5);
disp('End of Recording.');
% 回放录音数据
play(recObj);
% 获取录音数据
myRecording = getaudiodata(recObj);
% 绘制录音数据波形
figure(1)
plot(myRecording);
%画出语音信号频谱
figure(2)
Y=fft(myRecording);
plot(abs(Y))
