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

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?)

Answer 1

[David Goodmanson]

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.

[Michael Dinallo]

very helpful, I should have indicated in the question that we do interpolate to common frequencies,

but can we also shift the [-N/2, N/2] interval to an [0, N-1] interval and how might that change your steps in the above example (which I appreciate much - being an infrequent but grateful Matlab user)?

also, how might such a cable S21 compensation steps change if we use a DFT in Matlab? (on an internal beginning with N time points M(t)? would we have to still conjugate and image a symmetric spectral version of the S21 correction/compensation file?

[David Goodmanson]

[ In the answer I fixed a mistake in the line involving S12, so now the code runs]

For an N-point time array, whatever fourier method you use on discrete time points, it's a linear transform from N to N points. That implies N/2 complex amplitudes (or if you prefer, N/2 absolute amplitudes and N/2 phase shifts) from S12 at at N/2 frequencies.

Back to the fft, by shifting the [-N/2, N/2] interval to an [0, N-1] interval I assume you mean within the code, not using N measured frequencies of S12 because that does not work. You can use [0, N-1] in the code in a way, but I think all it does is obfuscate what is going on. For an fft, frequency m and frequency m-N are the same because of aliasing, which means in the N=8 case that

f = 0 1 2 3 4 5 6 7   is the same as   f_alias = 0 1 2 3 ny -3 -2 -1

where the nyquist point can be either 4 or -4. So you could use f, which corresponds to the fft without the use of fftshift. But you still have to get N/2 values from conjugating and flipping S12. Four lines of code would be different:

a = fft(m)                % no fftshift;
S12full = [S12(1:end)   conj(flip(S12(2:end-1)))];
aadj(N/2+1) = 0;          % nyquist point
madj = ifft((aadj));      % no ifftshift

Note that in the code, the line f = (-N/2:N/2-1)*delf is never actually used, so neither would be the line f = (0:N-1)*delf.

[Michael Dinallo]

getting me further along to my calculations - thanks very much - Mike