T=A（J+1）;
A（J+1）=A（I+1）;
A（I+1）=T;
end
K=NV2;
whileK<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
End
A
程序2：
%基2按时间抽选
forM=1:
3%将DFT做m次基2分解，从左到右，对每次分解作DFT运算
Nmr=2^M;
U=1;%旋转因子u初始化
WN=exp（-j*2*pi/Nmr）;%本次分解的基本DFT因子WN＝exp（-i*2*pi/Nmr）
forn=1:
Nmr/2%本次跨越间隔内的各次碟形运算
fork=n:
Nmr:
N%本次碟形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2;%确定碟形运算的对应单元下标（对称性）
t=A（kp）*U;%碟形运算的乘积项
A（kp）=A（k）-t;%碟形运算的加法项
A（k）=A（k）+t;
end
U=U*WN;%修改旋转因子，多乘一个基本DFT因子WN
end
end
（2）调试程序，运行结果：
程序1结果：
输入倒位序，输出自然序
程序2调试结果：
按时间抽选基2-FFT
分析与总结：
（1）输入是倒位序0、4、2、6、1、5、3、7，输出自然序1、2、3、4、5、6、7，在编写程序时注意换址的处理，否则程序不能成功。
（2）进行一次蝶形运算需要
次复数乘法和N次复数加法，因此一个人完整的基-2FFT算法需要
次复数乘法和
次复数加法，本来应该需要N^2次复数乘法和加法，这样就大大减少了运算量。
2.已知周期信号为
，试用DFT法分析其频谱。
（1）程序：
N=input（'N='）;T=0.025;n=1:
N;D=2*pi/（N*T）;%原始数据
xa=cos（20*n*T）;Xa=T*fftshift（fft（xa,N））;%计算分辨率
Xa
（1）;k=floor（-（N-1）/2:
（N-1）/2）;
TITLLE=sprintf（'N=%i,L=%i',N,N*T）;
plot（k*D,abs（Xa））;axis（[-20,20,0,max（abs（Xa））+2]）;
xlabel（'\omega'）;ylabel（'|X（j\omega）|'）;title（TITLE）;
（2）程序流图
（3）
输入原始数据
计算分辨率
求DFT，移位
绘制频谱图图
（3）运行结果
N=200,L=5
N=400,L=10
分析与总结：
此信号的最高频率为fh=20\（2*pi）,因此抽样频率fs>=2fh=20/pi，即取抽样时间间隔T=1/fs<=pi/20，由于截断的影响，频谱会扩散，所以需要采用更高的采样频率，取T=0.025。
因为抽样后的序列不是周期性序列，因而截断长度不可能恰好是周期的整数倍，所以不能得到准确的频谱峰值，只能去较大的L=NT，使得到数据更多一些，频谱更接近真实情况。
3、实验五
1.语音信号的采集
（1）程序：
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;figure
（1）;
plot（y）;title（'原时域波形'）;
ylabel（'amplitude'）;xlabel（'n'）;
figure
（2）;plot（abs（Y））;
axis（[0,4096,0,3]）;title（'原频谱特性'）;
ylabel（'amplitude'）;xlabel（'frequency\HZ'）;
（2）结果图：
分析与总结：
注意wavread函数的使用方法，读入wav格式的语音信号，该语音信号可由cooledit软件将mp3格式的文件转换为wav格式的，刚开始编译程序时老是命令窗口老是显示我的文件不能被打开，后来经过查找各种资料才知道，wav格式也可以分成好几种，matlab只能读取cooledit中windowspcm.wav格式的语音，如果保存成其他wav格式的则读取时会出现error。
编译出现波形的那一刹那我真的特别开心！
用滤波器对信号滤波，比较滤波前后语音信号波形频谱
2.1巴特沃斯低通滤波器滤波
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;
FS=8000;Fl=1000;Fh=1200;
wp=（Fl/FS）*2*pi;ws=（Fh/FS）*2*pi;
OmegaP=2*FS*tan（wp/2）;OmegaS=2*FS*tan（ws/2）;
[n,wn]=buttord（OmegaP,OmegaS,1,100,'s'）;
[b,a]=butter（n,wn,'s'）;[bz,az]=bilinear（b,a,FS）;
freqz（bz,az,4096,FS,'whole'）;
x=filter（b,a,y）;X=fft（x,4096）;figure
（2）;
sublpo（211）,plot（abs（Y））;axis（[0,4096,0,3]）;
title（'滤波器前信号'）;
sublpoy（211）,plot（abs（Y））;axis（[0,4096,0,3]）;
title（'双线性巴特沃斯低通滤波器信号'）;
2.2巴特沃斯低通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;
FS=8000;Fl=1000;Fh=1200;
wp=（Fl/FS）*2*pi;ws=（Fh/FS）*2*pi;
OmegaP=2*FS*tan（wp/2）;OmegaS=2*FS*tan（ws/2）;
[n,wn]=buttord（OmegaP,OmegaS,1,100,'s'）;
[b,a]=butter（n,wn,'s'）;[bz,az]=bilinear（b,a,FS）;
freqz（bz,az,4096,FS,'whole'）;
x=filter（b,a,y）;X=fft（x,4096）;figure
（2）;
subplot（211）,plot（abs（Y））;axis（[0,4096,0,3]）;
title（'双线性巴特沃斯低通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4096,0,3]）;
title（'双线性巴特沃斯低通滤波后的信号频谱'）;
2.3巴特沃斯高通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;FS=000;
Fl=4800;Fh=5000;
wp=（Fh/FS）*2*pi;ws=（Fl/FS）*2*pi;
OmegaP=2*FS*tan（wp/2）;OmegaS=2*FS*tan（ws/2）;
[n,wn]=buttord（OmegaP,OmegaS,1,100,'s'）;
[b,a]=butter（n,wn,'s'）;[bz,az]=bilinear（b,a,FS）;
freqz（bz,az,4096,FS,'whole'）;
x=filter（b,a,y）;X=fft（x,4096）;figure
（2）;
subplot（211）,plot（abs（Y））;axis（[0,4096,0,3]）;
title（'双线性巴特沃斯高通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4096,0,3]）;
title（'双线性巴特沃斯高通滤波后的信号频谱'）;sound（x,fs,bits）;
2.4切比雪夫带通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;
wp1=1200/8000*2*pi;wp2=3000/8000*2*pi;
ws1=1000/8000*2*pi;ws2=3200/8000*2*pi;
Rp=1;Rs=100;Wp1=tan（wp1/2）;Wp2=tan（wp2/2）;
Ws1=tan（ws1/2）;Ws2=tan（ws2/2）;
BW=Wp2-Wp1;W0=Wp1*Wp2;W00=sqrt（W0）;
WP=1;WS=WP*（W0^2-Ws1^2）/（Ws1*BW）;
[N,Wn]=cheb1ord（WP,WS,Rp,Rs,'s'）;
[B,A]=cheby1（N,Rp,Wn,'s'）;[BT,AT]=lp2bp（B,A,W00,BW）;
[num,den]=bilinear（BT,AT,0.5）;
figure
（1）;[h,omega]=freqz（num,den,32）;
subplot（221）;stem（omega/pi,abs（h））;
xlabel（'\omega/\pi'）;ylabel（'|H（z）|'）;
subplot（222）;stem（omega/pi,20*log10（abs（h）））;
xlabel（'\omega/\pi'）;ylabel（'增益.dB'）;
x=filter（B,A,y）;X=fft（x,4096）;figure
（2）;
subplot（211）,plot（abs（Y））;axis（[0,4096,0,3]）;
title（'切比雪夫带通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4096,0,3]）;
title（'切比雪夫带通滤波后的信号频谱'）;sound（x,fs,bits）;
2.5窗函数设计低通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;fs=8000;
fc=1200;wc=2*fc/fs;wn=kaiser（49）;
b=fir1（48,wc,wn）;freqz（b,1）
figure;x=fftfilt（b,y）;X=fft（x,4096）;
subplot（211）,plot（abs（Y））;axis（[0,4000,0,3]）;
title（'窗函数低通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4000,0,3]）;
title（'窗函数低通滤波后的信号频谱'）;sound（x,fs,bits）;
2.6凯泽窗函数高通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;fs=22050;
fc=4800;wc=2*fc/fs;wn=kaiser（49）;
b=fir1（48,wc,'high',wn）;
freqz（b,1）figure;
x=fftfilt（b,y）;X=fft（x,4096）;
subplot（211）,plot（abs（Y））;axis（[0,4000,0,3]）;
title（'窗函数高通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4000,0,3]）;
title（'窗函数高通滤波后的信号频谱'）;sound（x,fs,bits）;
2.7凯泽窗函数带通滤波器滤波
程序：
clear;closeall;
[y,fs,bits]=wavread（'C:
\Users\winodws7\Desktop\数字信号处理实验\陈增贤录音3.wav'）;
sound（y,fs,bits）;Y=fft（y,4096）;fs=8000;
fc1=1000;fc2=3200;wc1=2*fc1/fs;wc2=2*fc2/fs;
wn=kaiser（49）;b=fir1（48,[wc1wc2],wn）;
freqz（b,1）figure;
x=fftfilt（b,y）;X=fft（x,4096）;
subplot（211）,plot（abs（Y））;axis（[0,4000,0,3]）;
title（'窗函数带通滤波前的信号频谱'）;
subplot（212）,plot（abs（X））;axis（[0,4000,0,3]）;
title（'窗函数带通滤波后的信号频谱'）;sound（x,fs,bits）;
2.8FIR和IIR滤波效果的比较
2.8.1低通滤波比较：
分析：
由图可知，窗函数的低通滤波后，信号在阻带的波动非常小，效果相对较好一些。
2.8.2高通滤波比较：
分析与总结：
由图对比可知，窗函数高通滤波后信号在阻带的波动较大，所以效果不如巴特沃斯滤波器。
但是巴特沃斯高通滤波后信号竟然全被滤除了，刚开始还以为是因为程序逻辑问题，后来才知道不是这样的，分析大致原因是我录的音频不在老师给的高通频率参数范围。
前面的低通滤波也出现了同样的现象，原因都是一样的。
由于时间有限，不能再重新录音分析了，但是还是能够做到了比较判断滤波器的优劣程度。
2.8.3带通滤波比较：
分析与总结：
对比两图可知，切比雪夫滤波频谱比较尖锐，只能通过很小频段的波形，而窗函数能通过适当范围的频率波形，且滤波效果较为理想。
4、实验六
1．如何获取信号
（1）建立如图所示仿真模型：
（2）Sinewave参数设置如下图所示：
（3）仿真图：
2.设信号
，试利用Simulink计算
点的离散傅立叶变换。
（1）建立仿真模型：
（2）设置sinefromworkspace的参数：
（3）仿真结果：
分析与总结：
（1）在Simulink中可产生连续信号和离散信号，正弦连续信号可用模块sinewave，设置参数即可调节相角和频率，以及幅值。
在产生离散正弦信号时，要使用DSPsinewave模块和buffer模块。
还可以通过命令窗口输入产生任意序列，这就需要用到signalfromworkspace模块。
（2）在simulink中可对离散信号进行离散傅立叶变换，其中要用到buffer模块和fft模块，可通过设置buffer模块的参数改变离散傅立叶变换的点数。
（3）个人觉得simulink仿真比编程要简单，方便一些，特别是在进行信号傅立叶变换的时候，编程容易产生不能很快被查找出来的错误。
五、心得体会
《数字信号处理》这门课程是沈老师带的我们最后一门课了，这门课结合了之前我们所学的《信号与系统》、《MATLAB教程》的知识，使我们的动手能力和思维创新能力得到了新的提升，当然，这也是最难的一个实验，所花费的时间比以前我们做的《自动控制原理》实验所花的时间要多得多，从算法的流程设计，程序编写，调试和修改，到最后的仿真都需要我们认真仔细，不可马虎，因为代码非常长，一不小心就出现error，实验课三次就上完了，可到最终大报告的完成我一共花了三周多吧。
可能是我太追求完美了，在语音信号处理那一块花的时间比较多，老是想用cooledit录一段非常适合做滤波处理的语音，最后的结果还算差强人意。
总的来说，这次的实验效果还是非常不错的，整个过程，有疑惑，抱怨，也有惊喜，快乐，所以这是一段很神奇的体验，老师，以后我还会继续努力的！
升级会员 