功率谱分析由题目内容，设采样频率fs=1000HZ，数据长度为256，模型阶数为14，f1=200,f2=300、250。

（1）用最大熵法进行谱估计运行程序后，观察图像f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，这是因为AR模型的谱估计隐含着对数据和自相关函数的外推，其长度可能会超过给定长度，分辨率不受信源信号的限制。

（2）分别用Levinson递推法和Burg法进行功率谱分析①Levinson递推法运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，但本题中信号为正弦信号加白噪声，故图像观察不明显。

②Burg法运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率。

(3)改变信号的相位、频率、信噪比，上述谱分析结果有何变化如果正弦信号的频率过大，超过fs/2，会产生频率混叠现象，输入f1=600HZ，会在400HZ处产生一个波峰；降低信噪比会导致谱分辨率下降；信号起始相位的变动可导致谱线的偏移和分裂（我的图像观察不到）。

最大熵法估计N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');title('MEM f2/fs=0.3,Nfft=256,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*250*t); %0.25xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.25,Nfft=256,Oder=14');gridN=1024;Nfft=512; %修改数据长度512Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=512,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,24,Nfft,Fs); %修改阶数为24subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=24');GridBurg法估计N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*250*t); %0.25xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=250,Nfft=256, Oder=14');gridN=1024;Nfft=512; %修改数据长度512Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=512, Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,24,Nfft,Fs); %修改阶数为24subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=24');gridLevinson递推法N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*250*t); %0.25xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.25,Oder=14');gridN=1024;Nfft=512; %修改数据长度512Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=512,f2/fs=0.3,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,24,Nfft,Fs); %修改阶数为24subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=24');grid最大熵法改变信号的相位、频率、信噪比N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t+pi/6); %相位加了pi/6x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14，相位加pi/6');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比改为5[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14，性噪比=5');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*300*t);x2=sin(2*pi*400*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pmem(xn,14,Nfft,Fs);subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f1/fs=0.3,f2/fs=0.4,Nfft=256,Oder=14');gridBurg改变信号的相位、频率、信噪比N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t+pi/6); %相位加了pi/6x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14，相位加pi/6'); gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*300*t);x2=sin(2*pi*400*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f1/fs=300,f2/fs=400,Nfft=256, Oder=14'); gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比改为5[Pxx1,f]=pburg(xn,14,Nfft,Fs);subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');title('Burg f2/fs=300,Nfft=256, Oder=14,性噪比=5');gridLevinson法改变信号的相位、频率、信噪比N=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);subplot(4,1,1)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14');gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t+pi/6); %相位加了pi/6x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);subplot(4,1,2)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14，相位加pi/6'); gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*300*t);x2=sin(2*pi*400*t); %0.3xn=x1+awgn(x1,10)+x2+awgn(x2,10);[Pxx1,f]=pyulear(xn,14,Nfft,Fs);%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);subplot(4,1,3)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f1/fs=0.3,f2/fs=0.4,Oder=14'); gridN=1024;Nfft=256;Fs=1000;n=0:N-1;t=n/Fs;x1=sin(2*pi*200*t);x2=sin(2*pi*300*t); %0.3xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比位5[Pxx1,f]=pyulear(xn,14,Nfft,Fs);%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);subplot(4,1,4)plot(f,10*log10(Pxx1));xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14,性噪比=5'); grid。