数字信号处理实验报告二与三

实验二用FFT进行谱分析一.实验目的：1进一步加深DFT算法原理和根本性质的理解〔因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的根本性质〕。熟悉FFT程序结构及编程方法。2熟悉应用FFT对确定信号进行谱分析方法，熟悉FFT算法原理和FFT子程序的应用。3学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应有FFT。二.实验内容：用matlab编程产生并画出信号x1(n)、x2(n)、x3(n)、x4(n)、x5(n)。图1图1用matlab编制FFT函数对上述信号进行频谱分析，并画出上述信号谱图。三.实验结果〔1〕1.%Thisprogrammistogeneratesignalx1(n)=R4(n).k=-6:6;x=[zeros(1,6),ones(1,4),zeros(1,3)];stem(k,x);〔信号图如图1〕title('图1');2.n=-5:1:10;x=(n+1).*(n>=0&n<=3)+(8-n).*(n>=4&n<=7)+0;stem(n,x)；title('图2');3.n=-5:10;x=(4-n).*(n>=0&n<=3)+(n-3).*(n>=4&n<=7);stem(n,x);title('图3');4.n=-10:10;x=cos(pi/4*n);stem(n,x);title('图4');5.n=-10:10;x=sin(pi/8*n);stem(n,x);title('图5');实验结果〔2〕:FFT算法functiony=myditfft(x)%y=myditfft(x)%本程序对输入序列x实现DIT-FFT基2算法，点数取大于等于x长度的2的幂次%x为给定时间序列%y为x的离散傅立叶变换m=nextpow2(x);N=2^m;%求x的长度对应的2的最低幂次miflength(x)<N;%假设x的长度不是2的幂，补零到2的整数幂x=[x,zeros(1,N-length(x))];endnxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;%求1:2^m数列的倒序y=x(nxd);%将x倒序排列作为y的初始值formm=1:m;%将DFT作m次基2分解,从左到右，对每次分解作DFT运算Nmr=2^mm;u=1;%旋转因子u初始化为WN^0=1WN=exp(-i*2*pi/Nmr);%本次分解的根本DFT因子WN=exp(-i*2*pi/Nmr)forj=1:Nmr/2;%本次跨越间隔内的各次蝶形运算fork=j:Nmr:N;%本次蝶形运算的跨越间隔为Nmr=2^mmkp=k+Nmr/2;%确定蝶形运算的对应单元下标t=y(kp)*u;%蝶形运算的乘积项y(kp)=y(k)-t;%蝶形运算y(k)=y(k)+t;%蝶形运算endu=u*WN;%修改旋转因子,多乘一个根本DFT因子WNendend1.k=-6:6;x=[zeros(1,6),ones(1,4),zeros(1,3)];y=myditfft(x);k=-6:9;stem(k,y);xlabel('m');ylabel('X[M]');title('FFT图');2.n=-5:1:10;x=(n+1).*(n>=0&n<=3)+(8-n).*(n>=4&n<=7)+0;y=myditfft(x);stem(n,y);xlabel('n');ylabel('X[M]');title('FFT图');3.n=-5:10;x=(4-n).*(n>=0&n<=3)+(n-3).*(n>=4&n<=7);y=myditfft(x);stem(n,y);xlabel('n');ylabel('X[M]');title('FFT3');4.n=-10:10;x=cos(pi/4*n);y=myditfft(x);n=-10:21;stem(n,y);xlabel('n');ylabel('X[M]');title('FFT4');5.n=-10:10;x=sin(pi/8*n);y=myditfft(x);n=-10:21;stem(n,y);xlabel('n');ylabel('X[M]');title('FFT5');四.简要答复以下问题：①在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？答：不相同。当N=8时，序列x1(n)和x2(n)中相同的元素值对应的n值是不同的，所乘的旋转因子的值也不同，因而得到的最终结果也是不同的。同理，N=16时，所得的幅频特性也是不同的。②FFT在什么条件下也可以用来分析周期信号序列的频谱？如果正弦信号系统sin(2πf0k),f0=0.1Hz，用16点FFT来做DFT运算，得到的频谱是信号本身的真实谱吗？为什么？答：由于FFT算法对序列长度的要求是N=2^M，M为正整数。所以，当周期信号序列一个周期的长度满足N=2^M〔M为正整数〕的条件时，FFT可以用来分析周期信号的频谱。不是真实的频谱。因为序列的周期N=10不是2的整数次幂，所以不是真实的。实验三：用双线性变换法设计IIR数字滤波器一.实验目的：1.掌握用双线性变换法设计IIRDF的原理及具体设计方法，熟悉用双线性变换法设计IIRDF的计算机编程。2.观察用双线性变换法设计的DF的频响特性，了解双线性变换法的特点。3.熟悉用双线性变换法设计BW和CB型DF的全过程。二，实验原理：为了克服冲激响应不变法产生的频率混叠现象，这是从S平面到Z平面的标准变换z＝esT的多值对应关系导致的，为了克服这一缺点，产生了双线性变换法。下式就是模拟滤波器和经采样后的数字滤波器之间的变换关系：三.实验内容：1.读懂所给参考程序，熟悉程序的整体结构和功能。2.设计一个CBI型低通DF，通带截频fp=3000Hz，衰耗满足Apmax=3dB，阻带截频fS=3400Hz,衰耗ASmin=31dB,取样频率fs=8000Hz。画出其模拟滤波器及数字滤波器频响特性曲线。程序：FS=8000;Fl=3000;Fh=3400;%通带、阻带截止频率Rp=3;Rs=31;wp=Fl*2*pi;%临界频率采用角频率表示ws=Fh*2*pi;%临界频率采用角频率表示wp1=wp/FS;%求数字频率ws1=ws/FS;%求数字频率OmegaP=2*FS*tan(wp1/2);%频率预畸OmegaS=2*FS*tan(ws1/2);%频率预畸%选择滤波器的最小阶数[n,Wn]=cheb1ord(OmegaP,OmegaS,Rp,Rs,'s');%此处是代入经预畸变后获得的归一化模拟频率参数[bt,at]=cheby1(n,Rp,Wn,'s');%设计一个n阶的契比雪夫I型模拟滤波器subplot(2,1,1);[h,w]=freqz(bt,at);plot(w*FS/(2*pi),abs(h));grid;xlabel('频率/Hz');ylabel('幅值');title('模拟滤波器频响特性');[bz,az]=bilinear(bt,at,FS);%双线性变换为数字滤波器subplot(2,1,2);[H,W]=freqz(bz,az);%求解数字滤波器的频率响应plot(W*FS/(2*pi),abs(H));grid;xlabel('频率/Hz');ylabel('幅值');title('数字滤波器频响特性');3、设计一个BW型低通DF，满足：通带截频fp=100Hz，衰耗满足Apmax=3dB，阻带截频fS=400Hz,衰耗ASmin=15dB,取样频率fs=2000Hz。画出其模拟滤波器及数字滤波器频响特性曲线。FS=2000;Fl=100;Fh=400;%通带、阻带截止频率Rp=3;Rs=15;wp=Fl*2*pi;%临界频率采用角频率表示ws=Fh*2*pi;%临界频率采用角频率表示wp1=wp/FS;%求数字频率ws1=ws/FS;%求数字频率OmegaP=2*FS*tan(wp1/2);%频率预畸OmegaS=2*FS*tan(ws1/2);%频率预畸%选择滤波器的最小阶数[n,Wn]=buttord(OmegaP,OmegaS,Rp,Rs,'s');%此处是代入经预畸变后获得的归一化模拟频率参数[bt,at]=butter(n,Wn,'s');%设计一个n阶的契比雪夫I型模拟滤波器subplot(2,1,1);[h,w]=freqz(bt,at);plot(w*FS/(2*pi),abs(h));grid;xlabel('频率/Hz');ylabel('幅值');title('模拟滤波器频响特性');[bz,az]=bilinear(bt,at,FS);%双线性变换为数字滤波器[H,W]=freqz(bz,az);%求解数字滤波器的频率响应subplot(2,1,2);plot(W*FS/(2*pi),abs(H));grid;xlabel('频率/Hz');ylabel('幅值');title('数字滤波器频响特性');设计一个BW型高通DF，满足：通带截频fp=400Hz，衰耗满足Apmax=3dB，阻带截频fS=350Hz,衰耗ASmin=15dB,取样频率fs=1000Hz。画出其模拟滤波器及数字滤波器频响特性曲线。fp=400;fs=350;ft=1000;rp=3;rs=15;wp=fp/(ft/2);ws=fs/(ft/2);%利用Nyquist频率进行归一化[n,wc]=buttord(wp,ws,rp,rs,'s');%求模拟滤波器的最小阶数和截止频率[b,a]=butter(n,wc,'high','s');%设计高通模拟滤波器系数b，a%Omega1=linspace(0,1000,5);%[h,w]=freqs(b,a,Omega1);[h,w]=freqs(b,a);subplot(2,1,1);plot(w*ft/(2*pi),abs(h));grid;title('模拟滤波器');xlabel('频率/Hz');ylabel('幅值');[N,Wc]=buttord(wp,ws,rp,rs);%求数字滤波器的最小阶数和截止频率[B,A]=butter(n,wc,'high');%设计高通数字滤波器系数b，a[H,W]=freqz(B,A);%绘出频率响应曲线subplot(2,1,2);plot(W*ft/(2*pi),abs(H));grid;title('数字滤波器');xlabel('频率/Hz');ylabel('幅值');设计一个CB型带通DF，满足：通带边界频率为100Hz～500Hz，通带衰耗小于3dB，过渡带宽20Hz，阻带衰耗大于15dB，取样频率fs=2000Hz。画出其模拟滤波器及数字滤波器频响特性曲线。fp=[100,500];fs=[80,520];ft=2000;rp=3;rs=15;wp=fp/(ft/2);ws=fs/(ft/2);[n,wn]=cheb2ord(wp,ws,rp,rs,'s');[b,a]=cheby2(n,rs,wn,'s');[h,w]=freqs(b,a);subplot(2,1,1);plot(w*ft/(2*pi),abs(h));grid;title('模拟滤波器');xlabel('频率/Hz');ylabel('幅值');[N,Wn]=cheb2ord(wp,ws,rp,rs);[B,A]=cheby2(N,r