gongzhenq.m

gongzhenQ.m:计算双稳态随机共振系统输入输出信噪比增益

clear all
clc
fs=5;   %参样频率
f=0.01; %信号频率
Ts=1/fs;%参样时间
h=1/fs; %时间步长
t=0:Ts:4095*Ts;
%双稳态系统参数
a=1;
b=1;
s=0.3*sin(2*pi*0.01*t);         %信号
for kk=1:1000
    disp(kk);
D(kk,1)=kk/200;
Dp=D(kk,1); %噪声强度
Noise=zeros(1,length(t));
Noise=sqrt(2*Dp)*randn(size(t));
x1=s+Noise;  %噪声
%输入加噪信号傅立叶变换
y=fft(x1,4096); 
pyy=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
%输入信噪比计算
Bwin=0;
Swin=0;
for i=1:2049
    Bwin=pyy(i)+Bwin;
    if(pyy(i)>Swin)
        Swin=pyy(i);
    end
end
Bwin=(Bwin-Swin);
SNRin(kk,1)=10*log(Swin/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;
Swout=0;
for i=1:2049
    Bwout=py(i)+Bwout;
    if(py(i)>Swout)
        Swout=py(i);
    end
end
Bwout=(Bwout-Swout);
SNRout(kk,1)=10*log(Swout/Bwout)/log(10);
Q(kk,1)=(Swout/Bwout)/(Swin/Bwin);
end 
figure(1);
plot(D,SNRin);
hold on;
plot(D,SNRout,'red');
figure(2);
plot(D,Q);