数字信号处理实验.docx
《数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验.docx（27页珍藏版）》请在冰豆网上搜索。
数字信号处理实验
实验三快速傅里叶变换及其应用
实验中用到的信号序列：
高斯序列
0
Xa（n）=
1其他
衰减正弦序列：
0
Xb（n）=
1其他
三角波序列：
n0
Xc（n）=8-n4
0其他
反三角波序列：
4-n0
Xd（n）=n-44
0其他
1.观察高斯序列的时域和频域特性，固定信号中p=8,改变q值，使q=2,4,8，观察其时域频域特性，了解q取不同值时的影响；固定信号中q=8,改变p值，使p=8,13,14，观察其时域频域特性，了解q取不同值时的影响，注意p取何值时，会发生明显的泄漏现象，混叠是否也随之出现？
记录实验中观察到的现象，绘出相应的时域和频域特性曲线。
代码：
n=0:
15;
%p=8;q=2;
p=8;q=2;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,1）;
stem（n,x）;
title（'p=8,q=2,’时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,2）;
stem（n,y）;
title（'p=8,q=2,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
%p=8;q=4;
p=8;q=4;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,3）;
stem（n,x）;
title（'p=8,q=4,时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,4）;
stem（n,y）;
title（'p=8,q=4,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
%p=8;q=8;
p=8;q=8;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,5）;
stem（n,x）;
title（'p=8,q=8,时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,6）;
stem（n,y）;
title（'p=8,q=8,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
运行结果：
分析：
q越大，信号在时域上越宽，频域上越窄，不易混跌。
代码：
n=0:
15;
%p=8;q=8;
p=8;q=8;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,1）;
stem（n,x）;
title（'p=8,q=8,时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,2）;
stem（n,y）;
title（'p=8,q=8,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
%p=13;q=8;
p=13;q=8;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,3）;
stem（n,x）;
title（'p=13,q=8,时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,4）;
stem（n,y）;
title（'p=13,q=8,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
%p=14;q=8;
p=14;q=8;
x=exp（（-（n-p）.^2）/q）
y=abs（fft（x））;
subplot（3,2,5）;
stem（n,x）;
title（'p=14,q=8,时域特性’）;
xlabel（'n'）;
ylabel（'Xa'）;
subplot（3,2,6）;
stem（n,y）;
title（'p=14,q=8,频域特性’）;
xlabel（'w'）;
ylabel（'h'）;
运行结果：
分析：
P=14时，发生明显泄漏现象，在1/2频率处出现明显混叠。
2.观察衰减正弦序列的时域频域特性，a=0.1，f=0.625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375,0.5625，观察这两种情况下频谱的形状和谱峰的出现位置，有无混叠和泄漏现象？
说明产生现象的原因。
代码：
n=0:
15;
a=0.1;
%f=0.0625
f=0.0625;
x=exp（-a.*n）.*sin（2*pi.*f.*n）;
y=abs（fft（x））;
subplot（3,2,1）;
stem（n,x）;
xlabel（'n'）;
ylabel（'Xb'）;
title（'a=0.1;f=0.0625;时域特性'）;
subplot（3,2,2）;
stem（n,y）
xlabel（'w'）;
ylabel（'h'）;
title（'a=0.1;f=0.0625;频域特性'）;
%f=0.4375
f=0.4375;
x=exp（-a.*n）.*sin（2*pi.*f.*n）;
y=abs（fft（x））;
subplot（3,2,3）;
stem（n,x）;
xlabel（'n'）;
ylabel（'Xb'）;
title（'a=0.1;f=0.4375;时域特性'）;
subplot（3,2,4）;
stem（n,y）
xlabel（'w'）;
ylabel（'h'）;
title（'a=0.1;f=0.4375;频域特性'）;
%f=0.5625
f=0.5625;
x=exp（-a.*n）.*sin（2*pi.*f.*n）;
y=abs（fft（x））;
subplot（3,2,5）;
stem（n,x）;
xlabel（'n'）;
ylabel（'Xb'）;
title（'a=0.1;f=0.5625;时域特性'）;
subplot（3,2,6）;
stem（n,y）
xlabel（'w'）;
ylabel（'h'）;
title（'a=0.1;f=0.5625;频域特性'）;
运行结果：
分析：
f=0.4375时出现泄漏，原因：
加窗截断得过短；
f=0.5625时出现混叠和泄露，原因：
采样频率过低，加窗截断过短。
3.观察三角波和反三角波的时域频域特性，用N=8点FFT分析时域频域特性，观察两者的序列形状和频谱曲线有何异同？
在末尾补零，用N=32代入，观察发生的变化？
两种情况下频谱有相同之处吗？
说明了什么？
代码：
%三角波
N=8;
fori=1:
N
ifi<5
Xb（i）=i-1;
else
Xb（i）=8-i+1;
end
end
n=0:
N-1;
subplot（4,2,1）;
stem（n,Xb）;
xlabel（'n'）;
ylabel（'Xb'）;
title（'三角波时域特性'）;
subplot（4,2,3）;
Xk=fft（Xb,8）;
n=0:
N-1;
stem（n,abs（Xk））;
title（'8点频域特性'）;
xlabel（'w'）;
ylabel（'Xk'）;
N=8;
fori=1:
N
ifi<5
Xb（i）=i-1;
else
ifi<9
Xb（i）=8-i+1;
else
Xb（i）=0;
end
end
end
Xk=fft（Xb,32）;
n=0:
31;
subplot（4,2,4）;
stem（n,abs（Xk））;
title（'32点频域特性'）;
xlabel（'w'）;
ylabel（'Xk'）;
%反三角波
N=8;
fori=1:
N
ifi<5
Xb（i）=5-i;
else
Xb（i）=i-5;
end
end
n=0:
N-1;
subplot（4,2,5）;
stem（n,Xb）;
xlabel（'n'）;
ylabel（'Xb'）;
title（'反三角波时域特性'）;
subplot（4,2,7）;
Xk=fft（Xb,8）;
n=0:
N-1;
stem（n,abs（Xk））;
title（'8点频域特性'）;
xlabel（'w'）;
ylabel（'Xk'）;
N=8;
fori=1:
N
ifi<5
Xb（i）=5-i;
else
ifi<9
Xb（i）=i-5;
else
Xb（i）=0;
end
end
end
Xk=fft（Xb,32）;
n=0:
31;
subplot（4,2,8）;
stem（n,abs（Xk））;
title（'32点频域特性'）;
xlabel（'w'）;
ylabel（'Xk'）;
运行结果：
分析：
N=8时，由于采样点过少，出现栅栏效应，部分点被遮住；
对三角波影响不大，对反三角波影响很大，N=8时，很多重要频谱信息遗失。
适当补零可减少栅栏效应。
4.一个连续信号含两个频率分量，经采样得
X（n）=sin（2
）+cos（2
）n=0,1,……N-1
已知N=16,
分别为1/16,1/64，观察期频谱；当N=128，
，其结果有何不同？
为什么？
代码：
%N=16
N=16;
f=1/16;
fori=1:
N
Xb（i）=sin（2*pi*0.125*（i-1））+cos（2*pi*（0.125+f）*（i-1））;
end
Xk=fft（Xb,16）;
n=0:
15;
subplot（2,2,1）;
stem（n,abs（Xk））;
title（'16点频域特性,1/16'）;
xlabel（'w'）;
ylabel（'Xk'）;
N=16;
f=1/64;
fori=1:
N
Xb（i）=sin（2*pi*0.125*（i-1））+cos（2*pi*（0.125+f）*（i-1））;
end
Xk=fft（Xb,16）;
n=0:
15;
subplot（2,2,2）;
stem（n,abs（Xk））;
title（'16点频域特性,1/64'）;
xlabel（'w'）;
ylabel（'Xk'）;
%N=128
N=128;
f=1/16;
fori=1:
N
Xb（i）=sin（2*pi*0.125*（i-1））+cos（2*pi*（0.125+f）*（i-1））;
end
Xk=fft（Xb,128）;
n=0:
127;
subplot（2,2,3）;
stem（n,abs（Xk））;
title（'128点频域特性,1/16'）;
xlabel（'w'）;
ylabel（'Xk'）;
N=128;
f=1/64;
fori=1:
N
Xb（i）=sin（2*pi*0.125*（i-1））+cos（2*pi*（0.125+f）*（i-1））;
end
Xk=fft（Xb,128）;
n=0:
127;
subplot（2,2,4）;
stem（n,abs（Xk））;
title（'128点频域特性,1/64'）;
xlabel（'w'）;
ylabel（'Xk'）;
运行结果：
分析：
=1/16，N=16，N=64采样频谱均正确；
=1/64，N=16时采样频率过低，无法区分辩频率相差很小的两个信号，
N=64时采样频谱正确；频谱如实地反映了信号的特性。
5.用FFt计算Xa（n）（p=8,q=2）和Xb（n）（a=0.1,f=0.0625）的16点循环卷积和线性卷积。
代码：
p=8;
q=2;N=16;
fori=1:
N
Xa（i）=exp（-（（i-1-p）^2）/q）;
end
a=0.1;
f=0.0625;
fori=1:
N
Xb（i）=exp（-a.*i）.*sin（2.*pi.*f.*i）;
end
Xk=fft（Xa,16）;
Yk=fft（Xb,16）;
rm=real（ifft（Xk.*Yk））;
m=0:
N-1;
subplot（1,2,1）;
stem（m,rm）;
xlabel（'m'）;
ylabel（'rm'）;
title（'循环卷积'）;
Xk=fft（Xa,32）;
Yk=fft（Xb,32）;
rm=real（ifft（Xk.*Yk））;
m=0:
30;
subplot（1,2,2）;
stem（m,rm（1:
31））;
xlabel（'m'）;
ylabel（'rm'）;
title（'线性卷积'）;
运行结果：
6．产生一个512点的随机序列Xe（n），并用Xc（n）和Xe（n）线性卷积，观察卷积前后Xe（n）频谱的变化。
要求将Xe（n）分成8段，分别采用重叠相加法和重叠保留法。
重叠相加法：
代码：
Xe=randn（1,512）;
output（1:
700）=0;
temp2（1:
71）=0;
temp1（1:
71）=0;
N=8;
fori=1:
N
ifi<5
Xb（i）=i-1;
else
Xb（i）=8-i+1;
end
end
Xbk=fft（Xb,71）;
forj=1:
8
fork=1:
71
ifk<65
temp2（k）=Xe（k+64.*（j-1））;
else
temp2（k）=0;
end
end
tempk=fft（temp2,71）;
temp2=real（ifft（tempk.*Xbk））;
output（（64*（j-1）+1）:
64*j）=[temp1（8:
64）,temp1（65:
71）+temp2（1:
7）];
temp1=temp2;
ifj==8
output（（64*j+1）:
64*（j+1））=temp1（8:
71）;
end
end
output1（1:
519）=output（58:
576）;
m=0:
518;
Xek=fft（Xe,519）;
outputk=fft（output1,519）;
subplot（2,1,1）;stem（m,Xek（1:
519）,'.'）;title（'before'）;
subplot（2,1,2）;stem（m,outputk（1:
519）,'.'）;title（'after'）;
运行结果：
重叠保留法：
代码：
Xe=randn（1,512）;
output（1:
700）=0;
temp2（1:
71）=0;
temp1（1:
71）=0;
N=8;
fori=1:
N
ifi<5
Xb（i）=i-1;
else
Xb（i）=8-i+1;
end
end
Xbk=fft（Xb,71）;
forj=1:
8
fork=1:
71
ifk<8
ifj==1
temp2（k）=0;
else
temp2（k）=Xe（k+64.*（j-1）-8）;
end
else
ifj==9
temp2（k）=0;
else
temp2（k）=Xe（k-7+64.*（j-1））;
end
end
end
tempk=fft（temp2,71）;
temp2=real（ifft（tempk.*Xbk））;
output（（64*（j-1）+1）:
64*j）=[temp2（8:
71）];
end
output1（1:
519）=output（1:
519）;
m=0:
518;
Xek=fft（Xe,519）;
outputk=fft（output1,519）;plot（m,Xek（1:
519））
subplot（2,1,1）;stem（m,Xek（1:
519）,'.'）;title（'before'）;
subplot（2,1,2）;stem（m,outputk（1:
519）,'.'）;title（'after'）
运行结果：
7.用FFt计算Xa（n）（p=8,q=2）和Xb（n）（a=0.1,f=0.0625）的16点循环相关和线性相关，它们之间有何异同？
代码：
p=8;
q=2;
N=16;
fori=1:
N
Xa（i）=exp（-（（i-1-p）^2）/q）;
end
a=0.1;
f=0.0625;
fori=1:
N
Xb（i）=exp（-a.*i）.*sin（2.*pi.*f.*i）;
end
Xak=fft（Xa,N）;
Xbk=fft（Xb,N）;
rm=real（ifft（conj（Xak）.*Xbk））;
m=0:
N-1;
subplot（2,2,1）;
stem（m,rm）;
title（'循环相关x在前'）;
rm=real（ifft（conj（Xbk）.*Xak））;
m=0:
N-1;
subplot（2,2,2）;
stem（m,rm）;
title（'循环相关y在前'）;
Xak=fft（Xa,2*N）;
Xbk=fft（Xb,2*N）;
rm=real（ifft（conj（Xak）.*Xbk））;
rm=[rm（（N+2）:
2*N）,rm（1:
N）];
m=-N+1:
N-1;
subplot（2,2,3）;
stem（m,rm）;
title（'线性相关x在前'）;
rm=real（ifft（conj（Xbk）.*Xak））;
rm=[rm（（N+2）:
2*N）,rm（1:
N）];
m=-N+1:
N-1;
subplot（2,2,4）;
stem（m,rm）;
title（'线性相关y在前'）;
运行结果：
分析：
共有四种结果：
x在前与y在前结果左右颠倒；线性相关有部分零值区间。
8.用FFt计算Xa（n）（p=8,q=2）和Xb（n）（a=0.1,f=0.0625）的自相关函数。
代码：
p=8;
q=2;
N=16;
fori=1:
N
Xa（i）=exp（-（（i-1-p）^2）/q）;
end
a=0.1;
f=0.0625;
fori=1:
N
Xb（i）=exp（-a.*i）.*sin（2.*pi.*f.*i）;
end
Xak=fft（Xa,2*N）;
Xbk=fft（Xb,2*N）;
rm=real（ifft（conj（Xak）.*Xak））;
rm=[rm（（N+2）:
2*N）,rm（1:
N）];
m=-N+1:
N-1;
subplot（1,2,1）;
stem（m,rm）;
title（'Xa线性自相关'）;
p=8;
q=2;
N=16;
fori=1:
N
Xa（i）=exp（-（（i-1-p）^2）/q）;
end
a=0.1;
f=0.0625;
fori=1:
N
Xb（i）=exp（-a.*i）.*sin（2.*pi.*f.*i）;
end
Xak=fft（Xa,2*N）;
Xbk=fft（Xb,2*N）;
rm=real（ifft（conj（Xbk）.*Xbk））;
rm=[rm（（N+2）:
2*N）,rm（1:
N）];
m=-N+1:
N-1;
subplot（1,2,2）;
stem（m,rm）;
title（'Xb线性自相关'）;
运行结果：
思考题
1.实验中的信号序列Xc（n）、Xd（n），在单位圆上的z变换频谱会相同吗？
如果不同，说明那一个低频分量会多一点，为什么？
不同，三角波的低频分量多。
三角波在时域上信号没有突变，低频分量较多；
反三角波在n=0和n=7时有突变，因此其高频分量较多。
2.对一个有限长序列进行DFT等价于该序列周期延拓后进行DFS运算，因为DFS只是取其中一个周期进行计算，所以FFT在一定条件下也可以用以分析周期信号序列。
如果实正弦信号sin（2
n），f=0.1用16点FFt来做DFS运算，得到的频谱是其本身的的真实谱吗？
为什么？
不是。
f=0.1时，DFS运算在一个周期内应取10个频率点进行计算。
此时的FFT取16个点进行计算取出了一个周期内的10个点，但同时多取了6个有效点进行计算。
频谱不再正确。
应该取10个点，其余6个补零。
升级会员 