%%%%%諧波小波變換的二帶和小波包實現%%%%%%%%%%
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));