Select Git revision
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
one_pfb.m 29.89 KiB
%-----------------------------------------------------------------------------
%
% 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: W. Lubberhuizen, 2015 (top.m for LOFAR Station, uses classes)
% Eric Kooistra, 2016 (for Apertif)
%
% Purpose : Model one PFB = PFIR + PFFT with real input:
% 1) to show the transfer function of a polyphase filterbank with real
% input like the subband filterbank
% 2) to create reference data for HDL implementation in data/:
%
% . PFIR coeffcients: one_pfb_m_pfir_coeff_fircls1_16taps_128points_16b.dat
% . PFIR coeffcients, WG, PFIR, PFFT data: one_pfb_m_chirp_8b_16taps_128points_16b_16b.dat
%
% Description :
%
% * General
% The data path (DP) is modelled per block of data. The block of data
% are counted by the block sequence number that thus acts as a
% timestamp. The data path consist of DSP functions. Each DSP
% function expects normalized full scale input (range -1:1) and it also
% normalized full scale output. Dependend on the tb.model_quantization
% the DSP functions outputs floating point or fixed point data that is
% quantized to the specified number of bits. The output full scale
% range remains -1:1.
%
% * Subband filterbank transfer function
% The real input implies that ctrl_pfft_subband.complex = false, and
% that the filterbank outputs tb.nof_subbands, where tb.subband_fft_size
% = tb.nof_complex*tb.nof_subbands.
% The subband filterbank can use LOFAR settings (tb.model_filterbank =
% 'LOFAR' or 'free choice' settings. With 'free choice' the number of
% subbands (tb.nof_subbands) and the quantization (tb.model_quantization)
% are free to set.
% Using a chirp via tb.chirp = 1 and ctrl_wg.df > 0 shows the transfer
% function of the polyphase filterbank.
% If tb.plot_per_block == 1 then the filterbank subband spectrum is
% plotted per block, else it is only plotted as a subband spectrogram
% image. Used tb.plot_per_block == 1 to debug the indexing in the
% filterbank input, pfir and pfft chain. Especially the
% flipud(ctrl_pfir_subband.coeff) is essential to have the indexing
% correct (see pfir.m).
% The chirp frequency increments forever, so the filterbank output
% sweeps from subband 0 up to fs/2, down to 0, up to fs/2, etc. For
% chirp frequency between DC and fsub/2 the subband 0 output is irratic
% between 1 and 0. For chirp frequency between fs/2-fsub/2 and fs/2 the
% the subband 0 output is irratic but almost 0.
% The PFIR can be bypassed using ctrl_pfir_subband.bypass = 0 to only
% have the FFT in the filterbank. Without PFIR and frequency there is
% aliasing on all subbands and no particular difference for frequencies
% near DC or fs/2.
% Towards fs/2 the WG sinusiod amplitude gets modulated more and more
% due to that there are less and less samples per period. Hence within% the impulse response of the filterbank the WG sinusoid is not stable
% for higher frequencies.
%
% * Indexing
% Matlab indices are from 1:N, but most languages use 0:N-1. Typically
% time starts at 0 and bin 0 represents DC. Therefore here default use
% 0:N-1 and add +1 to the index whenever it is used as index in a Matlab
% vector. Use capital I in index name if it is already a Matlab index.
%
% * Block processing
% The data is processed in blocks and captured at each interface in a
% matrix of tb.nof_tsub rows. The block size before the subband
% filterbank is defined by the FFT size (tb.subband_fft_size). The block
% size after the subband filterbank is defined by the number of subbands
% (tb.nof_subbands).
%
% * Flipped order of FIR coefficients:
% The FIR filter in the PFB is defined by a prototype filter h[] with N_fft
% * N_taps = 1024 *16 = 16384 coefficients. The impulse response of this
% prototype filter will show the coefficients in h[0:16383] order. The
% FIR coefficients of h[] in the PFB are flipped per column of N_fft = 1024
% coefficients, for all N_taps = 16 columns. In VHDL simulation this shows
% with tb_verify_pfb_response.vhd and in the Matlab model this flip is also
% done by ctrl_pfir_subband.coeff = flipud(ctrl_pfir_subband.coeff).
% With the flip the one_pfb.m yields SNR subband / spurious max = 64.305310
% dB in subband 64 for an input WG at 64.015625. Without the flip the SNR
% subband / spurious max = 29.857938 dB, so much worse. Hence flipping the
% FIR coefficients is needed, and not due to an implementation detail in
% the VHDL. The run_pfb.m is equivalent to one_pfb.m and is used to create
% reference input and expected output data for the wpfb implementation in
% VHDL by wpfb_unit_wide.vhd. The tb_tb_wpfb_unit_wide.vhd verifies multiple
% sets of reference data including sinus, chirp, noise and real input or
% complex input and wideband (f_sample > f_clk) or same rate (f_sample =
% f_clk). This tb guarantees that the VHDL agrees with the Matlab model.
% In the Matlab model the data blocks and data time are defined as:
%
% WG FIR with four taps FFT
% index t t
% 0 -1023 1023 2047 3071 4095 --> + --> -1023
% . . . . .
% . . . . .
% -2 2 . . . -2
% -1 1 . . . -1
% 1023 0 0 1024 2048 3072 --> + --> 0
% d h h h h d
%
% The waveform generator (WG) generates a block of 1024 samples with index
% 0:1023, where sample at index 1023 is the newest and sample at index 0 is
% the oldest. In a filter the newest sample needs to be multiplied with h[0]
% and the older samples are multiplied by the subsequent coefficients,
% because it is a convolution. Therefore the FIR coefficients need to be
% flipped up/down per column, to allow doing the filter as a d .* h vector
% multiply in pfir.m. The FFT operates on blocks of data with same index and
% time range as the WG. The data output of the FIR filter fits this input
% range of the FFT. Therefore no data flipping is needed.
%
% In the tb_verify_pfb_response.vhd the input stimuli is a block of N_fft =
% 1024 ones followed by N_taps-1 blocks with zeros. In time the oldest data
% will appear first in the simulator Wave Window, so therefore the
% fil_re_scope signal in the tb_verify_pfb_response.vhd will show the FIR
% coefficients h[] in order 1023:0, 2047:1024, ..., so flipped per block.
%
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.4; % subband range 0:tb.nof_subbands-1, can be fraction to have any sinusoid frequency
%tb.subband_wg = 1.55;
%tb.subband_wg = 12;
tb.subband_wg = 32;
tb.subband_wg = 64 + 1/64;
% 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.overflow = 'clip';
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;
ctrl_wg.q_full_scale = 2^(ctrl_wg.data_w-1);
lsb = 1/2^ctrl_wg.data_w;
end
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.overshoot_w = 1; % account for factor 2^overshoot_w in PFIR output quantization
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.overshoot_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');
hfir_subband_coeff = hfir_subband_coeff/max(hfir_subband_coeff);
hfir_subband_coeff = hfir_subband_coeff'; % Use column vector, same format as by pfir_coeff()
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
% Try effect if FIR coefficients would be shiftedsh = 1024;
sh = 15*1024;
sh = 0;
%hfir_subband_coeff = [hfir_subband_coeff(1+sh:16384), hfir_subband_coeff(1:sh)];
% Try effect if FIR coefficients would be in reversed order per wideband
% (ro = 4), per polyphase (ro = 1024) or another factor. There occur peaks
% at subband index N_fft / ro, so for ro = 4 at subband +-256+-0.5 and
% +-511.5. For ro = 32 at subband +-32+-0.5, +-64+-0.5, etc. Using
% ro = 1024 is equivalent to reverse order per polyphase and this yields
% peaks at +-0.5, +-1.5+-0.25, +-2.5+-0.25, +-3.5+-0.5, etc.
ro = 16;
x = reshape(hfir_subband_coeff, ro, ctrl_pfir_subband.nof_polyphases*ctrl_pfir_subband.nof_taps/ro);
x = flipud(x);
%hfir_subband_coeff = x(:);
% Try effect if FIR coefficients would be in reversed per tap
x = reshape(hfir_subband_coeff, ctrl_pfir_subband.nof_polyphases, ctrl_pfir_subband.nof_taps);
x = fliplr(x);
%hfir_subband_coeff = x(:);
% Flip ctrl_pfir_subband.coeff per poly phase, because that is the order
% in which the pfir() model and HDL expect the FIR coefficients
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 full scale 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.overshoot_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.full_scale_w = 0;
ctrl_pfft_subband.data_w = 16;
ctrl_pfft_subband.data_q_full_scale = 2^(ctrl_pfft_subband.data_w-1);
end
%% 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_pfir_subband = fliplr(block_pfir_subband); % Try flipped FFT input
[ 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/one_pfb_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
%% Log subband statistics for subband_i
% bin indices for a subband spectrum after the impulse response has
% settled, so after ctrl_pfir_subband.nof_taps * tb.nof_subbands
tg = (ctrl_pfir_subband.nof_taps + 4) * tb.nof_subbands;
tgI = tg + 1;
sub_bins = tgI + [0:tb.nof_subbands-1]; % Matlab indices of subbands in block tg
sub_I_bin = tgI + tb.subband_i; % Matlab index of subband_I
sub_ampl = abs(data_pfft_subband);
data = sub_ampl.';
data = data(:);
sub_I_ampl = data(sub_I_bin);
sub_I_power = sub_I_ampl^2;
sub_spurious = data(sub_bins);
sub_spurious(tb.subband_I) = 0;
sub_spurious_max = max(sub_spurious);
sub_spurious_power = sum(sub_spurious.^2);
sub_I_snr = sub_I_power / sub_spurious_power;
sub_I_snr_max = sub_I_power / sub_spurious_max^2;
disp(sprintf(['\nStatistics for WG subband %f in subband bin %d:\n', ...
'. Subband amplitude = %10.8f\n', ...
'. Subband power = %10.8f\n', ...
'. Spurious amplitude max = %10.8f\n', ...
'. Spurious total power = %10.8f\n', ...
'. SNR subband / spurious total = %f dB\n', ...
'. SNR subband / spurious max = %f dB\n'], ...
tb.subband_wg, tb.subband_i, ...
sub_I_ampl, ...
sub_I_power, ...
sub_spurious_max, ...
sub_spurious_power, ...
10 * log10(sub_I_snr), ...
10 * log10(sub_I_snr_max)));
%% 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 indexhx = 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([-3 3]); % 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(3,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;
% Plot subrange of Tsub after impulse response at low level scale to zoom
% in on spurious levels, see also plot subband spectrogram
subplot(3,1,2);
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');
xlim([19.9 23.1])
ylim([-0.001 0.005])title(sprintf('Subband data - spurious levels (for WG subband = %6.3f)', 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(3,1,3);
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 - phase (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));