reconstruct.m

%reconstruction of signal after polyphase filterbank
% John Bunton created Dec 2005
%
% modified to add variables "Impulse Response' and "Correct"
% Impulse Response = 1 calculates the impulse response of the system
% the result is in variable "recon".  It also causes the transfer function
% correction filter to be calculated, Result in variable "correct"
% Correct = 1 applies the correction filter to the reconstructed signal
% "recon"
% John Bunton 2015

ImpulseResponse = 0;   % set to 1 to calculate correction filter
Correct = 1;           % set to 1 to allow correction filter to be applied


len=5000;
z=rand(len,1);
z=z>.5;
z=2*(z-.5);

%cl is the prototype lowpass filter

clear in
in(1:16:len*16)=z;
in=conv(in,ones(16,1));
%or
k=len*16;

%in = cos(((1:k).^2 )*5000/(k)^2); %comment out these two line to have random signal input
%in=cos(2*pi*(1:k)*0.031/16);

if (ImpulseResponse == 1)
    in=in.*0;   %for testing impulse response of system
    in(k/2)=1;  %and allows calculation of correction filter
end

% Set up parameters for the analysis and synthesis filterbank
firlength=12;      %FFT length is 1/firlength of filter length
block = length(cl)/firlength;
Oversampling_Ratio=32/24;
step  = fix(length(cl)/(12*Oversampling_Ratio));
clf

%%%%% these is the actual analysis and synthesis operations
out=polyphase_analysis(in,cl,block,step);
recon=polyphase_synthesis_b(out,cl,step );
%%%%%

% EK
Q = 12;  % interpolation factor in Gen_filter12.m
adjust = block / Q;  % normalize gain per bin for analysis LPF
adjust = adjust / Oversampling_Ratio;  % compensate for oversampling factor
adjust = adjust * block / Q;  % normalize gain per bin for synthesis LPF
% JB
%adjust = 256*step/block;

sprintf('sum(cl) = %f', sum(cl))  % LPF DC gain, = 12.062337
sprintf('adjust = %f', adjust)  % gain adjust, = 192

recon = adjust*recon;  %adjust gain

if (ImpulseResponse == 1)
    correct=real(ifft(1./fft(recon)));  %correction filter to fix frequency response error
end

% Choose parameter k to suit oversampling ratio
% oversampling ratio 32/23 then k=0
% oversampling ratio 32/24 then k=39
% oversampling ratio 32/25 then k=72
% oversampling ratio 32/26 then k=3
% oversampling ratio 32/27 then k=48
k=39;
corr=correct;
% Correct fo time shift for transfer function correction filter
corr(1+k:length(corr))=corr(1:length(corr)-k);

% Correct for errror in transfer function
if (Correct == 1)
    recon = ifft (fft(recon).* fft(fftshift(corr)));
end

% Take the difference to see the error
diff = recon-in(1:length(recon));

% SNR
% First do rng('default') to reset the random number generator, then
% run reconstruct. The original code by John Bunton then yields:
% . rng('default'); reconstruct, with Correct = 0: SNR = 30.56 dB
% . rng('default'); reconstruct, with Correct = 1: SNR = 62.24 dB
SNR = 10 * log10(sum(abs(in(10000:70000) .^ 2)) / sum(abs(diff(10000:70000) .^ 2)));
sprintf('SNR = %6.2f [dB]', SNR)

%plot(diff)
%plot(db(fft(recon))) %to see frequency errot
%plot(db(fft(diff(length(in)/8:7*length(in)/8))))


figure(3)
hold off
plot(db(fft(in(10000:70000)/30000)),'r')
hold on
plot(db(fft(diff(10000:70000)/30000)),'b')
grid
title('Spectrum of input and error')
figure(2)
plot(diff)
title('Difference between input and reconstructed')
figure(1)
hold off
plot(in(1:10000))
hold on
title('Input and Reconstructed signal')
if (ImpulseResponse == 1)
    plot(recon)
else
    plot(recon(1:10000),'r')
end