279 lines
4.7 KiB
Matlab
279 lines
4.7 KiB
Matlab
%% Test fracF_dpw using a physical LFM
|
|
%
|
|
% Parameters chosen to match previous FrFT validation work:
|
|
%
|
|
% Fs = 512 MHz
|
|
% T = 1 us
|
|
% B = 64 MHz
|
|
%
|
|
% Matched FrFT order:
|
|
%
|
|
% a = -(2/pi)*atan(Fs/(beta*T))
|
|
%
|
|
% where:
|
|
%
|
|
% beta = B/T
|
|
%
|
|
% Notes:
|
|
% - FFT is computed on the original (non-interpolated) signal.
|
|
% - FrFT is computed on the interpolated signal.
|
|
% - Power spectra are averaged across the entire DPW.
|
|
|
|
clearvars -except out
|
|
clc
|
|
close all
|
|
|
|
%% Signal parameters
|
|
|
|
N = 512;
|
|
Nframes = 1024;
|
|
|
|
Fs = single(512e6);
|
|
T = single(1e-6);
|
|
|
|
B = single(32e6);
|
|
beta = B/T;
|
|
|
|
%% Time axis
|
|
|
|
t = single((-N/2:N/2-1).') / Fs;
|
|
|
|
%% Generate LFM
|
|
|
|
x = exp(1j*pi*beta*(t.^2));
|
|
|
|
x = complex(single(real(x)), ...
|
|
single(imag(x)));
|
|
|
|
%% Create DPW
|
|
|
|
X = repmat(x,1,Nframes);
|
|
|
|
%% Interpolate exactly as Simulink
|
|
|
|
halfbandInterp = dsp.FIRHalfbandInterpolator;
|
|
|
|
Xint = halfbandInterp(X);
|
|
|
|
%% Matched FrFT order
|
|
|
|
aMatch = single(-(2/pi)*atan(Fs/(beta*T)));
|
|
|
|
fprintf('\n');
|
|
fprintf('Fs = %.3f MHz\n',double(Fs)/1e6);
|
|
fprintf('T = %.3f us\n',double(T)*1e6);
|
|
fprintf('B = %.3f MHz\n',double(B)/1e6);
|
|
fprintf('aMatch = %.6f\n',double(aMatch));
|
|
|
|
%% FFT reference
|
|
%
|
|
% FFT detector operates on the original non-interpolated signal.
|
|
|
|
FFTref = fftshift(fft(X,[],1),1)/N;
|
|
|
|
%% FrFT
|
|
%
|
|
% FrFT detector operates on the interpolated signal.
|
|
|
|
[Achirp,H,Cchirp,Aa] = fracF_init(aMatch);
|
|
|
|
Ffrft = fracF_dpw( ...
|
|
Xint,...
|
|
Achirp,...
|
|
H,...
|
|
Cchirp,...
|
|
Aa);
|
|
|
|
%% Mean power spectrum across the DPW
|
|
|
|
Pfft = mean(abs(FFTref).^2,2);
|
|
|
|
Pfrft = mean(abs(Ffrft).^2,2);
|
|
|
|
%% Peak comparison
|
|
|
|
peakFFT = max(Pfft);
|
|
peakFrFT = max(Pfrft);
|
|
|
|
gain_dB = 10*log10(double(peakFrFT/peakFFT));
|
|
|
|
fprintf('\n');
|
|
fprintf('FFT peak power : %.6f\n',double(peakFFT));
|
|
fprintf('FrFT peak power : %.6f\n',double(peakFrFT));
|
|
fprintf('Processing gain : %.3f dB\n',gain_dB);
|
|
|
|
%% Normalize spectra for display
|
|
|
|
Pfft_dB = 10*log10(Pfft/max(Pfft));
|
|
|
|
Pfrft_dB = 10*log10(Pfrft/max(Pfrft));
|
|
|
|
%% Display averaged spectra
|
|
|
|
figure
|
|
|
|
subplot(2,1,1)
|
|
|
|
plot(Pfft_dB)
|
|
|
|
grid on
|
|
|
|
ylim([-60 5])
|
|
|
|
title('FFT Mean Power Spectrum')
|
|
xlabel('FFT Bin')
|
|
ylabel('Power (dB)')
|
|
|
|
subplot(2,1,2)
|
|
|
|
plot(Pfrft_dB)
|
|
|
|
grid on
|
|
|
|
ylim([-60 5])
|
|
|
|
title(sprintf('FrFT Mean Power Spectrum (a = %.6f)', ...
|
|
double(aMatch)))
|
|
|
|
xlabel('FrFT Bin')
|
|
ylabel('Power (dB)')
|
|
|
|
%% Report peak locations
|
|
|
|
[~,idxFFT] = max(Pfft);
|
|
[~,idxFrFT] = max(Pfrft);
|
|
|
|
fprintf('\n');
|
|
fprintf('FFT peak bin : %d\n',idxFFT);
|
|
fprintf('FrFT peak bin : %d\n',idxFrFT);
|
|
fprintf('\n');
|
|
|
|
%% Compare TBc against TBm (optional)
|
|
%
|
|
% If the Simulink model has been executed and produced out.Fsim,
|
|
% compare both implementations.
|
|
|
|
if exist('out','var')
|
|
|
|
fprintf('\n');
|
|
fprintf('TBc vs TBm Comparison\n');
|
|
fprintf('---------------------\n');
|
|
|
|
Ftbm = out.Fsim;
|
|
|
|
%% Dimension check
|
|
|
|
fprintf('TBc size : [%d %d]\n', ...
|
|
size(Ffrft,1), size(Ffrft,2));
|
|
|
|
fprintf('TBm size : [%d %d]\n', ...
|
|
size(Ftbm,1), size(Ftbm,2));
|
|
|
|
assert(isequal(size(Ffrft),size(Ftbm)), ...
|
|
'TBc and TBm dimensions differ.');
|
|
|
|
%% Error metrics
|
|
|
|
err = Ftbm - Ffrft;
|
|
|
|
maxErr = max(abs(err(:)));
|
|
|
|
rmsErr = sqrt(mean(abs(err(:)).^2));
|
|
|
|
refPeak = max(abs(Ffrft(:)));
|
|
|
|
relErr = maxErr / refPeak;
|
|
|
|
%% Results
|
|
|
|
fprintf('\n');
|
|
fprintf('Reference peak : %.9g\n',double(refPeak));
|
|
fprintf('Maximum error : %.9g\n',double(maxErr));
|
|
fprintf('RMS error : %.9g\n',double(rmsErr));
|
|
fprintf('Relative error : %.9g\n',double(relErr));
|
|
|
|
if maxErr == 0
|
|
|
|
fprintf('\nPASS: Outputs are bit-identical.\n');
|
|
|
|
elseif relErr < 1e-5
|
|
|
|
fprintf('\nPASS: Outputs are numerically equivalent.\n');
|
|
|
|
else
|
|
|
|
fprintf('\nWARNING: Outputs differ.\n');
|
|
|
|
end
|
|
|
|
%% Visual comparison
|
|
|
|
frameIdx = 1;
|
|
|
|
figure
|
|
|
|
subplot(3,1,1)
|
|
|
|
plot(abs(Ffrft(:,frameIdx)))
|
|
|
|
grid on
|
|
|
|
title('TBc Output')
|
|
xlabel('Bin')
|
|
ylabel('|F|')
|
|
|
|
subplot(3,1,2)
|
|
|
|
plot(abs(Ftbm(:,frameIdx)))
|
|
|
|
grid on
|
|
|
|
title('TBm Output')
|
|
xlabel('Bin')
|
|
ylabel('|F|')
|
|
|
|
subplot(3,1,3)
|
|
|
|
plot(abs(Ftbm(:,frameIdx) - Ffrft(:,frameIdx)))
|
|
|
|
grid on
|
|
|
|
title('Absolute Error')
|
|
xlabel('Bin')
|
|
ylabel('|Error|')
|
|
|
|
%% Mean power spectrum comparison
|
|
|
|
Ptbc = mean(abs(Ffrft).^2,2);
|
|
|
|
Ptbm = mean(abs(Ftbm).^2,2);
|
|
|
|
Ptbc_dB = 10*log10(Ptbc/max(Ptbc));
|
|
|
|
Ptbm_dB = 10*log10(Ptbm/max(Ptbm));
|
|
|
|
figure
|
|
|
|
plot(Ptbc_dB,'LineWidth',1.5)
|
|
|
|
hold on
|
|
|
|
plot(Ptbm_dB,'--','LineWidth',1.5)
|
|
|
|
grid on
|
|
|
|
ylim([-60 5])
|
|
|
|
xlabel('Bin')
|
|
ylabel('Power (dB)')
|
|
|
|
title('TBc vs TBm Mean Power Spectrum')
|
|
|
|
legend('TBc','TBm')
|
|
|
|
else
|
|
|
|
fprintf('\n');
|
|
fprintf('TBm comparison skipped (out.Fsim not found).\n');
|
|
|
|
end |