Resolve RTSD-265

Closes RTSD-265

applications/lofar2/model/pfb_bunton_annotated/reconstruct.m

@@ -8,47 +8,65 @@
 % Correct = 1 applies the correction filter to the reconstructed signal
 % outData
 % John Bunton 2015
+%
+% Usage from README.txt:
+% . first run Gen_filter12 to create prototype lowpass filter cl
+% . then run reconstruct with ImpulseResponse = 1 and Correct = 1 to calculate
+%   an impulse response and a correction filter for the analysis PFB and
+%   synthesis PFB transfer function error
+% . then rerun reconstruct with ImpulseResponse = 0 and Correct = 0 or 1
+% . use rng('default') to check SNR results, comment rng() to experiment with
+%   new noise input evry rerun
 
 close all
 
+% Uncomment rng() line or run rng() line first in shell, to produce same random
+% numbers as with MATLAB restart, to verify the SNR results for Correct = 0, 1.
+rng('default')
+
 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 > 0.5;
-z = 2 * (z - 0.5);
-
 % cl is the prototype lowpass filter
 Ncoefs = length(cl);  % = 2304
 Ntaps = 12;
 Nfft = Ncoefs / Ntaps;  % = 192
 
+% Input data
+% 1) random block wave signal
 % Create random block wave signal with block size Nhold = 16 samples, first
 % insert 16-1 zeros between every sample in z. This yield length
 % (len - 1) * 16 + 1. Then convolve with ones(16, 1) to replicate each value
 % of z instead of the zeros. The length of conv(A, B) is length(A) +
 % length(B) - 1, so length of inData becomes ((len - 1) * 16 + 1) + (16 - 1) =
 % len * 16 = 80000.
+len = 5000;
+z = rand(len,1);
+z = z > 0.5;
+z = 2 * (z - 0.5);
 clear inData
 Nhold = 16;
 inData(1:Nhold:len*Nhold) = z;
 inData = conv(inData, ones(Nhold, 1));
-k = len * Nhold;  % = length(inData)
+k = length(inData);  % = len * Nhold
+
+% 2) random normal noise signal
+%inData = randn(1, k);  % yields similar SNR as random block wave signal
 
-%inData = cos(((1:k).^2 ) * len / (k)^2);  % comment out these two line to have random signal input
+% 3) sine wave signal
+%inData = cos(((1:k).^2 ) * len / (k)^2);
 %inData = cos(2*pi*(1:k) * 0.031 / 16);
 
+% 4) impulse signal
 if ImpulseResponse == 1
-    inData = inData .* 0;   % for testing impulse response of system
-    inData(k / 2) = 1;  % and allows calculation of correction filter
+    inData = inData .* 0;  % for testing impulse response of system
+    inData(k / 2) = 1;     % and allows calculation of correction filter
 end
 
 % Set up parameters for the analysis and synthesis filterbank
 Ros = 32 / 24;
 % fix() rounds to nearest integer towards zero.
 Nstep = fix(Nfft / Ros);  % = 144
-clf
 
 % these is the actual analysis and synthesis operations
 outBins = polyphase_analysis(inData, cl, Nfft, Nstep);
@@ -62,8 +80,8 @@ adjust = adjust * Nfft;  % normalize gain per bin for synthesis LPF
 outData = adjust * outData;
 
 if ImpulseResponse == 1
-     % correction filter to fix frequency response error, also
-     % corrects gain because input inData has unit level
+    % correction filter to fix frequency response error, also
+    % corrects gain because input inData has unit level
     correct = real(ifft(1./fft(outData)));
 end
 
@@ -78,7 +96,7 @@ corr = correct;
 % Correct for time shift for transfer function correction filter
 corr(1+k:length(corr)) = corr(1:length(corr)-k);
 
-% Correct for errror in transfer function
+% Correct for error in transfer function
 if Correct == 1
     outData = ifft(fft(outData) .* fft(fftshift(corr)));
 end
@@ -87,10 +105,16 @@ end
 diff = outData - inData(1:length(outData));
 
 % SNR
-% First do rng('default') to reset the random number generator, then
-% run reconstruct. The original code by John Bunton then yields:
+% First do rng('default') to reset the random number generator, then run
+% reconstruct. The original code by John Bunton for the random block wave wave
+% input signal then yields:
 % . rng('default'); reconstruct, with Correct = 0: SNR = 30.56 dB
 % . rng('default'); reconstruct, with Correct = 1: SNR = 62.24 dB
+% Note:
+% . With the impulse data input the SNR = 304.46 dB, so perfect.
+% . For other input data there is still some error, due to the limited stopband
+%   attenuation of the prototype LPF (can be seen using DEVS = 0.003 instead of
+%   0.0003 in Gen_filter12).
 SNR = 10 * log10(sum(abs(inData(10000:70000) .^ 2)) / sum(abs(diff(10000:70000) .^ 2)));
 sprintf('SNR = %6.2f [dB]', SNR)