數(shù)字信號(hào)處理MATLAB上機(jī)實(shí)驗(yàn)2014年

1、1一、實(shí)驗(yàn)?zāi)康?、熟悉MATLAB的編程特點(diǎn)。2、了解DFT的計(jì)算及其應(yīng)用。二、實(shí)驗(yàn)內(nèi)容及要求1、用三種不同的DFT程序計(jì)算的傅里葉變換 ，并比較三種程序計(jì)算機(jī)運(yùn)行時(shí)間。（1）編制用forloop語句的函數(shù)文件dft1.m，用循環(huán)變量逐點(diǎn)計(jì)算 。)()(8nRnx)(jeX)(kX上機(jī)實(shí)驗(yàn)一：離散傅里葉變換的計(jì)算和應(yīng)用2（）編寫用矩陣運(yùn)算的的文件dft2.m ，完成下列矩陣運(yùn)算：（）調(diào)用庫函數(shù)，直接計(jì)算 。（）分別用上述三種不同方式編寫的程序labdft.m:計(jì)算序列(自己構(gòu)建)的傅里葉變換，并畫出相應(yīng)的幅頻和相頻特性，再比較各個(gè)程序的計(jì)算機(jī)運(yùn)行時(shí)間。) 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、有一連續(xù)信號(hào) 以采樣頻率fs=3kHz對(duì)信號(hào)xa(t)進(jìn)行采樣，分析下列三種情況的幅頻特性。（）采集數(shù)據(jù)長度點(diǎn)，做點(diǎn)的DFT。（）采集數(shù)據(jù)長度點(diǎn)，補(bǔ)零到點(diǎn)，做點(diǎn)的DFT。（）采集數(shù)據(jù)長度點(diǎn)，做點(diǎn)的DFT。觀察三幅不同頻率特性圖，分析和比較它們的特點(diǎn)以及形成的原因。改變采樣頻率和數(shù)據(jù)長度比較頻譜分析結(jié)果，說明原因。3、提交以上上機(jī)頻譜分析過程及結(jié)果圖。以學(xué)號(hào)或姓名作

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一、實(shí)驗(yàn)?zāi)康?熟悉等波紋IIR、 FIR數(shù)字濾波器的設(shè)計(jì)方法。二、實(shí)驗(yàn)內(nèi)容及要求 分別編制采用remez交換算法設(shè)計(jì)等波紋FIR數(shù)字濾波器和采用橢圓低通模擬濾波器設(shè)計(jì)等波紋IIR的程序。上機(jī)實(shí)驗(yàn)五：等波紋數(shù)字濾波器的設(shè)計(jì)116要求的低通FIR數(shù)字濾波器指標(biāo)如下：通帶截止頻率 ，允許的幅度失真為 ，阻帶邊界頻率為 ，最小衰減

37、為 RS=40dB，抽樣頻率為 80KHz，求出濾波器的單位取樣響應(yīng)，幅頻和相頻響應(yīng)，繪出它們的圖形。2、編制采用橢圓低通模擬濾波器設(shè)計(jì)低通IIR數(shù)字濾波器的程序。指標(biāo)同上，比較濾波性能。3、提交以上上機(jī)結(jié)果圖及分析過程。以學(xué)號(hào)或姓名作為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=