% E03_7.m                  19oct15  29mar21
% reconstr_dft_dct.m; wht; ht
% 
close all
%% create or get data
clear all;
 x = [ones(1,38)  zeros(1,12)];            %  rectangular signal 
% x = linspace(0,1,50);                      %  ramp
% x = [linspace(0,1,25)  linspace(1,0,25)];   % double ramp
% load ecgSignals; x = ecgnl(1031:1030+210);  % ecg data
lx = max(size(x));
ome = [0:lx-1] ./ lx * 2 * pi;   % frequency axis
return;

%% DFT
X = fft(x);  Xm = abs(X);  Xp = angle(X);   Xpu = unwrap(Xp);
figure(1);
subplot(221); plot(x); axis([1,lx,-.5,1.5]);
subplot(222);  plot(ome,Xm);  axis tight;
subplot(223);  plot(ome,Xp);   axis tight; 
subplot(224);  plot(ome,Xpu);   axis tight;  
pause();

err = [];
for i=1:lx/2+1;        % reconstruct using an increasing no. of sinusoids
    Xmn = Xm;
    Xmn(i+1:end-i+1) = 0;
    xn = real(ifft( Xmn .* exp(j*Xp) ));
    subplot(411);  plot(x); axis([1,lx,-.5,1.5]);   % original data
    subplot(412);  plot(xn);   axis([1,lx,-.5,1.5]); % progressive reconstruction
         text(5,0,['DFT - ' num2str(i)]);
    subplot(413);  stem(ome,Xmn);  axis tight;  % frequency components used in reconstr.
	 err(i) = mean((xn-x).^2);
	subplot(414);  stem(err);  axis([1 lx/2 0 .2])  % mse of the reconstruction
	if i<15 pause(); else pause(.1); end;
end;
return;

%%  compare with DCT
figure(2)
Xdct = dct(x);
err = [];
for i=1:lx;        % reconstruct using an increasing no. of DCT coeffs.
    Xdctn = Xdct;
    Xdctn(i+1:end) = 0;
    xn = idct(Xdctn);
    subplot(411);  plot(x); axis([1,lx,-.5,1.5]);
    subplot(412);  plot(xn);   axis([1,lx,-.5,1.5]);  text(5,0,['DCT - ' num2str(i)]);
    subplot(413);  stem(ome,Xdctn);  axis tight;
	 err(i) = mean((xn-x).^2);
	subplot(414);  stem(err);  axis([1 lx 0 .2])  % mse of the reconstruction
	if i<15 pause(); else pause(.1); end;
end;

%%  compare with Walsh-Hadamard Transform
figure(3)
Xwht = fwht(x);     % if input length is not 2^n, FWHT does zero padding
Xwht = Xwht(1:lx);  % ...and then we need to truncate
err = [];
for i=1:lx;      % reconstruct using an increasing no. of WHT coeffs.
    Xwhtn = Xwht;
    Xwhtn(i+1:end) = 0;
    xn = ifwht(Xwhtn);  % again, a 2^n elements sequence
    xn = xn(1:lx);      % truncated (reconstruction may not be exact)
    subplot(411);  plot(x); axis([1,lx,-.5,1.5]);
    subplot(412);  plot(xn);   axis([1,lx,-.5,1.5]);  text(5,0,['WHT - ' num2str(i)]);
    subplot(413);  stem(ome,Xwhtn);  axis tight;
	 err(i) = mean((xn-x).^2);
	subplot(414);  stem(err);  axis([1 lx 0 .2])  % mse of the reconstruction
	if i<15 pause(); else pause(.1); end;
end;

%% compare with Haar transform [NOTE Matlab >= 2016b]
% note multilevel analysis makes sense only for long vector data like ECG

[a,d] = haart(x); % approximation (a) and detail (d) Haar coeffs.
d,a                   % a is a scalar only if length(x)=2^n
NLev = max(size(d))   % d is a Matlab cell formed by NLev vectors
                      % NLev depends on length(x)
figure(41);
for i=1:NLev;      % plot coeffs. of HT at each level (from fine to coarse)
	subplot(NLev,1,i); stem(d{1,i});
    end            % [note no. of data in d+a = length(x)]
figure(42);
for i=1:NLev-1;    % reconstruct using a and a subset of levels (from all to 1)
	subplot(NLev-1,1,i); plot(ihaart(a,d,i));
	end


%% compare with Haar transform   FOR OLD MATLAB VERSIONS
% https://people.sc.fsu.edu/~jburkardt/m_src/haar/haar.html
%
figure(4)
X = x;              % direct Haar transform---------------
s = sqrt(2);
w = zeros(1,lx);
k = 1;  while (k*2 <= lx), k=k*2; end
while (1<k)
  k = floor(k/2);
  w(  1:  k) = ( X(1:2:2*k-1) + X(2:2:2*k) ) / s;
  w(k+1:k+k) = ( X(1:2:2*k-1) - X(2:2:2*k) ) / s;
  X(1:2*k) = w(1:2*k);
end                        %------------------------------

err = [];
for i=1:lx;        % reconstruct using an increasing no. of HT coeffs.
    Xn = X;
    Xn(i+1:end) = 0;
    
  xn = Xn;         
  
  
  % inverse Haar transform------------
  w = zeros(1,lx);
  m = 1;
  while (m*2 <= lx)
    w(1:2:2*m-1) = ( xn(1:m) + xn(1+m:m+m) ) / s;
    w(2:2:2*m)   = ( xn(1:m) - xn(1+m:m+m) ) / s;
    xn(1:2*m) = w(1:2*m);
    m = m * 2; 
  end                 %-------------------------------------
    
  subplot(411);  plot(x); axis([1,lx,-.5,1.5]);
  subplot(412);  plot(xn);   axis([1,lx,-.5,1.5]);   text(5,0,['HT - ' num2str(i)]);
  subplot(413);  stem(ome,Xn);  axis tight;
   err(i) = mean((xn-x).^2);
  subplot(414);  stem(err);  axis([1 lx 0 .2])  % mse of the reconstruction
  if i<15 pause(); else pause(.1); end;
end;

