% 采用trous算法(undecimated wavelet transform)實(shí)現(xiàn)小波變換
% 分解兩層

clear;clc;

%% 1.生成信號(hào)

f=50; % 頻率
fs=800; % 采樣率
T=128; % 信號(hào)長(zhǎng)度
n=1:T;
y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 信號(hào)

%% 2.正變換

l1=wfilters('db4','l')*sqrt(2)/2; % 參考低通濾波器
l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通濾波器1
h1=wfilters('db4','h')*sqrt(2)/2; % 參考高通濾波器
h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通濾波器1

low1=ifft(fft(y).*fft(l1_zeros)); % 低頻分量1
high1=ifft(fft(y).*fft(h1_zeros)); % 高頻分量1

l2=dyadup(l1); % 原濾波器插值
l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通濾波器2
h2=dyadup(h1); % 原濾波器插值
h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通濾波器2

low2=ifft(fft(low1).*fft(l2_zeros)); % 低頻分量2
high2=ifft(fft(low1).*fft(h2_zeros)); % 高頻分量2

%% 3.反變換

lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重構(gòu)低通濾波器2
hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重構(gòu)高通濾波器2
lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重構(gòu)低通濾波器1
hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重構(gòu)高通濾波器1

lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重構(gòu)低頻分量1(lowr=low1)
r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重構(gòu)源信號(hào)(r_s=y)

%% 4.繪圖

figure(1);
plot(y);
title('源信號(hào)');

figure(2);
plot(low1,'r');
hold on;
plot(low2,'b');
legend('第一層低頻','第二層低頻');

figure(3);
plot(high1,'r');
hold on;
plot(high2,'b');
legend('第一層高頻','第二層高頻');

figure(4);
plot(low1,'r');
hold on;
plot(lowr,'b.');
legend('第一層低頻','重構(gòu)第一層低頻');

figure(5);
plot(y,'r');
hold on;
plot(r_s,'b.');
legend('源信號(hào)','重構(gòu)信號(hào)');

disp(norm(low1-lowr))
disp(norm(y-r_s))