Clarify gain adjust. Add SNR for reference.

applications/lofar2/model/pfb_bunton/reconstruct.m

@@ -10,7 +10,7 @@
 % 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;
@@ -46,8 +46,18 @@ 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
@@ -72,6 +82,14 @@ 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))))