%%%%%谐波小波变换的二带和小波包实现%%%%%%%%%%

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));