MATLAB经典功率谱估计法

一、作业内容：对两个正弦信号做叠加后，计算离散随机过程信号的功率谱函数，由功率谱，估计信号的频率。在matlab上实现之，并观察波形进行验证。二、实现步骤：（一）、构造环境：1、两个正弦波分别为A*sin(2*pi*f1*n+a)、B*sin(2*pi*f2*n+a)，规定取样点范围n=1~128；构造函数x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a)；2、在x1基础上加入加性高斯白噪声，取定信噪比为+3,来定义x2的函数为x2=x1+W（噪声）；3、对离散信号x2做非参数化谱估计，以傅里叶变换为基础，先对x2做傅里叶变换，求出其频谱；4、求x2的功率谱p(w)，用周期图法；用间接法；分别估计做出功率谱，并输出其功率谱波形。5、更改采样点数，验证功率谱波形的主瓣函数图形什么情况下有重叠程度、什么情况下能够很好的区分开来。（二）、在matlab中编写相应程序：clearall;%清除工作空间所有之前的变量closeall;%关闭之前的所有的figureclc;%清除命令行之前所有的文字n=1:1:128;%设定采样点n=1-128f1=0.2;%设定f1频率的值0.2f2=0.213;%设定f2频率的值0.213A=1;%取定第一个正弦函数的振幅B=1;%取定第一个正弦函数的振幅a=0;%设定相位为0x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a);%定义x1函数，不添加高斯白噪声x2=awgn(x1,3);%在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0;%定义临时值，并规定初始值为0temp=fft(x2,128);%对x2做快速傅里叶变换pw1=abs(temp).^2/128;%对temp做经典功率估计k=0:length(temp)-1;w=2*pi*k/128;figure(1);%输出x1函数图像plot(w/pi/2,pw1)%输出功率谱函数pw1图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1添加高斯白噪声后的,周期图法功率频谱分析');grid;%-------------------------------------------------------------------------pw2=temp.*conj(temp)/128;%对temp做向量的共轭乘积k=0:length(temp)-1;w=2*pi*k/128;figure(2);plot(w/pi/2,pw2);%输出功率谱函数pw2图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1自相关法功率谱估计');grid;三、在matlab中，输出的功率谱图像。1、用直接法，功率谱图像。2、用间接法（自相关函数）做出功率谱图像。由上面的图像发现，峰值主瓣有重合，不能区分开来，峰值点处的坐标只有一个峰值点x=0.2109。所以，不能有效的估计出f1，f2。3、更改采样点n=512后。clearall;%清除工作空间所有之前的变量closeall;%关闭之前的所有的figureclc;%清除命令行之前所有的文字n=1:1:512;%设定采样点n=1-128f1=0.2;%设定f1频率的值0.2f2=0.213;%设定f2频率的值0.213A=1;%取定第一个正弦函数的振幅B=1;%取定第一个正弦函数的振幅a=0;%设定相位为0x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a);%定义x1函数，不添加高斯白噪声x2=awgn(x1,3);%在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0;%定义临时值，并规定初始值为0temp=fft(x2,512);%对x2做快速傅里叶变换pw1=abs(temp).^2/512;%对temp做经典功率估计k=0:length(temp)-1;w=2*pi*k/512;figure(1);%输出x1函数图像plot(w/pi/2,pw1)%输出功率谱函数pw1图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1添加高斯白噪声后的,周期图法功率频谱分析');grid;%-------------------------------------------------------------------------pw2=temp.*conj(temp)/512;%对temp做向量的共轭乘积k=0:length(temp)-1;w=2*pi*k/512;figure(2);plot(w/pi/2,pw2);%输出功率谱函数pw2图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1自相关法功率谱估计');grid;（1）、采点为n=512时，功率谱图像由图像坐标可以清楚的得到峰值点的坐标，X1=0.1992，X2=0.2129，这样不会出现峰值的重合现象，能够清楚的估计信号的频率。4、更改采样点n=1024后。clearall;%清除工作空间所有之前的变量closeall;%关闭之前的所有的figureclc;%清除命令行之前所有的文字n=1:1:1024;%设定采样点n=1-128f1=0.2;%设定f1频率的值0.2f2=0.213;%设定f2频率的值0.213A=1;%取定第一个正弦函数的振幅B=1;%取定第一个正弦函数的振幅a=0;%设定相位为0x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a);%定义x1函数，不添加高斯白噪声x2=awgn(x1,3);%在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0;%定义临时值，并规定初始值为0temp=fft(x2,1024);%对x2做快速傅里叶变换pw1=abs(temp).^2/1024;%对temp做经典功率估计k=0:length(temp)-1;w=2*pi*k/1024;figure(1);%输出x1函数图像plot(w/pi/2,pw1)%输出功率谱函数pw1图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1添加高斯白噪声后的,周期图法功率频谱分析');grid;%-------------------------------------------------------------------------pw2=temp.*conj(temp)/1024;%对temp做向量的共轭乘积k=0:length(temp)-1;w=2*pi*k/1024;figure(2);plot(w/pi/2,pw2);%输出功率谱函数pw2图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x1自相关法功率谱估计');grid;（3）、采点为n=1024时，功率谱图像。由放大后的图像可以清楚的得到峰值点的坐标，X1=0.2002，X2=0.2129，这样也不会出现峰值的重合现象，能够清楚的估计信号的频率。5、分析原因当采样的点数为N=128时，此时采样的得到的图像分辨力很低，并且分辨率也比较低，这就导致了功率谱图像只能看到一个峰值点。采样点数为N=512时，此时，分辨力和分辨率比较高，可以清楚的区分到两个峰值点的横坐标，此时的横坐标就是信号的频率。采样点数N=1024时，分辨力和分辨率更高，得到的信号频率更接近真实的频率，但是这是以牺牲效率为代价的，采样的点数越多，所花的时间越长，这在实际的工程中是不切合实际的，因此，在我们估计随机信号的频率的时候，要合理的采取样本点数，尽可能的采取多的样点，来接近真实的信号频率，也要考虑实际的效率问题。四、总结通过这次作业，我加强了对傅里叶变换、离散随机过程求非参数化功率谱估计的认识，并初步学会了如何在matlab上对数学问题进行仿真实验，提高了动手能力。由仿真结果知道须知万事需实践，才能有发言权。感谢陈老师对我们的严格要求和期望，万分感谢！今后一定要严格要求自己，培养严谨的学术态度、求真的学习热情，以期在学术道路上能够走得更远。