《功率谱估计实验》由会员分享，可在线阅读，更多相关《功率谱估计实验（6页珍藏版）》请在金锄头文库上搜索。
1、Power Spectrum Estimation 上机实验姓名：戚永前 学号：10720885一：设输入音频信号 ，取 f=2KHz，f s=40KHz，N =128，W( n)为三角窗)2cos()fttXa序列。用计算机求出功率谱估值 及不分段的 。 kPx )(kBwx实验中使用 matlab 语言进行编程，在 MATLAB R2008b 上调试通过。所得到的结果如图 1，具体源代码参见附录程序源代码 1。图 1 实验 1 结果图上图左下图所示的是 所对应到实际频率的功率谱，右下图是不分段的)(kPx所对应到实际频率的功率谱。图中我们可以清楚的看到，在 2kHz 频率处由一个尖)(kB
2、wx峰，这与我们测试的音频信号的频率正好是吻合的，证明实验结果是正确的。二：设输入音频信号 ，取 f=2KHz，f s=40KHz，N =128，W( n)为三角窗)2cos()fttXa序列。用计算机求分段的 。用 2:1 覆盖分段，设各段的长度为 M=32。 kBwx具体测量流程图如图 2 所示：取样fs=20KHzXa(t) 2:1 分段M=32三角窗加权 规范运算 -/()-(2)jnMunxjeM/2 点ODFT取模平方再除 MU补奇数谱线内存 UBxx(k) L 个分段求平均 D(2k)图 2 2:1 分段覆盖测量流程图计算机编程流程图如图 3 所示。2:1 分原信号为 L 段，M
3、=32 ，三角窗加权每段 ，()()iixkWk并将其重组为矩阵 xn2逐段规范运算 -/()-(2)iiijnMunxje取实部、虚部，用 ODFT 得偶数谱线并求模平方除 N由 ODFT 谱线共轭奇偶对称性补奇数谱线Pxx(k)输入 N 及 Ts开始结束图 3 计算机测量流程图实验结果如图 4 所示,具体源程序参见附录程序源代码 2。图 4 2:1 分段覆盖所得的功率谱由图我们可以看到分段所得的功率谱分辨率较低，这与理论结果是一致的。附录：程序源代码 1(matlab)：%filename power_spectrum_estimation.m%version:10282009 pfx%t
4、his matlab program is used for power spectrum estimation. clear;clf;Fs=40000;%sample frequency 40kHz;N=128;%the number of samplesFc=2000;%the frequency of audio signal 2kHz;n=0:N-1;t=n/Fs;%time domainXt=cos(2*pi*Fc*t);%the audio signalTRGt=triang(N);%triangular windowU=sum(TRGt.2)/N;% the energy of
5、the windowYt=Xt.*TRGt;% the windowed signalfigure(1);subplot(2,2,1);plot(t,Xt,k);%display the original audio signal;title(Original Audio signal , Fs=40kHz,N=128,Fc=2kHz);xlabel(t/s);ylabel(X(t)/V);subplot(2,2,2);plot(t,Yt,k);%display the windowed audio signal;title(Windowed Audio signal, N=128,Trian
6、gular Window);xlabel(t/s);ylabel(Y(t)/V);for m=0:N/2-1%normalization ODFTUt(m+1)=(Xt(m+1)-j*Xt(N/2+m+1)*exp(-j*m*pi/N);Vt(m+1)=(Yt(m+1)-j*Yt(N/2+m+1)*exp(-j*m*pi/N);endDk1=fft(Ut);%N/2 points fft ODFTDk2=fft(Vt);%N/2 points fft ODFTPs1=abs(Dk1).2;%the power of the magnitude Ps1=Ps1/N;Ps2=abs(Dk2).2;
7、%the power of the magnitude Ps2=Ps2/N/U;for m=0:N/2-1pow1(2*m+1)=Ps1(m+1);%the even componentspow1(N-2*m)=Ps1(m+1);%add the odd componentspow2(2*m+1)=Ps2(m+1);%the even componentspow2(N-2*m)=Ps2(m+1);%add the odd componentsendpxx1=10*log10(pow1);%display as dBpxx2=10*log10(pow2);%display as dBf=n/N*
8、Fs;%the coordinate frequency subplot(2,2,3);plot(f,pxx1,k);%display the power spectrum of the original signaltitle(Original Signal Power Spectrum Pxx(k) N=128);xlabel(frequency/Hz);ylabel(Power Spectrum(dB);subplot(2,2,4);plot(f,pxx2,k);%display the power spectrum of the windowed signaltitle(Windowe
9、d Signal Power Spectrum Bxx(k) N=128);xlabel(frequency/Hz);ylabel(Power Spectrum(dB);程序源代码 2(matlab)：%filename: power_spectrum_estimation_2nd.m%version: 10302009 pfx%this matlab program is used for power spectrum estimation. clear;clf;Fs=40000;%sample frequency 40kHz;N=128;%the number of samplesM=32
10、;L=7;Fc=2000;%the frequency of audio signal 2kHz;n=0:N-1;t=n/Fs;% time domainXt=cos(2*pi*Fc*t);%the audio signalTRGt=triang(M);%triangular windowU=sum(TRGt.2)/M;% the energy of the windowXt1=Xt(1:32);Xt2=Xt(17:48);Xt3=Xt(33:64);Xt4=Xt(49:80);Xt5=Xt(65:96);Xt6=Xt(81:112);Xt7=Xt(97:128);Yt1=Xt1.*TRGt;
11、% the windowed signal1Yt2=Xt2.*TRGt;% the windowed signal2Yt3=Xt3.*TRGt;% the windowed signal3Yt4=Xt4.*TRGt;% the windowed signal4Yt5=Xt5.*TRGt;% the windowed signal5Yt6=Xt6.*TRGt;% the windowed signal6Yt7=Xt7.*TRGt;% the windowed signal7for m=0:M/2-1%normalization ODFTUt1(m+1)=(Yt1(m+1)-j*Yt1(M/2+m
12、+1)*exp(-j*m*pi/M); Ut2(m+1)=(Yt2(m+1)-j*Yt2(M/2+m+1)*exp(-j*m*pi/M);Ut3(m+1)=(Yt3(m+1)-j*Yt3(M/2+m+1)*exp(-j*m*pi/M);Ut4(m+1)=(Yt4(m+1)-j*Yt4(M/2+m+1)*exp(-j*m*pi/M);Ut5(m+1)=(Yt5(m+1)-j*Yt5(M/2+m+1)*exp(-j*m*pi/M);Ut6(m+1)=(Yt6(m+1)-j*Yt6(M/2+m+1)*exp(-j*m*pi/M);Ut7(m+1)=(Yt7(m+1)-j*Yt7(M/2+m+1)*e
13、xp(-j*m*pi/M);endDk1=fft(Ut1);%M/2 points fft ODFTDk2=fft(Ut2);%M/2 points fft ODFTDk3=fft(Ut3);%M/2 points fft ODFTDk4=fft(Ut4);%M/2 points fft ODFTDk5=fft(Ut5);%M/2 points fft ODFTDk6=fft(Ut6);%M/2 points fft ODFTDk7=fft(Ut7);%M/2 points fft ODFTPs1=abs(Dk1).2/M/U;%the power of the magnitude Ps2=a
14、bs(Dk2).2/M/U;%the power of the magnitude Ps3=abs(Dk3).2/M/U;%the power of the magnitude Ps4=abs(Dk4).2/M/U;%the power of the magnitude Ps5=abs(Dk5).2/M/U;%the power of the magnitude Ps6=abs(Dk6).2/M/U;%the power of the magnitude Ps7=abs(Dk7).2/M/U;%the power of the magnitude for m=0:M/2-1pow1(2*m+1
15、)=Ps1(m+1);%the even componentspow1(M-2*m)=Ps1(m+1);%add the odd componentspow2(2*m+1)=Ps2(m+1);%the even componentspow2(M-2*m)=Ps2(m+1);%add the odd componentspow3(2*m+1)=Ps3(m+1);%the even componentspow3(M-2*m)=Ps3(m+1);%add the odd componentspow4(2*m+1)=Ps4(m+1);%the even componentspow4(M-2*m)=Ps
16、4(m+1);%add the odd componentspow5(2*m+1)=Ps5(m+1);%the even componentspow5(M-2*m)=Ps5(m+1);%add the odd componentspow6(2*m+1)=Ps6(m+1);%the even componentspow6(M-2*m)=Ps6(m+1);%add the odd componentspow7(2*m+1)=Ps7(m+1);%the even componentspow7(M-2*m)=Ps7(m+1);%add the odd componentsendPOW=(pow1+pow2+pow3+pow4+pow5+pow6+pow7)/L;pxx=10*log10(POW);%display as dBf=0:M-1/M*F
