数字信号处理MATLAB上机实验2014年

1、1一、实验目的1、熟悉MATLAB的编程特点。2、了解DFT的计算及其应用。二、实验内容及要求1、用三种不同的DFT程序计算的傅里叶变换 ，并比较三种程序计算机运行时间。（1）编制用forloop语句的函数文件dft1.m，用循环变量逐点计算 。)()(8nRnx)(jeX)(kX上机实验一：离散傅里叶变换的计算和应用2（）编写用矩阵运算的的文件dft2.m ，完成下列矩阵运算：（）调用库函数，直接计算 。（）分别用上述三种不同方式编写的程序labdft.m:计算序列(自己构建)的傅里叶变换，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。) 1(.) 1 ()0(.) 1(.)

2、1 ()0(1)1()1(21012100000NxxxWWWWWWWWWWWWNXXXNNNNNNNNNNNNNNNNN)(kX)(nx)(jeX3%dft1.m: function Am,pha = dft1(x)N = length(x);w= exp(-j*2*pi/N); for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum);pha(k)=angle(sum);end 4%dft2.m: function Am,pha = dft2(x)N = length(x);n = 0:N-1; k =

3、 0:N-1; nk=n*k;w= exp(-j*2*pi/N); wnk = w.(nk); Xk = x*wnk; Am=abs(Xk);pha=angle(Xk);5%dft3.m: function Am,pha = dft3(x) Xk = fft(x); Am=abs(Xk);pha=angle(Xk);6%labdft.mx=ones(1,8),zeros(1,248);n=0:(length(x)-1);w=(2*pi/length(x)*n;7t=cputime;Am1,pha1=dft1(x);t1=cputime-tfigure(1)subplot(2,1,1);plot

4、(w,Am1,g);title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2);plot(w,pha1,r);grid;title(Phase Part);xlabel(frequency in radians);ylabel(arg(Xexp(jw)/radians);8t=cputime;Am2,pha2=dft2(x);t1=cputime-tfigure(2)subplot(2,1,1);plot(w,Am2,g);title(Magnitude part);xlabel(f

5、requency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2);plot(w,pha2,r);grid;title(Phase Part);xlabel(frequency in radians);ylabel(arg(Xexp(jw)/radians);9t=cputime;Am3,pha3=dft3(x);t1=cputime-tfigure(3)subplot(2,1,1);plot(w,Am3,g);title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|)

6、;subplot(2,1,2);plot(w,pha3,r);grid;title(Phase Part);xlabel(frequency in radians);ylabel(arg(Xexp(jw)/radians);1011、有一连续信号 以采样频率fs=3kHz对信号xa(t)进行采样，分析下列三种情况的幅频特性。（）采集数据长度点，做点的DFT。（）采集数据长度点，补零到点，做点的DFT。（）采集数据长度点，做点的DFT。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。改变采样频率和数据长度比较频谱分析结果，说明原因。3、提交以上上机频谱分析过程及结果图。以学号或姓名作

7、为word文件名。)100 . 92cos()100 . 72cos()105 . 62cos()(333ttttxa12%spectrum_com2% part 1% Spectrum based on the first 16 samples of x(n)fs=32000;figure(1)n=0:1:15;x=cos(2*pi/fs*6500*n)+cos(2*pi/fs*7000*n)+0.5*cos(2*pi/fs*9000*n);subplot(2,1,1);stem(n,x);title(signal x(n), 0 = n = 15);xlabel(n);axis(0,255

8、,-2.5,2.5);X=fft(x);magX=abs(X(1:1:16);k=0:1:15;f=fs/16*k;subplot(2,1,2);plot(f,magX);title(DTFT Magnitude);xlabel(frequency in Hz);axis(0,10000,0,25);13% part 2% High density spectrum (256 samples) based on the first 16 samples of x(n)figure(2)n=0:1:15;x=cos(2*pi/fs*6500*n)+cos(2*pi/fs*7000*n)+cos(

9、2*pi/fs*9000*n);y=x(1:1:16) zeros(1,240);k=0:1:255;subplot(2,1,1);stem(k,y);title(signal x(n),0 = n = 16+ 240 zeros);xlabel(n);axis(0,255,-2.5,2.5)X=fft(y);magX=abs(X(1:1:256);f=fs/256*k;subplot(2,1,2);plot(f,magX);title(DTFT Magnitude);xlabel(frequency in Hz);axis(0,10000,0,25);14% part 3% High res

10、olution spectrum based on 256 samples of the signal x(n)figure(3)n=0:1:255;x=cos(2*pi/fs*6500*n)+cos(2*pi/fs*7000*n)+cos(2*pi/fs*9000*n);subplot(2,1,1);stem(n,x);title(signal x(n), 0 = n 50dB% beta = 0.5842*(As-21)0.4+0.07886*(As-21) % kaiser;21dBAs50dBn=0:1:M-1;wc = (ws+wp)/275% Ideal LowPass filte

11、r computationalpha = (M-1)/2;m = n - alpha + eps;hd = sin(wc*m) ./ (pi*m);%wn=boxcar(M);%wn=triang(M);%wn=hanning(M);wn = hamming(M);%wn=blackman(M);% wn=kaiser(M,beta);% wn=chebwin(M+1,As); 76h = hd .* (wn);H,w = freqz(h,1,1000,whole);mag=abs(H);pha=angle(H);db=20*log10(mag+eps)/max(mag);delta_w =

12、2*pi/1000;Rp = -(min(db(1:1:wp/delta_w+1) % Passband RippleAs = -round(max(db(ws/delta_w+1:1:501) % Min Stopband attenuation77% plotsfigure(1);subplot(3,1,1); stem(n,hd); title(Ideal Impulse Response)axis(0 M-1 -0.1 0.3); xlabel(n); ylabel(hd(n);gridsubplot(3,1,2); stem(n,wn);title(Window)axis(0 M-1

13、 0 1.1); xlabel(n); ylabel(w(n);gridsubplot(3,1,3); stem(n,h);title(Actual Impulse Response)axis(0 M-1 -0.1 0.3); xlabel(n); ylabel(h(n);grid78figure(2);subplot(3,1,1); plot(w/pi,db);title(Actual Magnitude Response in dB);gridaxis(0 1 -100 10); xlabel(frequency in pi units); ylabel(|Hexp(jw) dB|)sub

14、plot(3,1,2);plot(w,mag);axis(0,3.15,0,1.5);title(Actual Magnitude Response);xlabel(w/rad); ylabel(|Hexp(jw)|); grid79subplot(3,1,3);plot(w,pha);hold on;n=0:7;x=zeros(8);plot(n,x,-); axis(0,3.15,-4,4);title(Actual Phase Response in radians);xlabel(w/rad); ylabel(arg(Hexp(jw); grid80n=0:256;fs=2000;x1

15、=sin(2*pi/fs*100*n);x2=sin(2*pi/fs*800*n);x=x1+x2;h=ifft(H);B=real(h)y=filter(h,1,x); % y = filter(B,1,x);B = fir1(M-1,wp,wn)81figure(3);subplot(2,1,1);plot(n,x); title(Souce signal:x=sin(2*pi*100*t)+sin(2*pi*800*t) );xlabel(time in s); ylabel(x); gridsubplot(2,1,2);plot(n,y);title(Filtered signal:y

16、);xlabel(time in s); ylabel(y); grid82838485% Frequency_Sample.m% Linear Phase Lowpass Filter Design % Digital Filter Specifications:% wp = 0.3pi,ws = 0.6pi, M =20% x(t)=sin(2*pi*100*t)+sin(2*pi*800*t),fs=2000;M = 20;wp = 0.3*pi;ws = 0.6*pi;alpha = -(M -1)/2;86kwp=fix(wp*M /(2*pi)+1kws=fix(ws*M /(2*

17、pi)+1k=0:1:M-1;wk = (2*pi/M)*k;Hrs=ones(1,kwp),zeros(1,M-2*kwp+1),ones(1,kwp-1)% Hrs = 1,1,1,1,zeros(1,13),1,1,1wd = 0,0.3,0.6,1;Hrd = 1,1,0,0; k1 = 0:floor(M-1)/2); k2 = floor(M-1)/2)+1:M-1;angH = -alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2);87H = Hrs.*exp(j*angH);h = real(ifft(H,M);Hm,wm = freqz(h,1,

18、512);mag=abs(Hm);pha=angle(Hm);db=20*log10(mag+eps)/max(mag)n=0:255;fs=2000;x1=sin(2*pi/fs*100*n);x2=sin(2*pi/fs*800*n);88x=x1+x2;y=filter(h,1,x);figure(1);subplot(2,2,1);plot(wk(1:11)/pi,Hrs(1:11),o,wd,Hrd);grid onaxis(0,1,-0.1,1.1); title(Frequency Samples: M=20)xlabel(frequency in pi units); ylab

19、el(Hr(k);set(gca,XTickMode,manual,XTick,0,0.3,0.6,1)set(gca,YTickMode,manual,YTick,0,1); 89subplot(2,2,2);stem(k,h);grid ontitle(Impulse Response); xlabel(n); ylabel(h(n);subplot(2,2,3);plot(wk(1:11)/pi,Hrs(1:11),o,wm/pi,mag); grid onaxis(0,1,-0.2,1.2); title(Amplitude Response)xlabel(frequency in p

20、i units); ylabel(Hr(w)set(gca,XTickMode,manual,XTick,0,0.3,0.6,1)set(gca,YTickMode,manual,YTick,0,1); 90subplot(2,2,4);plot(wm/pi,db);grid onaxis(0,1,-60,10);title(Magnitude Response); xlabel(frequency in pi units);ylabel(Decibels);set(gca,XTickMode,Manual,XTick,0;0.3;0.6;1);set(gca,YTickMode,Manual

21、,YTick,-18;0);set(gca,YTickLabelMode,manual,YTickLabels,18; 0)figure(2);subplot(2,1,1);plot(wm,mag);title(Designed Lowpass Filter Magnitude Response)xlabel(w/rad); ylabel(|Hexp(jw)|); grid91subplot(2,1,2);plot(wm,pha); title(Designed Lowpass Filter Phase Response in radians);xlabel(w/rad); ylabel(ar

22、g(Hexp(jw); grid%figure(3);subplot(2,1,1);plot(n,x); title(Souce signal:x=sin(2*pi*100*t)+sin(2*pi*800*t) );xlabel(time in s); ylabel(x); gridsubplot(2,1,2);plot(n,y);title(Filtered signal:y);xlabel(time in s); ylabel(y); grid92939495% Frequency_Sample1.m% Linear Phase Lowpass Filter Design % Digita

23、l Filter Specifications:% wp = 0.3pi,ws = 0.6pi,T=0.5,N=20% x(t)=sin(2*pi*100*t)+sin(2*pi*800*t),fs=2000;M = 20;wp = 0.3*pi;ws = 0.6*pi;alpha = -(M -1)/2;96kwp=fix(wp*M /(2*pi)+1kws=fix(ws*M /(2*pi)+1k=0:1:M-1;wk = (2*pi/M)*k;Hrs=ones(1,kwp),zeros(1,M-2*kwp+1),ones(1,kwp-1)% Hrs = 1,1,1,1,0.5,zeros(

24、1,11),0.5,1,1,1wd = 0,0.3,0.6,1;Hrd = 1,1,0,0; k1 = 0:floor(M-1)/2); k2 = floor(M-1)/2)+1:M-1;angH = -alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2);97H = Hrs.*exp(j*angH);h = real(ifft(H,M);Hm,wm = freqz(h,1,512);mag=abs(Hm);pha=angle(Hm);db=20*log10(mag+eps)/max(mag)n=0:255;fs=2000;x1=sin(2*pi/fs*100*n)

25、;x2=sin(2*pi/fs*800*n);x=x1+x2;98y=filter(h,1,x);figure(4);subplot(2,2,1);plot(wk(1:11)/pi,Hrs(1:11),o,wd,Hrd);grid onaxis(0,1,-0.1,1.1); title(Frequency Samples: M=20)xlabel(frequency in pi units); ylabel(Hr(k);set(gca,XTickMode,manual,XTick,0,0.3,0.4,0.6,1)set(gca,YTickMode,manual,YTick,0,0.5,1);

26、99subplot(2,2,2);stem(k,h);grid ontitle(Impulse Response); xlabel(n); ylabel(h(n);subplot(2,2,3);plot(wk(1:11)/pi,Hrs(1:11),o,wm/pi,mag); grid onaxis(0,1,-0.2,1.2); title(Amplitude Response)xlabel(frequency in pi units); ylabel(Hr(w)set(gca,XTickMode,manual,XTick,0,0.3,0.4,0.6,1)set(gca,YTickMode,ma

27、nual,YTick,0,0.5,1); 100subplot(2,2,4);plot(wm/pi,db);grid onaxis(0,1,-60,10);title(Magnitude Response); xlabel(frequency in pi units);ylabel(Decibels);set(gca,XTickMode,Manual,XTick,0;0.3;0.6;1);set(gca,YTickMode,Manual,YTick,-29;0);set(gca,YTickLabelMode,manual,YTickLabels,29; 0)%figure(5);subplot

28、(2,1,1);plot(wm,mag);101title(Designed Lowpass Filter Magnitude Response)xlabel(w/rad); ylabel(|Hexp(jw)|); gridsubplot(2,1,2);plot(wm,pha); title(Designed Lowpass Filter Phase Response in radians);xlabel(w/rad); ylabel(arg(Hexp(jw); grid%figure(6);subplot(2,1,1);plot(n,x); title(Souce signal:x=sin(

29、2*pi*100*t)+sin(2*pi*800*t) );xlabel(time in s); ylabel(x); gridsubplot(2,1,2);plot(n,y);title(Filtered signal:y);xlabel(time in s); ylabel(y); grid102103104105% Frequency_Sample2.m% Linear Phase Lowpass Filter Design % Digital Filter Specifications:% wp = 0.3pi,ws = 0.6pi,T=0.66,T=0.33,N=20% x(t)=s

30、in(2*pi*100*t)+sin(2*pi*800*t),fs=2000;M = 20;wp = 0.3*pi;ws = 0.6*pi;alpha = -(M -1)/2;kwp=fix(wp*M /(2*pi)+1kws=fix(ws*M /(2*pi)+1106k=0:1:M-1;wk = (2*pi/M)*k;Hrs=ones(1,kwp),zeros(1,M-2*kwp+1),ones(1,kwp-1)% Hrs = 1,1,1,1,0.66,0.33,zeros(1,9),0.33,0.66,1,1,1wd = 0,0.3,0.6,1;Hrd = 1,1,0,0; k1 = 0:

31、floor(M-1)/2); k2 = floor(M-1)/2)+1:M-1;angH = -alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2);H = Hrs.*exp(j*angH);h = real(ifft(H,M);107Hm,wm = freqz(h,1,512);mag=abs(Hm);pha=angle(Hm);db=20*log10(mag+eps)/max(mag)n=0:255;fs=2000;x1=sin(2*pi/fs*100*n);x2=sin(2*pi/fs*800*n);x=x1+x2;y=filter(h,1,x);108fig

32、ure(7);subplot(2,2,1);plot(wk(1:11)/pi,Hrs(1:11),o,wd,Hrd);grid onaxis(0,1,-0.1,1.1); title(Frequency Samples: M=20)xlabel(frequency in pi units); ylabel(Hr(k);set(gca,XTickMode,manual,XTick,0,0.3,0.4,0.5,0.6,1)set(gca,YTickMode,manual,YTick,0,0.33,0.66,1); subplot(2,2,2);stem(k,h);grid ontitle(Impu

33、lse Response); xlabel(n); ylabel(h(n);109subplot(2,2,3);plot(wk(1:11)/pi,Hrs(1:11),o,wm/pi,mag); grid onaxis(0,1,-0.2,1.2); title(Amplitude Response)xlabel(frequency in pi units); ylabel(Hr(w)set(gca,XTickMode,manual,XTick,0,0.3,0.4,0.5,0.6,1)set(gca,YTickMode,manual,YTick,0,0.33,0.66,1); subplot(2,

34、2,4);plot(wm/pi,db);grid on；axis(0,1,-60,10);title(Magnitude Response); xlabel(frequency in pi units);ylabel(Decibels);set(gca,XTickMode,Manual,XTick,0,0.3,0.4,0.5,0.6,1);set(gca,YTickMode,Manual,YTick,-31;0);set(gca,YTickLabelMode,manual,YTickLabels,31; 0)%110figure(8);subplot(2,1,1);plot(wm,mag);t

35、itle(Designed Lowpass Filter Magnitude Response)xlabel(w/rad); ylabel(|Hexp(jw)|); gridsubplot(2,1,2);plot(wm,pha); title(Designed Lowpass Filter Phase Response in radians);xlabel(w/rad); ylabel(arg(Hexp(jw); grid111%figure(9);subplot(2,1,1);plot(n,x); title(Souce signal:x=sin(2*pi*100*t)+sin(2*pi*8

36、00*t) );xlabel(time in s); ylabel(x); gridsubplot(2,1,2);plot(n,y);title(Filtered signal:y);xlabel(time in s); ylabel(y); grid112113114115一、实验目的 熟悉等波纹IIR、 FIR数字滤波器的设计方法。二、实验内容及要求 分别编制采用remez交换算法设计等波纹FIR数字滤波器和采用椭圆低通模拟滤波器设计等波纹IIR的程序。上机实验五：等波纹数字滤波器的设计116要求的低通FIR数字滤波器指标如下：通带截止频率 ，允许的幅度失真为 ，阻带边界频率为 ，最小衰减

37、为 RS=40dB，抽样频率为 80KHz，求出滤波器的单位取样响应，幅频和相频响应，绘出它们的图形。2、编制采用椭圆低通模拟滤波器设计低通IIR数字滤波器的程序。指标同上，比较滤波性能。3、提交以上上机结果图及分析过程。以学号或姓名作为word文件名。kHzfc1502. 0%,21kHzfs20117% Fir_remez.m% Lowpass FIR filter design using Parks McClellan remez algorithm% Digital Filter Specifications:% wp = 0.3pi, fc=1.5kHz; ws=0.6pi, de

38、lta1=0.0144, delta2=0.0101, Rp=0.25dB, As=40dB, Fs=2kHzfc=300;fs=600;Fs=2000;Rp=0.25;Rs=40;118f=fc,fs;m=1,0;delta1=(10(Rp/20)-1)/(10(Rp/20)+1)delta2=(1+delta1)*(10(-Rs/20)delta=delta1,delta2;%dev=0.0144,0.0101;N,f0,m0,w = remezord(f,m,delta,Fs)N=N+1; % modification119hn = remez(N,f0,m0,w);H,w=freqz(hn,1,512,Fs);%H=fft(hn,512);%w = 0:511; %w=2*w/512mag=abs(H);pha=angle(H);db=