欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

随机信号处理之功率谱估计中均方误差随信噪比变化的情况

程序员文章站 2022-03-29 20:00:29
...

本博客是在之前一篇博客Levinson-Durbin递推算法的基础上,另外编写了一个程序,以N=256点,P=128阶为例,分析频率估计的均方误差和信噪比之间的关系。由未加噪声的信号使用自相关法求出信号功率Power,分别取噪声功率为1-10000即(标准差a为1-100),转化为信噪比关系:
随机信号处理之功率谱估计中均方误差随信噪比变化的情况
每种噪声标准差情况下进行500次求功率谱并寻峰进行频率估计,求出对应的频率估计的均方误差并画出与SNR的关系曲线:
随机信号处理之功率谱估计中均方误差随信噪比变化的情况
MATLAB程序结果:
随机信号处理之功率谱估计中均方误差随信噪比变化的情况
Python程序结果:
随机信号处理之功率谱估计中均方误差随信噪比变化的情况
实验结果与理论情况相符合,MSE 都有一个随着 SNR 增大而减小的过程,f1相对于f2 ,受到噪声影响较大,随着SNR 的减小,其MSE升高的更早。
MATLAB程序:

clc
clear
%构造加入高斯白噪声的原函数
f1=0.1;f2=0.13;
N=256;
p=128;
n=1:N;
xn=2*sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4);

%求未加噪声的信号的自相关函数
Rx=xcorr(xn);
%求自相关的DFT,即功率谱密度,然后求得信号功率
Sx=abs(fft(Rx,N));
Power=sum(Sx);

msef1=zeros(1,100);
msef2=zeros(1,100);
SNR=zeros(1,100);
for sigma=1:100
    wucha1=0;
    wucha2=0;
    count=500;
    for cishu=1:count
        wn=sigma*randn(1,N);
        yn=xn+wn;
         %求自相关函数,取实际的Rx(0)-Rx(p)部分
        Rm=xcorr(yn,'biased');
        R=zeros(1,p+1);
        for i=1:p+1   %其中索引从1开始,故用R(1)--R(p+1)代表Rx(0)--Rx(p)
            R(i)=Rm(N+i-1);
        end
        %LevinsonDurbin递推
        a=zeros(p);
        k=zeros(1,p);
        rou=zeros(1,p);
        a(1,1)=-R(2)/R(1); 
        k(1)=a(1,1);
        rou(1)=R(1)*(1-k(1)^2);%1阶初始条件,实际初始条件是0阶,但是为了避免索引为0,从1阶开始使用
        for m=2:p
            a(m,m)=-R(m+1)/rou(m-1);%用作初始条件,下面加上累加项
            for i=1:m-1
                a(m,m)=a(m,m)-(a(m-1,i)*R(m+1-i))/rou(m-1);
            end
            k(m)=a(m,m);
            if k(m)>1
                break;
            end
            for j=1:m-1
                a(m,j)=a(m-1,j)+k(m)*a(m-1,m-j);
            end
            rou(m)=rou(m-1)*(1-(k(m))^2);
        end

        %得到p+1个参数,求频率响应
        G2=rou(p);
        ap=zeros(1,p);
        for k=1:p
            ap(k)=a(p,k);
        end
        [H,w] = freqz(G2^0.5,[1,ap],N);
        Hf = abs(H);
        Sx=Hf.^2;

        %估计f1,f2
        [pks,locs] = findpeaks(Sx,'SortStr','descend');%寻峰并降序排列
        f1_guji=locs(2)/N;
        f2_guji=locs(1)/N;
        wucha1=wucha1+(f1_guji-f1).^2;
        wucha2=wucha2+(f2_guji-f2).^2; %500次的误差平方累加
    end
    msef1(sigma)=wucha1/count;
    msef2(sigma)=wucha2/count;
    SNR(sigma)=10*log10(Power/(sigma*sigma));
end
plot(SNR,msef1,'-b');
hold on
plot(SNR,msef2,'-r');
xlabel('SNR');
ylabel('MSE');
legend('mse(f1)', 'mse(f2)');

Python程序:

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
import math

#自定义自相关函数
def xcorr(data):
    length=len(data)
    R=[]
    for m in range(0,length):
        sum = 0.0
        for n in range(0,length-m):
            sum=sum+data[n]*data[n+m]
        R.append(sum/length)
    return R

N=256
p=128

#构造原函数
f1=0.1
f2=0.13
n=np.linspace(0,N-1,N)
xn=2*np.sin(2*np.pi*f1*n+np.pi/3)+10*np.sin(2*np.pi*f2*n+np.pi/4)

#求自相关
Rx=xcorr(xn)
#求信号功率
Sx=abs(np.fft.fft(Rx))
Power=sum(Sx)

#加入高斯白噪声
msef1=[]
msef2=[]
SNR=[]
for sigma in range(1,100):
    wucha1=0.0
    wucha2=0.0
    count=500
    for cishu in range(1,count):
        wn=np.random.normal(0,sigma,len(xn))
        yn=[]
        for x,y in zip(xn,wn):
            z=x+y
            yn.append(z)

        #求自相关
        R=xcorr(yn)

        #递推初始条件
        a=np.zeros([p+1,p+1])
        rou=np.zeros(p+1)
        k=np.zeros(p+1)
        rou[0]=R[0]
        k[0]=1
        a[1][1]=-R[1]/rou[0]
        k[1]=a[1][1]
        rou[1]=rou[0]*(1-(k[1]**2))

        #递推
        for m in range(2,p+1):
            a[m][m]=-R[m]/rou[m-1]
            for i in range(1,m):
                a[m][m]=a[m][m]-(a[m-1][i]*R[m-i])/rou[m-1]
            k[m]=a[m][m]
            if k[m]>1:
                break
            for i in range(1,m):
                a[m][i]=a[m-1][i]+k[m]*a[m-1][m-i]
            rou[m]=rou[m-1]*(1-(k[m]**2))

        #得到p+1个参数,求频率响应
        G2=rou[p]
        ap=[]
        for k in range(1,p+1):
            ap.append(a[p][k])
        apk=np.insert(ap,0,[1])
        G=np.sqrt(G2)

        #计算频率响应
        w,h=scipy.signal.freqz(G,apk,worN=N)
        Hf=abs(h)
        Sx=Hf**2

        #估计f1,f2
        s=int(0.115*2*N)
        f1_guji=(np.argmax(Sx[0:s]))/(2*N)
        f2_guji=(np.argmax(Sx[s:int(N/2)])+s)/(2*N)
        wucha1=wucha1+(f1_guji-f1)**2
        wucha2=wucha2+(f2_guji-f2)**2
    msef1.append(wucha1/count)
    msef2.append(wucha2/count)
    snr=10*math.log(Power/(sigma*sigma),10)
    SNR.append(snr)

plt.plot(SNR,msef1,"-b")
plt.plot(SNR,msef2,"-r")
plt.xlabel("SNR")
plt.ylabel("MSE")
plt.title("blue:f1 red:f2")
plt.show()

其他几种算法也都有进行MSE随SNR曲线的求解和绘制,与此相类似,故不再赘述,可以查看以下资源:
现代法功率谱估计以及MSE随SNR曲线