求助电信号模极大滤波程序
主程序%____模极大值法去噪：
%二进制小波变换（离散平稳小波变换）并找到模极大序列,对极大值序列处理后重建信号;
%参数c, O_d4, O_d3可以调节。
function Denoise_w_Mod_sim_1(snr,aaa)
close all;
clc;
clear;
snr=5;
init=2055615866;
%=wnoise(3,10,snr,init);
load leleccum;
xref=leleccum(1:3920);
x = xref
signal=x;
points=1024; level=4; sr=360; num_inter=6; wf='db3';
%所处理数据的长度 分解的级数 抽样率 迭代次数 小波名称
offset=0;
%____进行二进制小波变换（离散平稳小波变换），并给出各级波形：
=wfilters(wf);
= swt(signal,level,Lo_D,Hi_D);
figure;
subplot(level,1,1); plot(real(signal)); grid on;axis tight;
for i=1:level
subplot(level+1,2,2*(i)+1);
plot(swa(i,:)); axis tight;grid on;xlabel('time');
ylabel(strcat('a ',num2str(i)));
subplot(level+1,2,2*(i)+2);
plot(swd(i,:)); axis tight;grid on;
ylabel(strcat('d ',num2str(i)));
end
%____求小波变换的模极大值及其位置,并按级给出小波变换模极大的波形:
% swa:小波概貌;swd:小波细节;
% ddw:局部极大位置; wpeak:小波变换的局部极大序列。
ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
posw=swd.*(swd>0);
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
negw=swd.*(swd<0);
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
ddw=pddw|nddw;
ddw(:,1)=1;
ddw(:,points)=1;
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;
figure;
subplot(level+1,1,1); plot(real(signal)); grid on;axis tight;
for i=1:level
subplot(level+1,1,i+1);
plot(wpeak(i,:)); axis tight;grid on;
ylabel(strcat('j= ',num2str(i)));
end
%____进行模极大值的处理：
C=152;
%此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。
D4_wpeak=wpeak(level,:);
M=max(D4_wpeak);
Thr=C*M/level; %阈值计算，可参考论文："3mm波段脉冲雷达系统研究和小波去噪分析"。
D4_wpeak=D4_wpeak.*(abs(D4_wpeak)>Thr);
%模极大值的处理方式：
%在尺度j上极大值点位置，构造一个搜索区域，
%在尺度j-1中，极大值点落在该区域的点保留，其他的置0；
D3_wpeak=wpeak(level-1,:);
D4_p=(D4_wpeak~=0);
O_d4=1;%该参数确定在上一级搜索极大值的范围，可以调整。
for P_d4=O_d4:(length(D4_wpeak)-O_d4);
if D4_p(P_d4)==1;
for i=1:O_d4-1;
D4_p(P_d4-i)=1;
end ;
end;
end;
D3_wpeak=D3_wpeak.*D4_p;
D2_wpeak=wpeak(level-2,:);
D3_p=(D3_wpeak~=0);
O_d3=1;%该参数确定在上一级搜索极大值的范围，可以调整。
for P_d3=O_d3:(length(D3_wpeak)-O_d3);
if D3_p(P_d3)==1;
for i=1:O_d3-1;
D3_p(P_d3-i)=1;
end ;
end;
end;
D2_wpeak=D2_wpeak.*D3_p;
%第一层单独处理，在第二层极大值点位置上，保留第一层相应极大值点：
D1_wpeak=wpeak(1,:);
D2_p=(D2_wpeak~=0);
D1_wpeak=D1_wpeak.*D2_p;
wpeak=;
wpeak=wpeak';
%____重构信号：
pswa=swa(level,:); % pswa: 为待重建的信号
wframe=(wpeak~=0); %迭代初始化
w0=zeros(1,points);
=swt(w0,level,Lo_D,Hi_D);
w2=d;% w2为待重建小波
for j=1:num_inter
w2=Py_Pgama(d,wpeak,wframe,1,sr); % 先进行Py投影和 Pgama投影
w0=iswt(pswa,w2,Lo_R,Hi_R); % 再进行Pv投影
=swt(w0,level,Lo_D,Hi_D); % Pv
end
pswa=iswt(swa(level,:),w2,Lo_R,Hi_R); % 计算重建信号
xcrr=xref-pswa; % 重建误差
figure,
subplot(411)
plot(xref(1:points),'r');
axis();
subplot(412)
plot(x(1:points),'r');
axis();
subplot(413)
plot(pswa(1:points));
axis();
subplot(414)
plot(xcrr(1:points));
axis();
子程序
function =P_gama(interval,lev,sr);
T=length(interval);
%该函数对一个区间进行Pgama投影，返回修正的区间
if T==2
inter=interval;
else
t=linspace(0,(T-1)/sr,T);
para=(()\')';
alpha=para(1);
beta=para(2);
inter=alpha.*exp(2^(-lev).*t)+beta.*exp(-2^(-lev).*t);
end
子程序
function pc3inte=P_y(interval,len);
% 该函数对区间进行裁减即Py投影，返回裁剪后的区间信号
if sign(interval(1))==sign(interval(len))
interval=interval.*(sign(interval)==sign(interval(1)));
inte=interp1(,,(1:len),'linear');
interval=sign(interval(1))*(abs(inte)-(abs(inte)-abs(interval)).*((abs(inte)-abs(interval))>0));
else
sgn=sign(interval(len)-interval(1));
intemax=max();
intemin=min();
for i=1:len-2
if sign(interval(i+1)-interval(i))~=sgn
interval(i+1)=interval(i);
end
if interval(i+1)>intemax
interval(i+1)=intemax;
end
if interval(i+1)<intemin
interval(i+1)=intemin;
end
end
end
pc3inte=interval;
子程序
functionw2=Py_Pgama(w1,wpeak,wframe,level,sr);
% 该函数用于进行 Pgama 和 Py 投影
err=wpeak-w1.*(wpeak~=0);
w2=zeros(size(wpeak));
=size(wpeak);
% 对每一级小波分别进行处理
for m=1:r
frame=find(wpeak(m,:));
num_interval=length(frame)-1;
% 先找到以模极大划分的区间, 然后对每一区间进行Py投影
for j=1:num_interval
interval=w1(m,frame(j):frame(j+1));
len=length(interval);
if len>2
w1(m,frame(j):frame(j+1))=P_y(interval,len);
end
end
% 再逐一区间进行Pgama投影
for j=1:num_interval
interval=err(m,frame(j):frame(j+1));
if r==1
err(m,frame(j):frame(j+1))=P_gama(interval,level,sr);
else
err(m,frame(j):frame(j+1))=P_gama(interval,m,sr);
end
end
w2(m,:)=w1(m,:)+err(m,:);
end
提问：怎么样把这个电信号弄进去，用这个滤波的程序把它滤出来。（原来的程序是给定一个信号，对它进行加噪然后滤波），我现在这个电信号本来就有噪声，用这个程序对它进行滤波，一共要出来三个图。最后一个图只要含噪的信号和滤波的信号。
有哪个大虾能指点下。。。本人不盛感激。。。
页:
[1]