Signal compensation of Time signal via FFT and cable loss frequency spectrum

1 view (last 30 days)
we apply an FFT to a measured time signal M(t) of N sampled points, the sampled points have a constant sampling rate of delta-t. M(t) becomes m(f) having magnitude and phase consisting of N points, the Nyquist frequency is f=N/(2*delta-t), above this frequency m(f) symmetrically repeats itself. The FFT of M(t) is in preparation for correction/compenstion of the measurement cable frequency spectrum magnitude losses (roll-off) and transit time compensation (phase-delay of cable spectrum);
we acquire the cable magnitude and phase response frequency spectrum via a vector network analyzer S-parameter, specifically S21, consistency of a number of points that meets or exceeds the Nyquist frequency of m(f) or f=N/(2*delta-t).
Next, we divide the FFT of the measured time signal by the cable mag/phase spectrum,
Problem - we have S21 parameter spectrum only up to the Nyquist frequency of the FFT measured time signal -
Question: do we need to also create an above Nyquist frequency spectrum of the cable S21 to also divide the symmetric above Nyquist frequencies of the measured time signal for a causal IFFT of the compensated/cable-loss corrected time signal?
If so, how?
(frequency decimation of the cable S21 beginning at Nyquist and creating a symmetric cable S21 for division into all N-points of the FFT?)

Answers (1)

David Goodmanson
David Goodmanson on 29 Jun 2022
Edited: David Goodmanson on 29 Jun 2022
Hi Michael,
Although the Nyquist frequency is N/(2*delta_t), you have to consider both positive and negative frequencies. Assuming N has an an even number of points, the appropriate frequency array for M(f) is (-N/2:N/2-1)*delta_f. The zero frequency point is at array position N/2+1 with positive frequencies above and negative frequencies below that point.
Here delta_f, the array spacing, equals fs/N where fs is the sampling frequency for the time signal. fs = 1/delta_t by definition, so the two deltas must necessarily satisfy delta_f*delta_t = 1/N for an fft.
The fftshift function swaps the two halves an arrray, meaing that fftshift(M(f)) matches the frequency array above.
For S21 you have to make sure that delta_f for that measurement matches that of the frequency array for the fft, so interpolation might be required. It may be necessary to extrapolate S21 down to zero frequency, but hopefully it is approximately constant down there. You have the positive S21 frequencies, and to find S21 for the negative frequencies, you can use S21(-f) = conj(S21(f)). So to match the lower half of the frequency array, S21 is flipped and complex conjugated.
An example is below. The function S12 is no attempt at all to create a realistic cable transfer function but merely to demonstrate an effect. Zooming in on the plot shows a slight time delay in the rise of the adjusted pulse compared to the original.
For even N there is a slight complication. For example N = eight points, assuming delta_f = 1 for simplicity,
the shifted frequency array looks like
-4 -3 -2 -1 0 1 2 3.
There is the dc point, and all the positive and negative frequencies pair up except for the nyquist point -4. The fft of a real function in the time domain makes the corresponding nyquist value real, but dividing by S12 will make the resulting value complex, which will give a complex result going back to the time domain. Obviosly that is not wanted, so the code below sets that value to zero. In most cases it is small anyway. Odd N does not have that problem. For N=7, for example, a properly set up freq array is the same as above only with -4 missing.
N = 1000;
delt = 1e-8;
t = (0:N-1)*delt;
m = (t-1e-6).*exp(-1e7*(t-1e-6));
m(m<0) = 0;
m = m/max(m);
a = fftshift(fft(m)); % m(f)
delf = 1/(delt*N);
f = (-N/2:N/2-1)*delf;
fadj = (0:N/2)*delf; % analyzer frequencies
S12 = (.99+3e-8*2*pi*i*fadj).^3;
S12full = [conj(flip(S12(2:end))) S12(1:end-1)];
aadj = a./S12full;
aadj(1) = 0; % nyquist point
% use ifftshift to put value corresponding to f = 0 back at the
% beginning of the array
madj = ifft(ifftshift(aadj));
figure(1); grid on
plot(t,m,t,madj)
A bit more about the nyquist point. For even N, it corresponds to both frequency N/2 and -N/2. The code starts with S12 values corresponding to positive frequencies 0:N/2 and makes adjustments to match the shifted frequency array. Nyquist ends up corresponding to -N/2.
  3 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!