diff --git a/applications/apertif/matlab/apertif_matlab_readme.txt b/applications/apertif/matlab/apertif_matlab_readme.txt
index 1e4164c3cfd841826325107281e36d22af5ab057..d4bc462b97be65376b92b7039710a316f7f39fb3 100644
--- a/applications/apertif/matlab/apertif_matlab_readme.txt
+++ b/applications/apertif/matlab/apertif_matlab_readme.txt
@@ -32,7 +32,9 @@ This directory contains Matlab code for modelling DSP aspects of Apertif:
arts_fir.png Figure FIR filter transfer function plot saved by FIR_LPF_ApertifCF.m
pfir_coeff.m Function Compute polyphase filterbank coefficients (same purpose as FIR_LPF_ApertifCF.m but with more options)
- run_pfir_coeff.m Script Run pfir_coeff.m
+ run_pfir_coeff.m Script Run pfir_coeff.m to show filter response in time and frequency domain
+
+ run_pfir.m Script Run the PFIR with real input to create HDL reference data (based on one_pfb.m)
4) FFT
@@ -50,7 +52,8 @@ This directory contains Matlab code for modelling DSP aspects of Apertif:
reorder_serial.m Function Select data from block of data
corner_turn.m Function Collect M blocks of in_data(1:K) and then transpose [M] and [K] for corner turn
- run_pfft.m Script Run data path model with only the subband PFFT (real input)
+ run_pfft.m Script Run the PFFT with real input to create HDL reference data (based on one_pfb.m)
+ run_pfft_complex.m Script Run the PFFT with complex input to create HDL reference data (based on one_pfb.m)
6) Polyphase filterbank
diff --git a/applications/apertif/matlab/run_pfir.m b/applications/apertif/matlab/run_pfir.m
new file mode 100644
index 0000000000000000000000000000000000000000..f26287a9fe301e2163a74f9d19b962231c7afb01
--- /dev/null
+++ b/applications/apertif/matlab/run_pfir.m
@@ -0,0 +1,462 @@
+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2016
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%-----------------------------------------------------------------------------
+% Author: Eric Kooistra, 2016 (for Apertif)
+%
+% Purpose : Run the PFIR with real input:
+% To create input and output reference data for HDL implementation in
+% data/, e.g.:
+%
+% . PFIR coeffcients: run_pfir_m_pfir_coeff_fircls1_16taps_128points_16b.dat
+% . PFIR coeffcients, WG, PFIR, PFFT data: run_pfir_m_chirp_8b_16taps_128points_16b_16b.dat
+%
+% Description :
+%
+% * See one_pfb.m for general description. This run_pfir.m is similar and
+% also still includes the PFFT, but its main purpose to verify the HDL
+% for the PFIR.
+clear all;
+close all;
+fig=0;
+
+%% Selectable test bench settings
+tb.model_filterbank = 'Free choice'; % 'LOFAR' = use fixed settings of LOFAR subband filterbank, else use free choice of tb.nof_subbands and PFIR coefficients
+%tb.model_filterbank = 'LOFAR';
+if strcmp(tb.model_filterbank, 'LOFAR')
+ tb.model_quantization = 'fixed point';
+ tb.nof_subbands = 512;
+else
+ tb.model_quantization = 'floating point';
+ tb.model_quantization = 'fixed point';
+ tb.nof_subbands = 64;
+end
+
+% Input signal 'sinusoid' carrier frequency
+tb.subband_wg = 47; % subband range 0:tb.nof_subbands-1, can be fraction to have any sinusoid frequency
+%tb.subband_wg = 1.55;
+%tb.subband_wg = 12;
+
+% Model a frequency sweep of the 'sinusoid'
+tb.chirp = 0; % 0 = use fixed tb.subband_wg frequency, else increment WG frequency every block to have chirp frequency sweep
+tb.chirp = 1;
+if tb.chirp
+ tb.nof_tsub = 200; % number of subband periods to simulate
+else
+ tb.nof_tsub = 100; % number of subband periods to simulate
+end
+
+tb.plot_per_block = 0; % 1 = plot spectrum for each block in time, else skip this plot to save time
+%tb.plot_per_block = 1;
+
+
+%% Derived test bench settings
+tb.nof_complex = 2;
+tb.subband_i = floor(tb.subband_wg); % natural subband index in range 0:tb.nof_subbands-1
+tb.subband_I = tb.subband_i+1; % equivalent Matlab indexing for range 1:tb.nof_subbands
+if tb.subband_i<tb.nof_subbands-1 % determine index of next neighbour subband
+ tb.subband_iplus1 = tb.subband_i+1;
+else
+ tb.subband_iplus1 = 0; % wrap subband index tb.nof_subbands to subband index 0
+end
+tb.subband_Iplus1 = tb.subband_iplus1+1; % equivalent Matlab index for next neighbour subband in range 1:tb.nof_subbands
+tb.subband_fft_size = tb.nof_complex*tb.nof_subbands; % subband filterbank real FFT
+
+fs = 1; % normalized sample frequency
+fsub = fs/tb.subband_fft_size; % subband frequency relative to fs
+
+disp(sprintf (['Test bench settings:\n', ...
+ '. Model filterbank = %s\n', ...
+ '. Model coefficients and data = %s\n', ...
+ '. Number of subbands = %d\n', ...
+ '. Number of subband periods to simulate = %d\n', ...
+ '. Subband filterbank real FFT size = %d\n', ...
+ '. WG subband frequency = %f\n'], ...
+ tb.model_filterbank, tb.model_quantization, tb.nof_subbands, tb.nof_tsub, tb.subband_fft_size, tb.subband_wg));
+
+
+%% Waveform generator
+ctrl_wg.signal = 'sinusoid';
+ctrl_wg.complex = false;
+ctrl_wg.block_nr = 0;
+ctrl_wg.block_size = tb.subband_fft_size;
+ctrl_wg.sop = 0;
+ctrl_wg.psk = 0;
+if strcmp(tb.model_quantization, 'floating point')
+ ctrl_wg.data_w = 0;
+ lsb = 1/2^8; % Assume 8 bit ADC to have reference for AGWN level
+else
+ ctrl_wg.data_w = 8;
+ lsb = 1/2^ctrl_wg.data_w;
+end
+ctrl_wg.q_full_scale = 2^(ctrl_wg.data_w-1);
+ctrl_wg.agwn_sigma = 2.6*lsb;
+ctrl_wg.agwn_sigma = 0; % AGWN sigma
+ctrl_wg.ampl = 1; % full scale maps to 1, so ampl > 1 causes analogue clipping
+if ctrl_wg.agwn_sigma>0
+ ctrl_wg.ampl = 1-5*ctrl_wg.agwn_sigma;
+ ctrl_wg.ampl = 0.9;
+end
+if tb.chirp
+ ctrl_wg.df = 0.5*fsub; % increment freq by df per block to create chirp
+else
+ ctrl_wg.df = 0;
+end
+if ctrl_wg.df == 0
+ ctrl_wg.freq = tb.subband_wg*fsub; % start WG at tb.subband_wg
+else
+ ctrl_wg.freq = (tb.subband_i-1)*fsub; % start WG about one subband before tb.subband_wg to account for impulse response of several tsub
+end
+ctrl_wg.phase = 0; % Phase offset normalized to 2pi
+if ctrl_wg.freq == 0 && ctrl_wg.df == 0
+ ctrl_wg.offset = 1; % DC offset
+else
+ ctrl_wg.offset = 0; % DC offset
+end
+
+%% Subband PFIR filter parameters
+ctrl_pfir_subband.bypass = 0;
+ctrl_pfir_subband.nof_polyphases = tb.subband_fft_size;
+
+if strcmp(tb.model_quantization, 'floating point')
+ ctrl_pfir_subband.data_w = 0; % no quantization
+ ctrl_pfir_subband.coeff_w = 0; % no quantization
+else
+ ctrl_pfir_subband.data_w = 16;
+ ctrl_pfir_subband.backoff_w = 1; % attenuate PFIR data output by 2^backoff_w before quantization to account for PFIR overshoot
+ ctrl_pfir_subband.coeff_w = 16; % nof bits per coefficient including 1 sign bit
+ ctrl_pfir_subband.coeff_q_max = 2^(ctrl_pfir_subband.coeff_w-1)-1;
+ ctrl_pfir_subband.coeff_q_full_scale = 2^(ctrl_pfir_subband.coeff_w-1);
+ ctrl_pfir_subband.data_q_full_scale = 2^(ctrl_pfir_subband.data_w-1-ctrl_pfir_subband.backoff_w);
+end
+
+if strcmp(tb.model_filterbank, 'LOFAR')
+ % Load the quantized PFIR coefficients of LOFAR station (fixed point)
+ ctrl_pfir_subband.nof_taps = 16; % Number of taps
+ ctrl_pfir_subband.nof_coefficients = ctrl_pfir_subband.nof_polyphases*ctrl_pfir_subband.nof_taps; % Number of filter coefficients (taps)
+ ctrl_pfir_subband.data_w = 16;
+ ctrl_pfir_subband.config.design = 'lofar file';
+ hfir_subband_coeff = load('data/Coeffs16384Kaiser-quant.dat');
+else
+ % Calculate free choice PFIR coefficients (can be floating point or fixed point)
+ ctrl_pfir_subband.nof_taps = 16; % Number of taps
+ ctrl_pfir_subband.nof_coefficients = ctrl_pfir_subband.nof_polyphases*ctrl_pfir_subband.nof_taps; % Number of filter coefficients (taps)
+ ctrl_pfir_subband.config.design = 'fir1';
+ ctrl_pfir_subband.config.design = 'fircls1'; % 'fir1', 'fircls1'
+ ctrl_pfir_subband.config.design_flag = ''; % only for fircls1: 'trace', ''
+ ctrl_pfir_subband.config.interpolate = 'interpft'; % only for fircls1: 'resample', 'fourier', 'interpft'
+ ctrl_pfir_subband.r_pass = 1e-3; % only for fircls1
+ ctrl_pfir_subband.r_stop = 1e-4; % only for fircls1
+ ctrl_pfir_subband.hp_factor = 1.050; % Adjust channel half power bandwidth
+ ctrl_pfir_subband.hp_factor = 1;
+ ctrl_pfir_subband.BWchan = ctrl_pfir_subband.hp_factor / tb.subband_fft_size; % Channel bandwidth
+ hfir_subband_coeff = pfir_coeff(ctrl_pfir_subband.nof_polyphases, ...
+ ctrl_pfir_subband.nof_taps, ...
+ ctrl_pfir_subband.BWchan, ...
+ ctrl_pfir_subband.r_pass, ...
+ ctrl_pfir_subband.r_stop, ...
+ ctrl_pfir_subband.coeff_w, ...
+ ctrl_pfir_subband.config);
+end
+ctrl_pfir_subband.coeff = reshape(hfir_subband_coeff, ctrl_pfir_subband.nof_polyphases, ctrl_pfir_subband.nof_taps);
+ctrl_pfir_subband.coeff = flipud(ctrl_pfir_subband.coeff);
+ctrl_pfir_subband.Zdelays = zeros(ctrl_pfir_subband.nof_polyphases, ctrl_pfir_subband.nof_taps-1);
+ctrl_pfir_subband.gain = 1; % no gain adjustment in PFIR, just apply the coefficients to the input data
+
+disp(sprintf(['Subband PFIR filter settings:\n', ...
+ '. Bypass = %d\n', ...
+ '. Number of polyphases = %d\n', ...
+ '. Number of taps = %d\n'], ...
+ ctrl_pfir_subband.bypass, ...
+ ctrl_pfir_subband.nof_polyphases, ...
+ ctrl_pfir_subband.nof_taps));
+if strcmp(tb.model_quantization, 'fixed point')
+ disp(sprintf(['. Coefficient width = %d\n', ...
+ '. Coefficient q_max = %d\n', ...
+ '. Data width = %d\n', ...
+ '. Data backoff width = %d\n', ...
+ ', Data q_full_scale = %d\n'], ...
+ ctrl_pfir_subband.coeff_w, ...
+ ctrl_pfir_subband.coeff_q_max, ...
+ ctrl_pfir_subband.data_w, ...
+ ctrl_pfir_subband.backoff_w, ...
+ ctrl_pfir_subband.data_q_full_scale));
+end
+
+%% Subband FFT parameters
+ctrl_pfft_subband.fft_size = tb.subband_fft_size;
+ctrl_pfft_subband.complex = false; % real input PFB
+ctrl_pfft_subband.gain = ctrl_pfft_subband.fft_size; % compensate default gain of FFT in PFFT
+if strcmp(tb.model_quantization, 'floating point')
+ ctrl_pfft_subband.data_w = 0;
+else
+ ctrl_pfft_subband.data_w = 16;
+end
+ctrl_pfft_subband.data_q_full_scale = 2^(ctrl_pfft_subband.data_w-1);
+
+
+%% Run the data path processing
+t_start = cputime;
+data_wg = zeros(tb.nof_tsub, tb.subband_fft_size);
+data_pfir_subband = zeros(tb.nof_tsub, tb.subband_fft_size);
+data_pfft_subband = zeros(tb.nof_tsub, tb.nof_subbands);
+for bi = 0:tb.nof_tsub-1
+ bI = bi+1;
+
+ % Data path (DP)
+ [ctrl_wg, block_wg] = wg( ctrl_wg);
+ [ctrl_pfir_subband, block_pfir_subband] = pfir( ctrl_pfir_subband, block_wg);
+ [ block_pfft_subband] = pfft( ctrl_pfft_subband, block_pfir_subband);
+
+ % Capture data at each DP interface
+ data_wg(bI, :) = block_wg;
+ data_pfir_subband(bI, :) = block_pfir_subband;
+ data_pfft_subband(bI, :) = block_pfft_subband;
+end
+t_stop = cputime;
+disp(sprintf('Total processing time: %f seconds\n', t_stop-t_start));
+
+
+%% Save coeffcients and data reference files for stimuli and verification of HDL implementation
+if strcmp(tb.model_quantization, 'fixed point')
+ file_name_prefix = 'data/run_pfir_m_';
+ L = ctrl_pfir_subband.nof_taps;
+ N = ctrl_pfir_subband.nof_polyphases; % = tb.subband_fft_size;
+
+ % Quantize subband PFIR filter coefficients
+ int_hfir = round(ctrl_pfir_subband.coeff_q_full_scale * hfir_subband_coeff); % map largest coefficient to q_max
+
+ % Quantize WG output
+ wg_q_data = round(ctrl_wg.q_full_scale * data_wg);
+
+ % Quantize PFIR subband output
+ pfir_subband_q_data = round(ctrl_pfir_subband.data_q_full_scale * data_pfir_subband);
+
+ % Quantize PFFT subband output real and imaginary
+ pfft_subband_q_data = round(ctrl_pfft_subband.data_q_full_scale * data_pfft_subband);
+
+ % Save the quantized integer subband PFIR filter coefficients in a separate file
+ file_name = sprintf([file_name_prefix, 'pfir_coeff_', ctrl_pfir_subband.config.design, '_%dtaps_%dpoints_%db.dat'], L, N, ctrl_pfir_subband.coeff_w);
+ fid = fopen(file_name, 'w');
+ fprintf(fid,'%d\n', int_hfir);
+ fclose(fid);
+
+ % Save PFIR filter coefficients again together with the data path signals for WG, PFIR and PFFT
+ tbegin = 1;
+ tend = tb.nof_tsub; % <= tb.nof_tsub
+ file_name = sprintf([file_name_prefix, '%s_%db_%dtaps_%dpoints_%db_%db.dat'], ctrl_wg.name, ctrl_wg.data_w, L, N, ctrl_pfir_subband.data_w, ctrl_pfft_subband.data_w);
+ fid = fopen(file_name, 'w');
+ fprintf(fid,'Nof lines PFIR coefficients: %d\n', length(int_hfir)); % save PFIR filter coefficients again also in the data file to have a complete set
+ fprintf(fid,'Nof lines WG output: %d\n', length(wg_q_data(:)));
+ fprintf(fid,'Nof lines PFIR output: %d\n', length(pfir_subband_q_data(:)));
+ fprintf(fid,'Nof lines PFFT output: %d\n', length(pfft_subband_q_data(:)));
+ % Subband PFIR coefficients
+ fprintf(fid,'%d\n', int_hfir);
+ % WG output
+ for bI = tbegin:tend
+ for bJ = 1:N
+ fprintf(fid,'%d\n', wg_q_data(bI, bJ));
+ end
+ end;
+ % PFIR subband output
+ for bI = tbegin:tend
+ for bJ = 1:N
+ fprintf(fid,'%d\n', pfir_subband_q_data(bI, bJ));
+ end
+ end
+ % PFFT subband output real and imaginary
+ for bI = tbegin:tend
+ for bJ = 1:tb.nof_subbands
+ fprintf(fid,'%10d%10d\n', real(pfft_subband_q_data(bI, bJ)), imag(pfft_subband_q_data(bI, bJ)));
+ end
+ end
+ fclose(fid);
+end
+
+
+%% Plot data path results
+xfig = 300;
+yfig = 200;
+xfigw = 1000;
+yfigw = 800;
+dfig = 20;
+
+ts = (0:tb.subband_fft_size*tb.nof_tsub-1)/tb.subband_fft_size; % time of ADC / WG samples in subband periods
+tsub_all = (0: tb.nof_subbands*tb.nof_tsub-1)/tb.nof_subbands; % time in subband periods for block of subbands
+tsub_one = (0: tb.nof_tsub-1); % time in subband periods for one subband
+
+sub_I = tb.subband_I + [0: tb.nof_subbands: tb.nof_subbands*tb.nof_tsub-1]; % get Matlab indices of all subband_I
+sub_Iplus1 = tb.subband_Iplus1 + [0: tb.nof_subbands: tb.nof_subbands*tb.nof_tsub-1]; % get Matlab indices of all next neighbour subband_I+1
+
+%% Plot WG output
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+data = data_wg';
+plot(ts, data(:))
+ylim([-1.3 1.3]);
+title(sprintf('WG output data (WG subband %6.3f)', tb.subband_wg));
+xlabel(sprintf('Time 0:%d [Tsub]', tb.nof_tsub-1));
+ylabel('Voltage');
+grid on;
+
+%% Plot FIR-filter coefficients
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+h = hfir_subband_coeff;
+NL = ctrl_pfir_subband.nof_coefficients;
+N = ctrl_pfir_subband.nof_polyphases;
+L = ctrl_pfir_subband.nof_taps;
+hi = 1:NL; % coefficients index
+hx = hi / N; % coefficients axis in tap units
+if ctrl_pfir_subband.coeff_w == 0
+ plot(hx, h);
+ title(['FIR filter coefficients for ', ctrl_pfir_subband.config.design, ' (floating point)']);
+else
+ plot(hx, h);
+ title(['FIR filter coefficients for ', ctrl_pfir_subband.config.design, ' (', num2str(ctrl_pfir_subband.coeff_w), ' bits)']);
+end
+ylo = min(h)*1.2;
+yhi = max(h)*1.2;
+ylim([ylo yhi]);
+grid on;
+
+%% Plot subband PFIR transfer function
+h = hfir_subband_coeff;
+NL = ctrl_pfir_subband.nof_coefficients;
+L = ctrl_pfir_subband.nof_taps;
+hf_abs = abs(fftshift(fft(h / sum(h), NL)));
+hi = 1:NL; % coefficients index
+fx = (hi - NL/2-1) / L; % frequency axis in subband units
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+plot(fx, db(hf_abs)); % db() = 20*log10 for voltage
+%xlim([-3 3]);
+grid on;
+title(['Subband PFIR filter transfer function for ', ctrl_pfir_subband.config.design, ' (nof taps = ', int2str(ctrl_pfir_subband.nof_taps), ')']);
+xlabel('Frequency [channels]');
+ylabel('Magnitude [dB]');
+
+%% Plot subband PFIR output
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+data = data_pfir_subband';
+plot(ts, data(:))
+title(sprintf('Subband PFIR filter output - FFT input data (WG subband %6.3f)', tb.subband_wg));
+ylim([-2 2]); % Delay tracking step when tb.subband_wg is .5 causes double range
+xlabel(sprintf('Time 0:%d [Tsub]', tb.nof_tsub-1));
+ylabel('Voltage');
+grid on;
+
+%% Plot PFFT subbands spectrum and phase for all tb.nof_tsub in one plot
+sub_ampl = abs(data_pfft_subband);
+sub_ampl_max = max(sub_ampl(:));
+sub_phase = angle(data_pfft_subband)*180/pi;
+x = sub_ampl < 0.1*sub_ampl_max;
+sub_phase(x) = 0; % force phase of too small signals to 0
+
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+subplot(2,1,1);
+data = sub_ampl';
+data = data(:);
+plot(tsub_all, data, 'k', tsub_all(sub_I), data(sub_I), 'ko', tsub_all(sub_Iplus1), data(sub_Iplus1), 'kx');
+title(sprintf('Subband data - amplitude (o,x = subband %d,%d for WG subband = %6.3f)', tb.subband_i, tb.subband_iplus1, tb.subband_wg));
+xlabel(sprintf('Subbands 0:%d at time 0:%d [Tsub]', tb.nof_subbands-1, tb.nof_tsub-1));
+ylabel('Voltage');
+grid on;
+subplot(2,1,2);
+data = sub_phase';
+data = data(:);
+plot(tsub_all, data, 'k', tsub_all(sub_I), data(sub_I), 'ko', tsub_all(sub_Iplus1), data(sub_Iplus1), 'kx');
+ylim([-180 180])
+title(sprintf('Subband data - amplitude (o,x = subband %d,%d for WG subband = %6.3f)', tb.subband_i, tb.subband_iplus1, tb.subband_wg));
+xlabel(sprintf('Subbands 0:%d at time 0:%d [Tsub]', tb.nof_subbands-1, tb.nof_tsub-1));
+ylabel('Phase [degrees]');
+grid on;
+
+%% Plot phase for the subband that is set in the WG
+wg_sub_phase = angle(data_pfft_subband(:, tb.subband_I))*180/pi;
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+plot(tsub_one, wg_sub_phase, '-o')
+ylim([-180 180])
+title(sprintf('Subband phase for subband %d in range 0:%d', tb.subband_i, tb.nof_subbands-1));
+xlabel(sprintf('Time 0:%d [Tsub]', tb.nof_tsub-1));
+ylabel('Phase [degrees]');
+grid on;
+
+%% Plot PFFT subbands spectrum and phase for all tb.nof_tsub in separate plots
+sub_ampl = abs(data_pfft_subband);
+sub_ampl_max = max(sub_ampl(:));
+sub_phase = angle(data_pfft_subband)*180/pi;
+x = sub_ampl < 0.1*sub_ampl_max;
+sub_phase(x) = 0; % force phase of too small signals to 0
+
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+
+if tb.plot_per_block
+ for bi = 0:tb.nof_tsub-1
+ bI = bi+1;
+ subplot(2,1,1);
+ data = sub_ampl(bI, :);
+ plot(0:tb.nof_subbands-1, data, 'k', tb.subband_i, data(tb.subband_I), 'ko', tb.subband_iplus1, data(tb.subband_Iplus1), 'kx');
+ ylim([0 sub_ampl_max])
+ title(sprintf('Subband data - amplitude (o,x = subband %d,%d for WG subband = %6.3f)', tb.subband_i, tb.subband_iplus1, tb.subband_wg));
+ xlabel(sprintf('Subbands 0:%d at time %d [Tsub]', tb.nof_subbands-1, bi));
+ ylabel('Voltage');
+ grid on;
+ subplot(2,1,2);
+ data = sub_phase(bI, :);
+ plot(0:tb.nof_subbands-1, data, 'k', tb.subband_i, data(tb.subband_I), 'ko', tb.subband_iplus1, data(tb.subband_Iplus1), 'kx');
+ ylim([-180 180])
+ title(sprintf('Subband data - amplitude (o,x = subband index %d,%d for WG subband = %6.3f)', tb.subband_i, tb.subband_iplus1, tb.subband_wg));
+ xlabel(sprintf('Subbands 0:%d at time %d [Tsub]', tb.nof_subbands-1, bi));
+ ylabel('Phase [degrees]');
+ grid on;
+ drawnow;
+ %pause(0.3)
+ end
+end
+
+%% Plot subband spectrogram for all tb.nof_tsub in one plot
+data = db(sub_ampl); % no need to scale data, range is already normalized
+if strcmp(tb.model_quantization, 'floating point')
+ db_low = -150;
+else
+ db_low = -20 - floor(6.02 * ctrl_pfft_subband.data_w);
+end
+%db_low = floor(min(data(data~=-Inf)))
+data(data<db_low) = db_low;
+mymap = jet(-db_low);
+%mymap(1,:) = [0 0 0]; % force black for dB(0)
+colormap(mymap);
+imagesc(data',[db_low 0]);
+colorbar;
+title(sprintf('Subband spectogram (max value = %f = %.2f dB)', sub_ampl_max, db(sub_ampl_max)));
+xlabel(sprintf('Time 0:%d [Tsub]', tb.nof_tsub-1));
+ylabel(sprintf('Subbands 0:%d', tb.nof_subbands-1));
+