gongzhensnr5.m

产生随机共振现象的输入输出信噪比曲线

clear all
clc
A=0.3;  %信号幅度
fs=10;  %采样频率
f=0.01; %信号频率
Ts=1/fs;%采样时间
h=1/fs; %时间步长
t=0:Ts:4095*Ts;
%双稳态系统参数
a=1;
b=1;
s=A*sin(2*pi*f*t);         %信号
for kk=1:1500
    disp(kk);
SNR(kk)=kk/50-20;
g=SNR(kk);
x1=AWGN(s,g);  %噪声
%输入无噪信号傅立叶变换
y=fft(s,4096); 
pyy=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
max=1;
for i=1:2049
    if(pyy(i)>pyy(max))
        max=i;
    end
end
%输入加噪信号傅立叶变换
y=fft(x1,4096); 
pyy=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
%输入信噪比计算
Bwin=0;
for i=1:2049
    Bwin=pyy(i)+Bwin;
end
Bwin=(Bwin-pyy(max));
SNRin(kk,1)=10*log(pyy(max)/Bwin)/log(10);
%四阶龙格库塔法对双稳态输出信号求解
x=sr(a,b,h,x1);
%输出信号求傅立叶变换
y=fft(x,4096); 
py=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
%输出信噪比计算
Bwout=0;
for i=1:2049
    Bwout=py(i)+Bwout;
end
Bwout=(Bwout-py(max));
SNRout(kk,1)=10*log(py(max)/Bwout)/log(10);
Q(kk,1)=SNRout(kk,1)-SNRin(kk,1);
end 
figure(1);
plot(SNR,SNRin);grid on;
hold on;
plot(SNR,SNRout,'red');
figure(2);
plot(SNR,Q);grid on;