rsd1.m

单频加白噪声的随机信号分析模型

clc;

f0=300;   %信号频率
fs=5000; %采样频率
Ts=1/fs;
len=512;  %采样点数 
n=1:len;
p_max=100; %AR模型的最高阶数
snr=0; %信噪比

x0=sin(2*pi*f0*Ts*n); 
x=awgn(x0,snr);   %加噪信号源

%%%%%%%画出原始加噪信号源波形%%%%%%%
figure(1)
plot(n/fs,x)  
legend('SNR=0')
xlabel('时间/秒')
title('加噪正弦信号源')

 %%%%%%%自相关图%%%%%%%
[c_x,lags]=xcorr(x,x);           
c_x=c_x/len;
figure(2)
plot(lags,c_x)
xlabel('n')
title('信号自相关')

%%%%汉明窗周期图%%%%%%%
figure(3)
periodogram(x,hamming(len),'onesided',len,fs)
title('汉明窗周期图')

%%%%%%%Levinson公式和AIC信息准则定阶数，求AR（p）系数%%%%%%%
e=zeros(1,p_max);
a=zeros(p_max+1,p_max+1);
P=1:p_max;
for order=1:p_max
   [a(order,1:(order+1)),e(order)]=levinson(c_x(len:2*len-1),order);
   if e(order)<0
       printf('e<0,error!');
       break;
   end
end
AIC=len*log(e)+2*P;
order=find(AIC==min(AIC)) %寻找AIC最小信息量所对应的阶数order
e1=e(order)               %order所对应的模型误差
a=a(order,1:(order+1))    %order所对应的模型系数

%%%%%%%画出AR(p)模型时域估计信号波形%%%%%%%
noise=randn(1,len);
noise_dr=sqrt(e1)*noise;
x_filter=filter(1,a,noise_dr);
figure(4)
plot(n/fs,x_filter)
xlabel('时间/秒')
title('AR模型时域信号波形')

 %%%%%%%AR模型时域信号自相关图%%%%%%%
[c_x_AR,lags_AR]=xcorr(x_filter,x_filter);           
c_x_AR=c_x_AR/len;
figure(5)
plot(lags_AR,c_x_AR)
xlabel('n')
title('信号自相关')

%%%%%%%画出AR(p)模型功率谱%%%%%%%
figure(6)
[h,w]=freqz(e1,a,len*4);
plot(w/(2*pi)*fs,abs(h))
xlabel('频率/Hz')
title('AR模型频谱图')