數(shù)字信號(hào)處理實(shí)驗(yàn)1-信號(hào)的頻譜分析與計(jì)算(附詳細(xì)程序)

PAGEPAGE16實(shí)驗(yàn)一信號(hào)的頻譜分析與計(jì)算一、用DFT進(jìn)行信號(hào)的譜分析1、用MATLABx(n

sin( 8

DFTN＝16，32程序：clc;clearallN1=16;N2=32;n1=0:N1-1;n2=0:N2-1;xn1=sin(pi/8*n1);xn2=sin(pi/8*n2);Xk1=fft(xn1,N1);Xk2=fft(xn2,N2);magXk1=abs(Xk1);magXk2=abs(Xk2);fd1=n1/N1;fd2=n2/N2;subplot(221);stem(n1,xn1);xlabel('n');ylabel('x(n)');title('x(n)N=16');subplot(223);stem(fd1,magXk1);xlabel('fd');ylabel('|X(k)|');title('|X(k)|N=16');gridsubplot(222);stem(n2,xn2);xlabel('n');ylabel('x(n)');title('x(n)N=32');subplot(224);stem(fd2,magXk2);xlabel('fd');ylabel('|X(k)|');title('|X(k)|N=32');grid圖像：N變了幅值。2、先用MATLAB產(chǎn)生出如下模擬信號(hào)：x(t)=cos(8πt)+cos(16πt)+cos(20πt)Fs=64HzDFTDFTN＝16，32，60DFT3個(gè)頻率，DFT長(zhǎng)度至少取多大？程序：clc;clearallFs=64;N1=16;%4,8,102fa=Fs/N,N32,N16時(shí)分離不出結(jié)果n1=0:N1-1;t1=n1/Fs;x1=cos(8*pi*t1)+cos(16*pi*t1)+cos(20*pi*t1);figure(1);subplot(211);plot(n1,x1);xlabel('n');ylabel('x(t)');title('N=16,x(t)');gridXk1=abs(fft(x1,N1));fa1=n1/N1*Fs;subplot(212);stem(fa1,Xk1);xlabel('f');ylabel('X(k)');title('N=16,X(k)');gridFs=64;N2=32;n2=0:N2-1;t2=n2/Fs;x2=cos(8*pi*t2)+cos(16*pi*t2)+cos(20*pi*t2);figure(2);subplot(211);plot(n2,x2);xlabel('n');ylabel('x(t)');title('N=32,x(t)');gridXk2=abs(fft(x2,N2));fa2=n2/N2*Fs;subplot(212);stem(fa2,Xk2);xlabel('f');ylabel('X(k)');title('N=32,X(k)');gridFs=64;N3=60;n3=0:N3-1;t3=n3/Fs;x3=cos(8*pi*t3)+cos(16*pi*t3)+cos(20*pi*t3);figure(3);subplot(211);plot(n3,x3);xlabel('n');ylabel('x(t)');title('N=60,x(t)');gridXk3=abs(fft(x3,N3));fa3=n3/N3*Fs;subplot(212);stem(fa3,Xk3);xlabel('f');ylabel('X(k)');title('N=60,X(k)');grid圖像：分析：由實(shí)驗(yàn)圖像可以看出當(dāng)N=32時(shí)，可以將信號(hào)包含的三個(gè)頻率分離表現(xiàn)出來。N=164Hz,8Hz,10Hz2Hz，頻率分辨率Δfa=Δfd*Fs=Fs/N2,N32,才能將三個(gè)頻率分量分離出來，N16時(shí)分離不出結(jié)果。N=60時(shí)，無法分別。模擬信號(hào)為x(t),抽樣后| ，截取N長(zhǎng)的x(n)，tnT Nx(n)x(n)WN

(n)X(ejw)1/2πX(ejw)WN

(ejw)W(ejw)

W(n)ejwnN1ejwn1ejwN

ejwNjwN/2ejwN

sin(wN/2)jw()N Nn

n0

1ejw ejwjw/2ejwsin(w/2)用matlab畫出幅譜圖clc;clearalln=(-100:0.1:100)N=30;x=abs((sin(0.005*pi*n*N))./(sin(0.005*pi*n)));plot(0.005*pi*n,x);axis([-0.550.55031]);grid可以看出：當(dāng)wN/2=k可以看出：當(dāng)wN/2=kπ，w=2πk/N時(shí)，|W(ejw)|0N當(dāng)w=0時(shí)，|W(ejw)|NNx(n)為余弦信號(hào)組合，X(ejw)是像1題中的取樣信號(hào)組成。X(ejw)1/2πX(ejw)WN

(ejw)相當(dāng)于對(duì)WN

(ejw)進(jìn)行頻譜搬移。XN(k)XN(ejw)|w2πk/Nk=fd*N=(4/64 8/64 10/64)*NN=32k=(245)N=60k=(3.757.59.735)非整數(shù)主瓣寬度Δw=4π/N,kw0w0+-2π/N，只采到了零w0附近的其他非零值。即產(chǎn)生了頻譜泄露。3、CZTx(t=cos(8πt+cos(16πtcos(20πt)的頻譜，設(shè)時(shí)域采樣頻率Fs=64HzN=64。CZT0~π/2M＝16，30打印信號(hào)的頻譜圖，并分析結(jié)果。程序：clc;clearallFs=64;N=64;n=0:N-1;t=n/Fs;x=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);subplot(311);plot(n,x);xlabel('n');ylabel('x(t)');title('N=64,x(t)');gridM1=16;M2=30;r1=0:M1-1;r2=0:M2-1;A=1;W1=exp(-1j*pi/2/M1);W2=exp(-1j*pi/2/M2);Xr2=abs(czt(x,M2,W2,A));fr1=r1/M1/4*Fs;fr2=r2/M2/4*Fs;subplot(312);stem(fr1,Xr1);xlabel('fr');ylabel('Xr');title('M=16,Xr');gridsubplot(313);stem(fr2,Xr2);xlabel('fr');ylabel('Xr');title('M=30,Xr');grid圖像： CZTMN頻率細(xì)化的目的。當(dāng)采樣范圍是0~π/2時(shí)，Wr=πr/2M=2πf fdr=r/4M Δfdr=1/4M采樣點(diǎn)數(shù)與頻率分辨率:Δfa=Fs/4M M=8,分辨率為2Hz,可以將三個(gè)信號(hào)分辨來。M=30時(shí)，同2題，采樣偏移產(chǎn)生了虛假諧波分量。二、DFTDCT的應(yīng)用FFT計(jì)算線性卷積（p9內(nèi)容）兩個(gè)序列：x1=[123456];x2=[3211234]，用DFT求兩者的循環(huán)卷積，F(xiàn)FTL=10，問哪些點(diǎn)上的循環(huán)卷積＝線性卷積？FFTL=12程序：clc;clearallx1=[123456];x2=[3211234];disp('Linearconvolution');y1=conv(x1,x2)subplot(321);stem(x1);subplot(323);stem(x2);subplot(322);stem(y1);gridtitle('Linearconvolution')disp('circularconvolutionusingFFTN=10');N=10;x1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];Xk1=fft(x1,N);Xk2=fft(x2,N);Yk=Xk1.*Xk2;y=ifft(Yk,N)subplot(324);stem(y);gridtitle('circularconvolutionusingFFTN=10');disp('circularconvolutionusingFFTN=12');N=12;x1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];Xk1=fft(x1,N);Xk2=fft(x2,N);Yk=Xk1.*Xk2;y=ifft(Yk,N)subplot(326);stem(y);gridtitle('circularconvolutionusingFFTN=12');結(jié)果：Linearconvolutiony1=3 8 14 21 30 42 37 36 40 43circularconvolutionusingFFTN=10y=41 32 14 21 30 42 37 36 40 43circularconvolutionusingFFTN=12y=Columns1through93.0000 8.0000 14.0000 21.0000 30.0000 42.000040.0000Columns10through1243.0000 38.0000 24.0000

38 2437.0000 36.0000圖像：Linearconvolution10 50500 2 4 64200 2 4 6 8

00 5 10 15circularconvolutionusingFFTN=105000 5 10circularconvolutionusingFFTN=125000 5 10 15分析：(1)M+N-1=12>L=10，會(huì)發(fā)生混疊，M+N-1=QQ-L<n<=L-10~9上2<n<9上的點(diǎn)的循環(huán)卷積等于線性卷積。(2) L=12M+N-1,不會(huì)發(fā)生混疊，結(jié)果與線性卷積相同。探地雷達(dá)信號(hào)譜分析附件一是采樣的一組探地雷達(dá)數(shù)據(jù)“radar-data.txt”，探地雷達(dá)探測(cè)的是水泥板的厚度，探地雷達(dá)接收的信號(hào)中含有水泥板上表面和下表面的反射波，采樣率為7000Hz，讀取這組數(shù)據(jù)，畫出其時(shí)域波形；DFT對(duì)這組數(shù)據(jù)進(jìn)行譜分析，畫出信號(hào)的幅譜圖，橫軸采用模擬頻率；讀取和顯示數(shù)據(jù)文件參考程序段如下：d=load('radar-data.txt');plot(d);grid程序：%MATLABProgramreadadata.md=load('radar-data.txt');%Plotusefuldatasegmentintimedomainsubplot(211),plot(d(550:680));gridtitle('radarsignalintimedomain')xlabel('Samplingfrequency=7KHz')ylabel('amplitude')Fs=7000;%ComputetheDFTofradarsignalD=abs(fft(d,2000));N=length(D);k=0:N/2-1;freq=k/N*Fs;subplot(212),plot(freq,D(1:N/2));title('ABSFFTofradardata')xlabel('frequency(Hz)')ylabel('amplitude')grid%MethodI%findthefirstMax.pointDmax1=0;Dkmax1=1;%========================================fork=1:N/2ifD(k)>Dmax1Dmax1=D(k);Dkmax1=k;endend%========================================%ComputetheFreqencyofthefirstMax.pointf1max=Dkmax1/N*Fs%findthefirstMin.pointafterthefirstMax.pointfork=Dkmax1:N/2ifD(k)<D(k+1)breakendend%findthesecondMax.pointafterthefirstMinpointDmax2=0;form=k:N/2ifD(m)>Dmax2Dmax2=D(m);Dkmax2=m;endend%ComputetheFreqencyofthesecondMax.pointf2max=Dkmax2/N*Fs%ComputetheFreq.differencebetweenthefirstandthesecondMax.pointsf21=f2max-f1max%MethodIIUsingMATLABfunctionMAX%findthefirstMax.pointDk=D(1:N/2);[D1,k1]=max(Dk);%ComputetheFreqencyofthefirstMax.pointf1=k1/N*Fs%findthefirstMin.pointafterthefirstMax.pointfork=k1:N/2ifDk(k)<Dk(k+1)breakendend%findthesecondMax.pointafterthefirstMinpointDk2=Dk(k+1:N/2);[D2,k2]=max(Dk2);%ComputetheFreqencyofthesecondMax.pointf2=(k+k2)/N*Fs%ComputetheFreq.differencebetweenthefirstandthesecondMax.pointsf=f2-f1圖像：edpma

10.50-0.5

radarsignalintimedomain-10 20 40 60 80 Samplingfrequency=7KHzABSFFTofradardata

120 14020edpma

1510500 500 1000 1500 2000frequency(Hz)

2500 3000 3500f1max=1638f2max=1806f21=168f1=1638f2=1806f=168DCT的圖象壓縮256*256“l(fā)ena.bmDCT（1）將Lena(256*256)分成8*8的塊；8*8DCT；MaskDCT系數(shù)去掉；IDCT恢復(fù)壓縮后的圖象；Mask（6）16*16（2）～（5）圖像壓縮參考程序如下：I=imread('lena256_256.bmp','bmp');Mask1=[11111110;11111100;11111000;11110000;11100000;11000000;10000000;00000000;];fori=1:32forj=1:32J1((i-1)*8+1:i*8,(j-1)*8+1:j*8)=idct2(dct2(I((i-1)*8+1:i*8,(j-1)*8+1:j*8)).*Mask1);endendimshow(I,[]);title('OriginalImage')；figure;imshow(J1,[]);title('compressImage’)程序：I=imread('Lena256.bmp','bmp');Mask1=[11111110;11111100;11111000;11110000;11100000;11000000;10000000;00000000;];fori=1:32forj=1:32J1((i-1)*8+1:i*8,(j-1)*8+1:j*8)=idct2(dct2(I((i-1)*8+1:i*8,(j-1)*8+1:j*8)).*Mask1);endendfigure(1);imshow(I,[]);title('OriginalImage');figure(2);imshow(J1,[]);title('compressImage');OriginalImage compressImage改變maskMask1=[11100000;11000000;10000000;00000000;00000000;00000000;00000000;00000000;];compressImage216*16分割時(shí)I=imread('Lena256.bmp','bmp');Mask1=[1111111111111110;1111111111111100;1111111111111000;