Resolve RTSD-265

Closes RTSD-265

applications/lofar2/model/pfb_bunton/reconstruct.m

-%reconstruction of signal after polyphase filterbank 
+%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 
+% 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 = 0;           % set to 1 to allow correction filter to be applied
+Correct = 1;           % set to 1 to allow correction filter to be applied
 
 
 len=5000;
@@ -36,9 +36,9 @@ end
 
 % Set up parameters for the analysis and synthesis filterbank
 firlength=12;      %FFT length is 1/firlength of filter length
-block = length(cl)/firlength; 
+block = length(cl)/firlength;
 Oversampling_Ratio=32/24;
-step  = fix(length(cl)/(12*Oversampling_Ratio));  
+step  = fix(length(cl)/(12*Oversampling_Ratio));
 clf
 
 %%%%% these is the actual analysis and synthesis operations
@@ -46,9 +46,19 @@ out=polyphase_analysis(in,cl,block,step);
 recon=polyphase_synthesis_b(out,cl,step );
 %%%%%
 
-recon = (256*step/block)*recon;  %adjust gain
-   
-    
+% 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
@@ -62,16 +72,24 @@ end
 k=39;
 corr=correct;
 % Correct fo time shift for transfer function correction filter
-corr(1+k:length(corr))=corr(1:length(corr)-k);  
+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))); 
+    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))))
@@ -96,4 +114,4 @@ if (ImpulseResponse == 1)
     plot(recon)
 else
     plot(recon(1:10000),'r')
-end
\ No newline at end of file
+end