数字信号处理实验1-信号的频谱分析与计算(附详细程序)

PAGEPAGE16实验一信号的频谱分析与计算一、用DFT进行信号的谱分析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产生出如下模拟信号：x(t)=cos(8πt)+cos(16πt)+cos(20πt)Fs=64HzDFTDFTN＝16，32，60DFT3个频率，DFT长度至少取多大？程序：clc;clearallFs=64;N1=16;%4,8,102fa=Fs/N,N32,N16时分离不出结果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图像：分析：由实验图像可以看出当N=32时，可以将信号包含的三个频率分离表现出来。N=164Hz,8Hz,10Hz2Hz，频率分辨率Δfa=Δfd*Fs=Fs/N2,N32,才能将三个频率分量分离出来，N16时分离不出结果。N=60时，无法分别。模拟信号为x(t),抽样后| ，截取N长的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可以看出：当wN/2=k可以看出：当wN/2=kπ，w=2πk/N时，|W(ejw)|0N当w=0时，|W(ejw)|NNx(n)为余弦信号组合，X(ejw)是像1题中的取样信号组成。X(ejw)1/2πX(ejw)WN

(ejw)相当于对WN

(ejw)进行频谱搬移。XN(k)XN(ejw)|w2πk/Nk=fd*N=(4/64 8/64 10/64)*NN=32k=(245)N=60k=(3.757.59.735)非整数主瓣宽度Δw=4π/N,kw0w0+-2π/N，只采到了零w0附近的其他非零值。即产生了频谱泄露。3、CZTx(t=cos(8πt+cos(16πtcos(20πt)的频谱，设时域采样频率Fs=64HzN=64。CZT0~π/2M＝16，30打印信号的频谱图，并分析结果。程序：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频率细化的目的。当采样范围是0~π/2时，Wr=πr/2M=2πf fdr=r/4M Δfdr=1/4M采样点数与频率分辨率:Δfa=Fs/4M M=8,分辨率为2Hz,可以将三个信号分辨来。M=30时，同2题，采样偏移产生了虚假谐波分量。二、DFTDCT的应用FFT计算线性卷积（p9内容）两个序列：x1=[123456];x2=[3211234]，用DFT求两者的循环卷积，FFTL=10，问哪些点上的循环卷积＝线性卷积？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');结果：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，会发生混叠，M+N-1=QQ-L<n<=L-10~9上2<n<9上的点的循环卷积等于线性卷积。(2) L=12M+N-1,不会发生混叠，结果与线性卷积相同。探地雷达信号谱分析附件一是采样的一组探地雷达数据“radar-data.txt”，探地雷达探测的是水泥板的厚度，探地雷达接收的信号中含有水泥板上表面和下表面的反射波，采样率为7000Hz，读取这组数据，画出其时域波形；DFT对这组数据进行谱分析，画出信号的幅谱图，横轴采用模拟频率；读取和显示数据文件参考程序段如下：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“lena.bmDCT（1）将Lena(256*256)分成8*8的块；8*8DCT；MaskDCT系数去掉；IDCT恢复压缩后的图象；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分割时I=imread('Lena256.bmp','bmp');Mask1=[1111111111111110;1111111111111100;1111111111111000;