clear;
% Computes 1D DWT of a sequence x[n] up to 2 levels
% and plots the DWT coefficients
% specify the name of the wavelet in ''WaveletName''
%
WaveletName = 'db8';
%WaveletName='bior3.1';
% compute level 1 DWT of x[n], which consists of the
% approximation coefficients ''a1'' and detail
% coefficients ''d1.''
L=100;
n = 0:L-1;% sequence length
A1 = 2.0; A2 = 1.5;% amplitudes of the two sinusoids
F1 = 0.1; F2 = 0.2;% frequencies normalized to sampling
%frequency
dn=[20,35,67];
nfilter=32;
x1 = A1*sin(2*pi*F1*n) + A2*cos(2*pi*F2*n);
x2= zeros(1,L);
x2(dn)=32;
b=fir1(nfilter,0.5);
% insert three low-pass filtered discrete impulses at positions 20,35,67
x2=conv(x2,b);
x2=x2(1:L);

% Make plots of x1 and x2 separately !
x=x1+x2;

figure
subplot(1,2,1),
NFFT=8;
spectrogram(x,hann(NFFT),0,NFFT,1),title('STFT with N_{FFT}=8');
subplot(1,2,2),
NFFT=32;
spectrogram(x,hann(NFFT),0,NFFT,1),title('STFT with N_{FFT}=32');
%%


[a1,d1] = dwt(x,WaveletName);
figure,subplot(3,1,1),plot(n,x,'k','LineWidth',1);
title('Input signal x[n]')
subplot(3,1,2),plot(a1,'k','LineWidth',1)
title('DWT of x[n]:Level-1 scaling coefficients')
subplot(3,1,3),plot(d1,'k','LineWidth',1)
title('Level-1 detail coefficients')
% Next, compute level 2 DWT of a1, the level 1 approximation
[a2,d2] = dwt(a1,WaveletName);
figure,subplot(2,1,1),plot(a2,'k','LineWidth',1)
title('level-2: scaling coefficients')
subplot(2,1,2),plot(d2,'k','LineWidth',1)
title('Level-2 detail coefficients');
%%
[LO_D,HI_D,LO_R,HI_R] = wfilters(WaveletName);
figure;
subplot(2,1,1),plot(LO_D),title('Low-pass analysis filter impulse response');
subplot(2,1,2),plot(HI_D),title('High-pass analysis filter impulse response');
figure;
subplot(2,1,1),plot(LO_R),title('Low-pass synthetis filter impulse response');
subplot(2,1,2),plot(HI_R),title('High-pass synthetis filter impulse response');
%%
Nfft=256;
fs=1;
fplot=0:fs/Nfft:fs;
fplot=fplot(1:Nfft/2);
LO_DF=fft(LO_D,Nfft);
HI_DF=fft(HI_D,Nfft);
LO_RF=fft(LO_R,Nfft);
HI_RF=fft(HI_R,Nfft);
LO_DF=20*log10(abs(LO_DF(1:Nfft/2)));
HI_DF=20*log10(abs(HI_DF(1:Nfft/2)));
LO_RF=20*log10(abs(LO_RF(1:Nfft/2)));
HI_RF=20*log10(abs(HI_RF(1:Nfft/2)));

figure;
subplot(2,1,1),plot(fplot,LO_DF),title('Low-pass analysis filter, freq. resp. magnitude');
subplot(2,1,2),plot(fplot,HI_DF),title('High-pass analysis filter, freq. resp. magnitude');
figure;
subplot(2,1,1),plot(fplot,LO_RF),title('Low-pass synthetis filter, freq. resp. magnitude');
subplot(2,1,2),plot(fplot,HI_RF),title('High-pass synthetis filter, freq. resp. magnitude');
xlabel('Normalized frequency, f_s=1');
%%

lev   = 8;
wname = WaveletName
nbcol = 64;
[c,l] = wavedec(x,lev,wname);
%Expand the discrete wavelet coefficients for visualization
len = length(x);
cfd = zeros(lev,len);
for k = 1:lev
    d = detcoef(c,l,k);
    d = d(:)';
    d = d(ones(1,2^k),:);
    cfd(k,:) = wkeep1(d(:)',len);
end
cfd =  cfd(:);
I = find(abs(cfd)<sqrt(eps));
cfd(I) = zeros(size(I));
cfd = reshape(cfd,lev,len);
cfd = wcodemat(cfd,nbcol,'row');
%Plot the discrete coefficients.
%colormap(pink(nbcol));
image(cfd);
tics = 1:lev;
labs = int2str((1:lev)');
ax = gca;
ax.YTickLabelMode = 'manual';
ax.YDir = 'normal';
ax.Box = 'On';
ax.YTick = tics;
ax.YTickLabel = labs;
title('Discrete Wavelet Transform, Absolute Coefficients.');
xlabel('Time ');
ylabel('Level');

figure
scales = (1:64); % Levels 1 to 5 correspond to scales 2, 4, 8, 16 and 32.
cwt(x,scales,wname,'plot');
colormap(pink(nbcol));
title('CWT of signal')
