% Code accompanying fall 2025 mini project #2 to apply two level DWT
% (discrete wavelet transform) and inverse
% Originally written by Dan Jacobellis 10/24/2025

% The coefficients below for a dyadic perfect reconstruction filterbank
% are based on the Cohen–Daubechies–Feauveau wavelet (bior6.8)
% If the wavelet toolbox is installed, you can also use the following code:
% [LA, HA, LS, HS] = wfilters('bior6.8');
% LA = LA(2:end)/sqrt(2); HA(2:end) = HA/sqrt(2);
% LS= sqrt(2)*LS(2:end); HS = sqrt(2)*HS(2:end);
coeffs = [
       0.00134974786501001                         0                         0      -0.00269949573002003
      -0.00135360470301001                         0                         0      -0.00270720940602003
       -0.0120141966670801        0.0102009221870399        0.0204018443740798        0.0240283933341602
       0.00843901203981008       -0.0102300708193699        0.0204601416387398        0.0168780240796202
        0.0351664733065404       -0.0556648607799594        -0.111329721559919       -0.0703329466130807
       -0.0546333136825205        0.0285444717151497       -0.0570889434302994        -0.109266627365041
       -0.0665099006248407         0.295463938592917         0.590927877185834         0.133019801249681
         0.297547906345713        -0.536628801791565          1.07325760358313         0.595095812691426
         0.584015752240756         0.295463938592917         0.590927877185834         -1.16803150448151
         0.297547906345713        0.0285444717151497       -0.0570889434302994         0.595095812691426
       -0.0665099006248407       -0.0556648607799594        -0.111329721559919         0.133019801249681
       -0.0546333136825205       -0.0102300708193699        0.0204601416387398        -0.109266627365041
        0.0351664733065404        0.0102009221870399        0.0204018443740798       -0.0703329466130807
       0.00843901203981008                         0                         0        0.0168780240796202
       -0.0120141966670801                         0                         0        0.0240283933341602
      -0.00135360470301001                         0                         0      -0.00270720940602003
       0.00134974786501001                         0                         0      -0.00269949573002003
];
LA = coeffs(:,1); % Lowpass Analysis
HA = coeffs(:,2); % Highpass Analysis
LS = coeffs(:,3); % Lowpass Synthesis
HS = coeffs(:,4); % Highpass Synthesis

% Example audio file built into matlab.
audiodata = load('handel.mat'); 
x = audiodata.y(1:end-1); fs = audiodata.Fs;

% Analysis (level 1)
L1 = downsample(conv(x,LA),2);
H1 = downsample(conv(x,HA),2);

% Analysis (level 2)
L1L2 = downsample(conv(L1,LA),2);
L1H2 = downsample(conv(L1,HA),2);

% Synthesis
delay = length(coeffs)-1;

L1rec = conv(upsample(L1L2,2),LS) + conv(upsample(L1H2,2),HS);
L1rec = L1rec(delay:end-delay-1);

xrec = conv(upsample(L1rec,2),LS) + conv(upsample(H1,2),HS);
xrec = xrec(delay:end-delay-1);

% verify that xrec is the same as x (accounting for small numerical error)
assert( max(abs(x-xrec)) < 1e-10 )

% create spectrograms for original signal and components
figure; spectrogram(x,1024,0,1024,'yaxis',fs);
clim([-70,0]); title('Original signal')

highest_octave = conv(upsample(H1,2),HS);
figure; spectrogram(highest_octave,1024,0,1024,'yaxis',fs);
clim([-70,0]); title('Highest octave component');

middle_octave = conv(upsample(L1H2,2),HS);
middle_octave = conv(upsample(middle_octave,2),LS);
figure; spectrogram(middle_octave,1024,0,1024,'yaxis',fs);
clim([-70,0]); title('Middle octave component');

low_component = conv(upsample(L1,2),LS);
figure; spectrogram(low_component,1024,0,1024,'yaxis',fs);
clim([-70,0]); title('Low frequency component');



