谐波小波包分解与重构附程序

taodede

%%%%%谐波小波变换的二带和小波包实现%%%%%%%%%%
clear
fs=4000;
fn=fs/2;%奈奎斯特频率
dett=1/fs;%采样周期
nx=400;%采样点数
tp=nx*dett;%采样时间长度
detf=1/tp;%频率分辨率
t=0:dett:(nx-1)*dett;%时间向量
f=0:detf:(nx-1)*detf;%频率向量
%%%%信号部分%%%%
x=sin(100*pi*t);
x(90:170)=0.98*x(90:170);
subplot(411);plot(t,x);
fft_x=fft(x);
subplot(412);stem(f,abs(fft_x));
%%%%%%%%%%%%%%%%%%%%%%二带实现%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%k=6;%最高分解层数，本部分中低层对应高频带
%%%%二带频带划分%%%%
%for j=1:k   
%fmax(j,:)=fs/2^j;
%fmin(j,:)=fs/2^(j+1);
%end
%m=fmin;
%n=fmax;
%p=floor(fmin/detf);%谐波小波频域最低点
%q=ceil(fmax/detf);
%%%构造谐波小波频域序列采样点%%%
%for i=1:k
%c1{i}=zeros(1,p(i,:));
%c2{i}=ones(1,q(i,:)-p(i,:));
%c3{i}=zeros(1,nx-q(i,:));
%c{i}=[c1{i},c2{i},c3{i}];
%end
%%%构造频谱序列%%%%
%for ii=1:k
%har_mol(ii,:)=1/((fmax(ii,:)-fmin(ii,:))*2*pi);%谐波小波频域幅值
%har_fre(ii,:)=har_mol(ii,:)*c{ii};%谐波小波频域特性
%end
%subplot(413);stem(f,har_fre);
%%%%计算小波系数%%%%
%for iii=1:k
%aw(iii,:)=fft_x.*har_fre(iii,:);%频域小波系数
%a(iii,:)=ifft(aw(iii,:));
%signal(iii,:)=real(a(iii,:))*4*pi*(n(iii,:)-m(iii,:));
%end
%subplot(413);plot(t,signal(6,:));
%%%%%%%%%%%%%%%%%%%%小波包实现%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=3;%最高分解层数
k1=2*k;%频带个数，小波包频带数为层数的2倍
detf_pack=fn/k1;%小波包中频率间隔
%%%小波包频带划分，大k1值对应高频成分，与二带情况相反%%%
for j=1:k1   
fmax(j,:)=j*detf_pack;
fmin(j,:)=(j-1)*detf_pack;
end
m=fmin;
n=fmax;
p=floor(fmin/detf);%谐波小波频域最低点
q=ceil(fmax/detf);
%%%构造谐波小波频域序列采样点%%%
for i=1:k1
c1{i}=zeros(1,p(i,:));
c2{i}=ones(1,q(i,:)-p(i,:));
c3{i}=zeros(1,nx-q(i,:));
c{i}=[c1{i},c2{i},c3{i}];
end
%%%%构造频谱序列%%%%
for ii=1:k1
har_mol(ii,:)=1/((fmax(ii,:)-fmin(ii,:))*2*pi);%谐波小波频域幅值
har_fre(ii,:)=har_mol(ii,:)*c{ii};%谐波小波频域特性
end
%subplot(413);stem(f,har_fre);
%%%%计算小波系数%%%%
for iii=1:k1
aw(iii,:)=fft_x.*har_fre(iii,:);%频域小波系数
a(iii,:)=ifft(aw(iii,:));
signal(iii,:)=real(a(iii,:))*4*pi*(n(iii,:)-m(iii,:));
end
subplot(413);plot(t,signal(6,:));

[ca,cd]=dwt(x,'db3');
subplot(414);plot(abs(cd));

taodede

x(90:170)=0.98*x(90:170);

signal(iii,:)=real(a(iii,:))*4*pi*(n(iii,:)-m(iii,:));


这俩行究竟是什么意思呢？大神们谁能告诉我，在线等。

taodede

谁知道呢？

雨夜

1. x(90:170)=0.98*x(90:170);

复制代码

这一句应该只是用于构造原始信号的,即从第90个数据开始,其幅值比前面的略小一点

1. signal(iii,:)=real(a(iii,:))*4*pi*(n(iii,:)-m(iii,:));

复制代码

这一句应该是重构信号的纵坐标

taodede

谢谢你的帮助