diff --git a/applications/apertif/matlab/apertif_matlab_readme.txt b/applications/apertif/matlab/apertif_matlab_readme.txt
index 09b6fdc6eee4e81f0de9338bdc7b8e33c2f70d8c..c56eabf7760d0ef814c30ca98a0c717267f42fdc 100644
--- a/applications/apertif/matlab/apertif_matlab_readme.txt
+++ b/applications/apertif/matlab/apertif_matlab_readme.txt
@@ -36,15 +36,13 @@ This directory contains Matlab code for modelling DSP aspects of Apertif:
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 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
sinusoids_and_fft_frequency_bins.m Script Code on DFT and DTFT from Matlab blog on internet
try_fft.m Script Try DFT and DTFT
-5) Data path HDL functions
+5) Data path functions
wg.m Function Waveform generator per block of data
dt.m Function Delay tracking per block of data
bsn_source.m Function Block Sequence Number source per block of data
@@ -54,10 +52,6 @@ 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 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
one_pfb.m Script Run data path model with subband polyphase filterbank (real input)
@@ -66,3 +60,11 @@ This directory contains Matlab code for modelling DSP aspects of Apertif:
two_pfb.m Script Run data path model with subband polyphase filterbank followed by a channel polyphase filterbank
+7) HDL reference data
+ run_pfir.m Script Run the PFIR with real input to create HDL reference data (based on one_pfb.m)
+ 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)
+ run_pfb_complex.m Script Run the PFB with complex input to create HDL reference data (based on one_pfb.m)
+
+
+
diff --git a/applications/apertif/matlab/run_pfb_complex.m b/applications/apertif/matlab/run_pfb_complex.m
new file mode 100644
index 0000000000000000000000000000000000000000..b9cb265e0c23186baca2a9c1de0b4160e6073a67
--- /dev/null
+++ b/applications/apertif/matlab/run_pfb_complex.m
@@ -0,0 +1,488 @@
+%-----------------------------------------------------------------------------
+%
+% 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 PFB = PFIR + PFFT with complex input:
+% To create input and output reference data for HDL implementation in
+% data/, e.g.:
+%
+% . PFIR coeffcients: run_pfb_complex_m_pfir_coeff_fircls1_16taps_128points_16b.dat
+% . PFIR coeffcients, WG, PFIR, PFFT data: run_pfb_complex_m_chirp_8b_16taps_128points_16b_16b.dat
+%
+% Description :
+%
+% * See one_pfb.m for general description.
+% This run_pfb_complex.m is almost the same as run_pfft_complex.m,
+% except that it includes the PFIR.
+
+clear all;
+close all;
+rng('default');
+fig=0;
+
+%% Selectable test bench settings
+% Input signal
+tb.model_signal = 'sinusoid'; % Use sinusoid to check the frequency response of the FFT with frequency set via tb.channel_wg
+tb.model_signal = 'phasor'; % Use phasor to check the frequency response of the FFT with frequency set via tb.channel_wg
+%tb.model_signal = 'noise'; % Use noise to check all aspects in one signal
+%tb.model_signal = 'square'; % Use square to create a square wave with period set via tb.channel_wg
+%tb.model_signal = 'impulse'; % Use impulse to check the impulse response for different ctrl_wg.sop
+
+tb.model_quantization = 'floating point';
+tb.model_quantization = 'fixed point';
+tb.nof_channels = 32;
+
+% Carrier frequency
+tb.channel_wg = 1; % channel range -tb.nof_channels/2:tb.nof_channels/2-1, can be fraction to have any sinusoid frequency
+%tb.channel_wg = 1.55;
+%tb.channel_wg = 12;
+
+tb.phase = 0; % initial 'sinusoid' phase offset normalized to 2pi
+tb.sop = 1; % initial 'impulse' start index in range ctrl_wg.block_size = tb.channel_fft_size
+
+% Model a frequency sweep of the 'sinusoid'
+tb.chirp = 0; % 0 = use fixed tb.channel_wg frequency or pulse period equal to block_size
+tb.chirp = 1; % else increment WG frequency every block to have chirp frequency sweep or slide the pulse
+if strcmp(tb.model_signal, 'noise')
+ tb.nof_tchan = 10 ; % number of channel periods to simulate
+elseif tb.chirp
+ tb.nof_tchan = 200; % number of channel periods to simulate
+else
+ tb.nof_tchan = 5; % number of channel 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.channel_i = floor(tb.channel_wg + tb.nof_channels/2); % natural channel index in range 0:tb.nof_channels-1
+tb.channel_I = tb.channel_i+1; % equivalent Matlab indexing for range 1:tb.nof_channels
+if tb.channel_i<tb.nof_channels-1 % determine index of next neighbour channel
+ tb.channel_iplus1 = tb.channel_i+1;
+else
+ tb.channel_iplus1 = 0; % wrap channel index tb.nof_channels to channel index 0
+end
+tb.channel_Iplus1 = tb.channel_iplus1+1; % equivalent Matlab index for next neighbour channel in range 1:tb.nof_channels
+tb.channel_fft_size = tb.nof_channels; % channel complex FFT
+
+fs = 1; % normalized sample frequency
+fchan = fs/tb.channel_fft_size; % channel frequency relative to fs
+
+disp(sprintf (['Test bench settings:\n', ...
+ '. Model coefficients and data = %s\n', ...
+ '. Number of channels = %d\n', ...
+ '. Number of channel periods to simulate = %d\n', ...
+ '. Channel filterbank complex FFT size = %d\n', ...
+ '. WG input signal = %s\n', ...
+ '. WG channel frequency = %f\n'], ...
+ tb.model_quantization, tb.nof_channels, tb.nof_tchan, tb.channel_fft_size, tb.model_signal, tb.channel_wg))
+
+
+%% Waveform generator
+% Common WG settings
+ctrl_wg.complex = false;
+if strcmp(tb.model_signal, 'phasor') || strcmp(tb.model_signal, 'noise')
+ ctrl_wg.complex = true;
+end
+ctrl_wg.block_nr = 0;
+ctrl_wg.block_size = tb.channel_fft_size;
+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 in floating point
+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.5*lsb;
+ctrl_wg.agwn_sigma = 0; % AGWN sigma
+ctrl_wg.ampl = 1; % full scale is 1, 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;
+ %ctrl_wg.ampl = 1*lsb;
+end
+ctrl_wg.offset = 0; % DC offset
+%ctrl_wg.offset = 0.1;
+
+% Signal specific WG settings
+if strcmp(tb.model_signal, 'sinusoid') || strcmp(tb.model_signal, 'phasor')
+ ctrl_wg.signal = 'sinusoid';
+ ctrl_wg.psk = 0; % no PSK modulation
+ if tb.chirp
+ ctrl_wg.df = 0.49*fchan; % increment freq by df per block to create chirp
+ %ctrl_wg.df = 0.5*fchan;
+ else
+ ctrl_wg.df = 0;
+ end
+ if ctrl_wg.df == 0
+ ctrl_wg.freq = tb.channel_wg*fchan; % start WG at tb.channel_wg
+ else
+ ctrl_wg.freq = (tb.channel_i-1)*fchan; % start WG about one channel before tb.channel_wg to account for impulse response of several tchan
+ end
+ ctrl_wg.phase = tb.phase; % Phase offset normalized to 2pi
+ if ctrl_wg.freq == 0 && ctrl_wg.df == 0
+ ctrl_wg.offset = 1; % Force DC offset for zero Hz
+ end
+elseif strcmp(tb.model_signal, 'noise')
+ ctrl_wg.signal = 'noise';
+ ctrl_wg.ampl = 0;
+ ctrl_wg.agwn_sigma = 1/3; % AGWN sigma is 1/3 full scale
+else
+ ctrl_wg.signal = 'pulse';
+ ctrl_wg.sop = tb.sop; % Initial start of pulse index normalized to ctrl_wg.block_size = tb.channel_fft_size
+ ctrl_wg.period = ctrl_wg.block_size; % Pulse period
+ if tb.chirp
+ ctrl_wg.period = ctrl_wg.block_size+1; % Pulse period for pulse moving in block like a sop chirp
+ end
+ if strcmp(tb.model_signal, 'impulse')
+ ctrl_wg.width = 1; % Pulse width for FFT impulse response
+ else % default to 'square')
+ ctrl_wg.ampl = 2*ctrl_wg.ampl;
+ ctrl_wg.offset = -ctrl_wg.ampl/2;
+ ctrl_wg.period = ctrl_wg.block_size/tb.channel_wg; % Pulse period equal to channel_wg period
+ ctrl_wg.width = ctrl_wg.period/2; % Pulse width for square wave
+ end
+ ctrl_wg.width_remaining = 0;
+end
+
+%% Channel PFIR filter parameters
+ctrl_pfir_channel.bypass = 0;
+ctrl_pfir_channel.nof_polyphases = tb.channel_fft_size;
+
+if strcmp(tb.model_quantization, 'floating point')
+ ctrl_pfir_channel.data_w = 0; % no quantization
+ ctrl_pfir_channel.coeff_w = 0; % no quantization
+else
+ ctrl_pfir_channel.data_w = 16;
+ ctrl_pfir_channel.backoff_w = 1; % attenuate PFIR data output by 2^backoff_w before quantization to account for PFIR overshoot
+ ctrl_pfir_channel.coeff_w = 16; % nof bits per coefficient including 1 sign bit
+ ctrl_pfir_channel.coeff_q_max = 2^(ctrl_pfir_channel.coeff_w-1)-1;
+ ctrl_pfir_channel.coeff_q_full_scale = 2^(ctrl_pfir_channel.coeff_w-1);
+ ctrl_pfir_channel.data_q_full_scale = 2^(ctrl_pfir_channel.data_w-1-ctrl_pfir_channel.backoff_w);
+end
+
+% Calculate free choice PFIR coefficients (can be floating point or fixed point)
+ctrl_pfir_channel.nof_taps = 16; % Number of taps
+ctrl_pfir_channel.nof_coefficients = ctrl_pfir_channel.nof_polyphases*ctrl_pfir_channel.nof_taps; % Number of filter coefficients (taps)
+ctrl_pfir_channel.config.design = 'fir1';
+ctrl_pfir_channel.config.design = 'fircls1'; % 'fir1', 'fircls1'
+ctrl_pfir_channel.config.design_flag = ''; % only for fircls1: 'trace', ''
+ctrl_pfir_channel.config.interpolate = 'interpft'; % only for fircls1: 'resample', 'fourier', 'interpft'
+ctrl_pfir_channel.r_pass = 1e-3; % only for fircls1
+ctrl_pfir_channel.r_stop = 1e-4; % only for fircls1
+ctrl_pfir_channel.hp_factor = 1.050; % Adjust channel half power bandwidth
+ctrl_pfir_channel.hp_factor = 1;
+ctrl_pfir_channel.BWchan = ctrl_pfir_channel.hp_factor / tb.channel_fft_size; % Channel bandwidth
+hfir_channel_coeff = pfir_coeff(ctrl_pfir_channel.nof_polyphases, ...
+ ctrl_pfir_channel.nof_taps, ...
+ ctrl_pfir_channel.BWchan, ...
+ ctrl_pfir_channel.r_pass, ...
+ ctrl_pfir_channel.r_stop, ...
+ ctrl_pfir_channel.coeff_w, ...
+ ctrl_pfir_channel.config);
+ctrl_pfir_channel.coeff = reshape(hfir_channel_coeff, ctrl_pfir_channel.nof_polyphases, ctrl_pfir_channel.nof_taps);
+ctrl_pfir_channel.coeff = flipud(ctrl_pfir_channel.coeff);
+ctrl_pfir_channel.Zdelays = zeros(ctrl_pfir_channel.nof_polyphases, ctrl_pfir_channel.nof_taps-1);
+ctrl_pfir_channel.gain = 1; % no gain adjustment in PFIR, just apply the coefficients to the input data
+
+disp(sprintf(['Channel PFIR filter settings:\n', ...
+ '. Bypass = %d\n', ...
+ '. Number of polyphases = %d\n', ...
+ '. Number of taps = %d\n'], ...
+ ctrl_pfir_channel.bypass, ...
+ ctrl_pfir_channel.nof_polyphases, ...
+ ctrl_pfir_channel.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_channel.coeff_w, ...
+ ctrl_pfir_channel.coeff_q_max, ...
+ ctrl_pfir_channel.data_w, ...
+ ctrl_pfir_channel.backoff_w, ...
+ ctrl_pfir_channel.data_q_full_scale));
+end
+
+%% Channel FFT parameters
+ctrl_pfft_channel.fft_size = tb.channel_fft_size;
+ctrl_pfft_channel.complex = true; % complex input FFT
+ctrl_pfft_channel.gain = ctrl_pfft_channel.fft_size;
+if strcmp(tb.model_quantization, 'floating point')
+ ctrl_pfft_channel.data_w = 0;
+else
+ ctrl_pfft_channel.data_w = 16;
+end
+ctrl_pfft_channel.data_q_full_scale = 2^(ctrl_pfft_channel.data_w-1);
+
+
+%% Run the data path processing
+t_start = cputime;
+data_wg = zeros(tb.nof_tchan, tb.channel_fft_size);
+data_pfir_channel = zeros(tb.nof_tchan, tb.channel_fft_size);
+data_pfft_channel = zeros(tb.nof_tchan, tb.nof_channels);
+for bi = 0:tb.nof_tchan-1
+ bI = bi+1;
+
+ % Data path (DP)
+ [ctrl_wg, block_wg] = wg( ctrl_wg);
+ [ctrl_pfir_channel, block_pfir_channel] = pfir( ctrl_pfir_channel, block_wg);
+ [ block_pfft_channel] = pfft( ctrl_pfft_channel, block_pfir_channel);
+
+ % Capture data at each DP interface
+ data_wg(bI, :) = block_wg;
+ data_pfir_channel(bI, :) = block_pfir_channel;
+ data_pfft_channel(bI, :) = block_pfft_channel;
+end
+t_stop = cputime;
+disp(sprintf('Total processing time: %f seconds', t_stop-t_start));
+
+
+%% Save coeffcients and data reference files for stimuli and verification of HDL implementation
+if strcmp(tb.model_quantization, 'fixed point')
+ % Quantize channel PFIR filter coefficients
+ int_hfir = round(ctrl_pfir_channel.coeff_q_full_scale * hfir_channel_coeff); % map largest coefficient to q_max
+
+ % Quantize WG output
+ wg_q_data = round(ctrl_wg.q_full_scale * data_wg);
+
+ % Quantize PFIR channel output
+ pfir_channel_q_data = round(ctrl_pfir_channel.data_q_full_scale * data_pfir_channel);
+
+ % Quantize PFFT channel output real and imaginary
+ pfft_channel_q_data = round(ctrl_pfft_channel.data_q_full_scale * data_pfft_channel);
+
+ % Save the quantized integer channel PFIR filter coefficients in a separate file
+ L = ctrl_pfir_channel.nof_taps;
+ N = ctrl_pfir_channel.nof_polyphases; % = tb.channel_fft_size;
+ file_name_prefix = 'data/run_pfb_complex_m_';
+ file_name = sprintf([file_name_prefix, 'pfir_coeff_', ctrl_pfir_channel.config.design, '_%dtaps_%dpoints_%db.dat'], L, N, ctrl_pfir_channel.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_tchan; % <= tb.nof_tchan
+ file_name = sprintf([file_name_prefix, '%s_%db_%dtaps_%dpoints_%db_%db.dat'], ctrl_wg.name, ctrl_wg.data_w, L, N, ctrl_pfir_channel.data_w, ctrl_pfft_channel.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_channel_q_data(:)));
+ fprintf(fid,'Nof lines PFFT output: %d\n', length(pfft_channel_q_data(:)));
+ % Channel PFIR coefficients
+ fprintf(fid,'%d\n', int_hfir);
+ % WG output
+ for bI = tbegin:tend
+ for bJ = 1:N
+ fprintf(fid,'%8d%8d\n', real(wg_q_data(bI, bJ)), imag(wg_q_data(bI, bJ)));
+ end
+ end;
+ % PFIR channel output
+ for bI = tbegin:tend
+ for bJ = 1:N
+ fprintf(fid,'%12d%12d\n', real(pfir_channel_q_data(bI, bJ)), imag(pfir_channel_q_data(bI, bJ)));
+ end
+ end
+ % PFFT channel output real and imaginary
+ for bI = tbegin:tend
+ for bJ = 1:tb.nof_channels
+ fprintf(fid,'%10d%10d\n', real(pfft_channel_q_data(bI, bJ)), imag(pfft_channel_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.channel_fft_size*tb.nof_tchan-1)/tb.channel_fft_size; % time of ADC / WG samples in channel periods
+tchan_all = (0: tb.nof_channels*tb.nof_tchan-1)/tb.nof_channels; % time in channel periods for block of channels
+tchan_one = (0: tb.nof_tchan-1); % time in channel periods for one channel
+
+chan_I = tb.channel_I + [0: tb.nof_channels: tb.nof_channels*tb.nof_tchan-1]; % get Matlab indices of all channel_I
+chan_Iplus1 = tb.channel_Iplus1 + [0: tb.nof_channels: tb.nof_channels*tb.nof_tchan-1]; % get Matlab indices of all next neighbour channel_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, real(data(:)), 'r', ts, imag(data(:)), 'b')
+ylim([-1.3 1.3]);
+title(sprintf('WG output data (WG channel %6.3f)', tb.channel_wg));
+xlabel(sprintf('Time 0:%d [Tchan]', tb.nof_tchan-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_channel_coeff;
+NL = ctrl_pfir_channel.nof_coefficients;
+N = ctrl_pfir_channel.nof_polyphases;
+L = ctrl_pfir_channel.nof_taps;
+hi = 1:NL; % coefficients index
+hx = hi / N; % coefficients axis in tap units
+if ctrl_pfir_channel.coeff_w == 0
+ plot(hx, h);
+ title(['FIR filter coefficients for ', ctrl_pfir_channel.config.design, ' (floating point)']);
+else
+ plot(hx, h);
+ title(['FIR filter coefficients for ', ctrl_pfir_channel.config.design, ' (', num2str(ctrl_pfir_channel.coeff_w), ' bits)']);
+end
+ylo = min(h)*1.2;
+yhi = max(h)*1.2;
+ylim([ylo yhi]);
+grid on;
+
+%% Plot channel PFIR transfer function
+h = hfir_channel_coeff;
+NL = ctrl_pfir_channel.nof_coefficients;
+L = ctrl_pfir_channel.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 channel 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(['Channel PFIR filter transfer function for ', ctrl_pfir_channel.config.design, ' (nof taps = ', int2str(ctrl_pfir_channel.nof_taps), ')']);
+xlabel('Frequency [channels]');
+ylabel('Magnitude [dB]');
+
+%% Plot channel PFIR output
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+data = data_pfir_channel.';
+plot(ts, real(data(:)), 'r', ts, imag(data(:)), 'b')
+title(sprintf('Channel PFIR filter output - FFT input data (WG channel %6.3f)', tb.channel_wg));
+ylim([-2 2]); % Delay tracking step when tb.channel_wg is .5 causes double range
+xlabel(sprintf('Time 0:%d [Tchan]', tb.nof_tchan-1));
+ylabel('Voltage');
+grid on;
+
+%% Plot PFFT channels spectrum and phase for all tb.nof_tchan in one plot
+chan_ampl = abs(data_pfft_channel);
+chan_ampl_max = max(chan_ampl(:));
+chan_phase = angle(data_pfft_channel)*180/pi;
+x = chan_ampl < 0.1*chan_ampl_max;
+chan_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 = chan_ampl.';
+data = data(:);
+plot(tchan_all, data, 'k', tchan_all(chan_I), data(chan_I), 'ko', tchan_all(chan_Iplus1), data(chan_Iplus1), 'kx');
+title(sprintf('Channel data - amplitude (o,x = channel %d,%d for WG channel = %6.3f)', tb.channel_i, tb.channel_iplus1, tb.channel_wg));
+xlabel(sprintf('Channels 0:%d at time 0:%d [Tchan]', tb.nof_channels-1, tb.nof_tchan-1));
+ylabel('Voltage');
+grid on;
+subplot(2,1,2);
+data = chan_phase.';
+data = data(:);
+plot(tchan_all, data, 'k', tchan_all(chan_I), data(chan_I), 'ko', tchan_all(chan_Iplus1), data(chan_Iplus1), 'kx');
+ylim([-180 180])
+title(sprintf('Channel data - amplitude (o,x = channel %d,%d for WG channel = %6.3f)', tb.channel_i, tb.channel_iplus1, tb.channel_wg));
+xlabel(sprintf('Channels 0:%d at time 0:%d [Tchan]', tb.nof_channels-1, tb.nof_tchan-1));
+ylabel('Phase [degrees]');
+grid on;
+
+%% Plot phase for the channel that is set in the WG
+wg_chan_phase = angle(data_pfft_channel(:, tb.channel_I))*180/pi;
+fig=fig+1;
+figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
+figure(fig);
+plot(tchan_one, wg_chan_phase, '-o')
+ylim([-180 180])
+title(sprintf('Channel phase for channel %d in range 0:%d', tb.channel_i, tb.nof_channels-1));
+xlabel(sprintf('Time 0:%d [Tchan]', tb.nof_tchan-1));
+ylabel('Phase [degrees]');
+grid on;
+
+%% Plot PFFT channels spectrum and phase for all tb.nof_tchan in separate plots
+chan_ampl = abs(data_pfft_channel);
+chan_ampl_max = max(chan_ampl(:));
+chan_phase = angle(data_pfft_channel)*180/pi;
+x = chan_ampl < 0.1*chan_ampl_max;
+chan_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_tchan-1
+ bI = bi+1;
+ subplot(2,1,1);
+ data = chan_ampl(bI, :);
+ plot(0:tb.nof_channels-1, data, 'k', tb.channel_i, data(tb.channel_I), 'ko', tb.channel_iplus1, data(tb.channel_Iplus1), 'kx');
+ ylim([0 chan_ampl_max])
+ title(sprintf('Channel data - amplitude (o,x = channel %d,%d for WG channel = %6.3f)', tb.channel_i, tb.channel_iplus1, tb.channel_wg));
+ xlabel(sprintf('Channels 0:%d at time %d [Tchan]', tb.nof_channels-1, bi));
+ ylabel('Voltage');
+ grid on;
+ subplot(2,1,2);
+ data = chan_phase(bI, :);
+ plot(0:tb.nof_channels-1, data, 'k', tb.channel_i, data(tb.channel_I), 'ko', tb.channel_iplus1, data(tb.channel_Iplus1), 'kx');
+ ylim([-180 180])
+ title(sprintf('Channel data - amplitude (o,x = channel index %d,%d for WG channel = %6.3f)', tb.channel_i, tb.channel_iplus1, tb.channel_wg));
+ xlabel(sprintf('Channels 0:%d at time %d [Tchan]', tb.nof_channels-1, bi));
+ ylabel('Phase [degrees]');
+ grid on;
+ drawnow;
+ %pause(0.3)
+ end
+end
+
+%% Plot channel spectrogram for all tb.nof_tchan in one plot
+data = db(chan_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_channel.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('Channel spectogram (max value = %f = %.2f dB)', chan_ampl_max, db(chan_ampl_max)));
+xlabel(sprintf('Time 0:%d [Tchan]', tb.nof_tchan-1));
+ylabel(sprintf('Channels 0:%d', tb.nof_channels-1));
+