gonglvpu.asv

小波变换用于半导体激光器可靠性分析的全部程序。我只完成了其中一部分

%产生fbm信号，H=0.3
%初始化randn发生器
randn('state',1);
%设置参数H和采样长度lg
H=0.3;
lg=1000;
x1=wfbm(H,lg);
%进行零均值处理
m1=mean(x1);
x=x1-m1;
t=1:lg;
%画出1/f信号时域波形图
figure(1);
subplot(2,1,1);
plot(t,x);
%画出1/f信号的功率谱
fs=1000;
nfft=100;
xdata=(x-mean(x))/std(x,1);
subplot(2,1,2);
y=abs(psd(xdata,nfft));
loglog((0:nfft/2-1)/nfft*fs,y(1:nfft/2));
axis([10,600,0.05,40]);
%生成白噪声信号
N=1000;
t=1:N;
y1=5*(rand(1,N)-0.5);
m2=mean(y1);
y=y1-m2;
%在1/f信号上加上白噪声
s=x+y;
t=1:1000;
%画出含有白噪声的1/f信号的时域波形图及其功率普
figure(2);
subplot(2,1,1);
plot(t,s,'LineWidth',2);
ydata=(y-mean(y))/std(y,1);
subplot(2,1,2);
y1=abs(psd(ydata,nfft));
loglog((0:nfft/2-1)/nfft*fs,y1(1:nfft/2));
%用db10小波对合成信号进行5层分解
[c,l]=wavedec(s,5,'db10');;
%从小波分解结构[c,l]中提取低频概貌系数
ca5=appcoef(c,l,'db10',5);
%从小波分解结构[c,l]中提取高频细节系数
cd5=detcoef(c,l,5);
cd4=detcoef(c,l,4);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
%求合成信号小波系数的样本方差
v5=sum((cd5).^2)/length(cd5);
v4=sum((cd4).^2)/length(cd4);
v3=sum((cd3).^2)/length(cd3);
%估计谱参数r的值
r=log2(abs((v5-v4)/(v4-v3)));
vx=(v5-v4)/((2.^(4*r))*((2.^r)-1));
vy=v4-(vx*(2.^(4*r)));
%重构小波系数
ra5=abs((vx*(2^(5*r)))/((vx*(2^(5*r))+vy)))*ca5;
rd5=abs((vx*(2^(5*r)))/((vx*(2^(5*r))+vy)))*cd5;
rd4=abs((vx*(2^(4*r)))/((vx*(2^(4*r))+vy)))*cd4;
rd3=abs((vx*(2^(3*r)))/((vx*(2^(3*r))+vy)))*cd3;
rd2=abs((vx*(2^(2*r)))/((vx*(2^(2*r))+vy)))*cd2;
rd1=abs((vx*(2^(1*r)))/((vx*(2^(1*r))+vy)))*cd1;
%
c1=[ra5,rd5,rd4,rd3,rd2,rd1];
rx=waverec(c1,l,'db10');
figure(3);
subplot(2,1,1);
plot(t,rx);
rxdata=(rx-mean(rx))/std(rx,1);
subplot(2,1,2);
y2=abs(psd(rxdata,nfft));
loglog((0:nfft/2-1)/nfft*fs,y2(1:nfft/2));r
axis([10,300,0.01,80]);