From 269637ccae7624b570ff7b039922aab761d338fd Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Wed, 14 Aug 2024 09:53:40 +0200
Subject: [PATCH] Clarify input data and expected SNR.
---
.../model/pfb_bunton_annotated/reconstruct.m | 54 +++++++++++++------
1 file changed, 39 insertions(+), 15 deletions(-)
diff --git a/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m b/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m
index 22c31fcdca..c4f1200dcd 100644
--- a/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m
+++ b/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)
--
GitLab