Clarify adjust. Support running reconstruct directly after Gen_filter12.m.

applications/lofar2/model/pfb_bunton_annotated/reconstruct.m

@@ -11,14 +11,18 @@
 %
 % 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
+% . then run reconstruct with ImpulseResponse = 0 and Correct = 0 to filter
+%   the inData
+% . to correct the overall response, rerun 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 = 1 to also
+%   apply the correction filtering
 % . use rng('default') to check SNR results, comment rng() to experiment with
-%   new noise input evry rerun
+%   new noise input every rerun
 
 close all
+format longE
 
 % 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.
@@ -57,6 +61,7 @@ k = length(inData);  % = len * Nhold
 %inData = cos(((1:k).^2 ) * len / (k)^2);
 %inData = cos(2*pi*(1:k) * 0.031 / 16);
 %inData = cos(2*pi*(0:k-1) * 1.5 / Nfft);  % use integer for CW at bin
+%inData = (0:k-1);
 
 % 4) impulse signal
 if ImpulseResponse == 1
@@ -74,10 +79,10 @@ outBins = polyphase_analysis(inData, cl, Nfft, Nstep);
 outData = polyphase_synthesis_b(outBins, cl, Nstep);
 
 % adjust gain
-% . unit gain in prototype LPF Gen_filter12.m
-adjust = Nfft;  % normalize gain per bin for analysis LPF
-adjust = adjust / Ros;  % compensate for oversampling factor
-adjust = adjust * Nfft;  % normalize gain per bin for synthesis LPF
+% . assume unit gain in prototype LPF Gen_filter12.m
+% . adjust for 1/Nfft in synthesis IFFT and adjust for 1/Nup in synthesis
+%   interpolation, with analysis Ndown = synthesis Nup = Nstep.
+adjust = Nfft * Nstep;
 outData = adjust * outData;
 
 if ImpulseResponse == 1
@@ -86,19 +91,18 @@ if ImpulseResponse == 1
     correct = real(ifft(1./fft(outData)));
 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 for time shift for transfer function correction filter
-corr(1+k:length(corr)) = corr(1:length(corr)-k);
-
 % Correct for error in transfer function
 if Correct == 1
+    % 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 for time shift for transfer function correction filter
+    corr(1+k:length(corr)) = corr(1:length(corr)-k);
     outData = ifft(fft(outData) .* fft(fftshift(corr)));
 end
 
@@ -121,6 +125,7 @@ diff = outData - inData(1:length(outData));
 %   With DEVP = 0.005 and W0 = 1 / 15.05 then db(0.005) = -46.0 dB and SNR =
 %   38.8 dB, so increased.
 SNR = 10 * log10(sum(abs(inData(10000:70000) .^ 2)) / sum(abs(diff(10000:70000) .^ 2)));
+SNR = 20 * log10(std(inData(10000:70000)) / std(diff(10000:70000)));
 sprintf('SNR = %.5f [dB]', SNR)
 
 %plot(diff)