Calculate an Input Signal using tfestimate, FFT and IFFT

Overview:

tfestimate is used to calculate the Transfer Function (TF) of a system by inputting a known signal (swept sine wave) and measuring the output signal. This Transfer Function must then be used to correct other output (measured) signals in order to obtain the original input signals.

1. INPUT1 -> SYSTEM [TF ? ] -> OUTPUT1 ------> Calculate TF using tfestimate
2. INPUT2 ? -> SYSTEM [TF] -> OUTPUT2 -------> Calculate Input 2,...,n

I have broken the problem into 4 main steps and have provided sample code and asked questions accordingly:

(Assistance or advice on any of the steps or an alternative solution would be greatly appreciate)

1. Using tfestimate to calculate the Transfer Function of a system
2. Calculating the FFT of the output signal
3. TF correction in the frequency domain
4. Performing IFFT on corrected signals

1. Using tfestimate to calculate the Transfer Function of a system

a. Input1 Signal: Swept sine wave

(0.2Hz – 2kHz in practice but for demonstrative purposes a signal with components at 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500 Hz was used to represent the input)

b. Output1 Signal:

• Sampled at 4kHz (fs=4000)
• Approximately 60000 samples were acquired (N~60000)

Questions:

1. Is it possible to calculate the Transfer Function of a system using tfestimate and the above input and output signals?
2. With [Txy,F] = TFESTIMATE(X,Y,WINDOW,NOVERLAP,NFFT,Fs); is it correct to say that the values of Txy are calculated at the discrete frequencies of F. (ie. Txy(1) is the TF of the system at F(1))
3. Are the parameters used to calculate tfestimate acceptable?

I have attached the matlab script I wrote and have included some sample code that should explain the procedure that I have been following.

The input signal comprises equal time steps of sine waves at varying frequencies, each with equal magnitude and no phase delay. Magnitude and phase modifications are made to the input signals and at each frequency to produce the output signal. (see Construct Signals in the script)

If my assumptions are correct then the sample code below should correctly calculate the TF of the system using tfestimate (see tfestimate in the script)

NFFT = 2^nextpow2(length(Vinsweep)); %Next power of 2 after number of samples in Input

[Txy,F] = tfestimate(Vinsweep, Voutsweep, [], [], NFFT, fs); %Calculate tfestimate

The results are plotted and compared to magnitude and phase bode plots calculated from the input and output signals (see Gain and Phase calculations for Bode plots and Compare plots in the script)

While similar at the lower frequencies there seems to be a lot of error in the TF from tfestimate at higher frequencies, this has however been noted before and seems to be common for tfestimate.

Linear interpolation at fixed discrete frequencies was used to get a smooth transfer function that was used for further calculations and steps in the process (see Linear Interpolation of the TF in the script)

fixedTxy = interp1(F,Txy,fixedF); %Interpolate TF at fixed freq

.

2. Calculating the FFT of the output signal

Firstly hypothetical output signals (output2, … , outputn) were produced for simulation:

a. Sampled at 1000Hz (fs=1000)

b. Approximately 60000 samples (N≈60000)

In practice these signals are EMG signals and would have a frequency range of approximately 0.5Hz – 500Hz.

The Matlab FFT function was used to calculate the Fourier Transform of the time domain signals to the frequency domain. This allows the magnitude and phase to be calculated at discrete frequencies.

Questions:

1. Is it possible to determine at which frequencies the Fourier Transform has been calculated? For instance can the frequencies be calculated as follows: |fax_Hz = [0:fs/N:(fs/N)*(N-1)]|Where N is the number of samples in the signal and fs is the sample frequency.
2. What is the implication of an Odd or Even N value and how should it be accounted for?
3. Should an alternative method be used to convert the time domain signals into the frequency domain in order to be able to correct for the transfer function?

For demonstrative purposes the following signals were produced: (see Hypothetical Output Signals in the script)

output2 = similar to input1 (Signal with components at 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500 Hz)

output3 = similar to output1 (Output2 with modifications made to both magnitude and phase at each frequency)

output4 = Swept sinewave with fixed amplitude and phase at frequencies 5; 10; 15; 20 Hz

output5 = Swept sinewave with fixed amplitude and phase at frequencies 0.1; 1; 10; 100 Hz

output6 = output4 with varying amplitude (0.5, 1, 2, 3) and varying phase (0, 30, 45, 60 degrees)

output7 = output5 with varying amplitude (0.5, 1, 2, 3) and varying phase (0, 30, 45, 60 degrees)

output8 = fixed frequency sine wave (f=10Hz) with no phase delay

output9 = fixed frequency sine wave (f=100Hz) with 30 degree phase delay

.

The FFT of output 2 – 9 were then calculated using the matlab FFT function (see Calculate FFT of Hypothetical Outout signals in the script)

fftSig=fft(outputSig);

The magnitude (in Decibels) in the frequency domain can be calculated by,

20*log10(abs(fftSig)))

and phase (in Degrees) in the frequency domain, can be calculated by

(180/pi).*angle(fftSig))

.

3. TF correction in the frequency domain

This was completed in 2 main steps:

1. Interpolation of the TF at the frequencies at which the FFT was calculated
2. Implementing magnitude and phase correction in the frequency domain

Questions:

1. Is it acceptable to interpolate the TF to the frequencies at which the FFT was calculated in order to do correction? And are the frequencies used correct?
2. Can this correction be done while the FFT signal and the TF are both in polar coordinate form?

correctedSig = fftSig2./interpTF

The discrete frequencies used for smoothing and the previously smoothed TF were used for interpolation to ensure that the TF remained smooth across all frequencies. (see Interpolate TF at frequencies of FFT in the script)

TF correction was implemented with both the TF and the interpolated FFT signals were in polar coordinate form. This could also be done by calculating the magnitude and phase at each frequency and then dividing the magnitude and subtracting the phase. (Code snippet above)

.

4. Performing IFFT on corrected signals

This is where some errors start to surface but I have been unable to find a way to do this correctly or an alternative solution. All the research I have done on how to perform this correctly has raised more questions than it has answered. I have laid out the procedure I followed but I have not yet been able to get it working correctly. While the reconstructed time signals, well the real part at least, are similar to the original signals with what seems to be the correct magnitude and phase correction there is some high frequency noise and discontinuities introduced. The time domain signals also possess both real and imaginary components which is incorrect.

2 main procedures were followed to complete this phase:

• The signals were reconstructed in the frequency domain by taking the first N/2 samples and then calculating the conjugate pairs, flipping those and recombining the samples with the DC sample and the sample at the Nyquist frequency.

%==========================================================================
%%Reconstructing Signals
%==========================================================================
function [reconSig] = reconstruction(correctedSig, N)

    if mod(N,2) == 0 % Even

        halfSig = correctedSig(2:N/2);  %Get Half
        conjHalf = conj(halfSig);       %Complex Conjugate of half
        flipHalf = conjHalf(end:-1:1);  %Flip left to right

        %Reconstruct
        reconSig = [correctedSig(1);halfSig(:);correctedSig((N/2)+1);flipHalf(:)];

    else %Odd
        halfSig = correctedSig(2:(N+1)/2);  %Get Half
        conjHalf = conj(halfSig);           %Complex Conjugate of half
        flipHalf = conjHalf(end:-1:1);      %Flip left to right

        %Reconstruct
        reconSig = [correctedSig(1);halfSig(:);flipHalf(:)];

    end
end

• Performing the IFFT on the reconstructed samples

Questions:

1. How to account for both ODD and EVEN number of samples (N) in the signal when calculating the complex conjugates?
2. What do you do with the samples at DC and at Nyquist frequency?
3. Can you perform the IFFT, as is, on the reconstructed signals?
4. How to construct frequency domain signals that do not have imaginary components once they are converted back to the time domain?

Firstly the signals were reconstructed using the function reconstruction which is shown above (see Reconstructing Signals in the script)

Then the Matlab IFFT function was used to convert the reconstructed signals back to the time domain (see Perform IFFT of the Reconstructed Signal in the script)

These signals were then compared to the original output signals (see Plot the Reconstructed Output vs Original Output Signals on sub plots in the script)

I understand that this is a complex problem but I have tried to lay it out as simply as possible and any assistance or advice with any aspect of the problem would be greatly appreciated.

Answer 1

[Jeremy]

I recently went through a similar exercise and so I can provide some answers..

1. Is it possible to calculate the Transfer Function of a system using tfestimate and the above input and output signals?

Yes but it will somewhat meaningless with those signals since they do not have content at any of the same frequencies. If you have a damped harmonic oscillator function it will give you a very meaningful results with a swept sine (use chirp function) input.

2. With [Txy,F] = TFESTIMATE(X,Y,WINDOW,NOVERLAP,NFFT,Fs); is it correct to say that the values of Txy are calculated at the discrete frequencies of F. (ie. Txy(1) is the TF of the system at F(1))

yes, assuming NFFT = Fs.

3. Are the parameters used to calculate tfestimate acceptable?

depends on how discrete your input signals are and how much frequency resolution you need, I usually use 1Hz bins as my starting point.

1. Is it possible to determine at which frequencies the Fourier Transform has been calculated? For instance can the frequencies be calculated as follows: fax_Hz = [0:fs/N:(fs/N)*(N-1)] Where N is the number of samples in the signal and fs is the sample frequency.

Yes, with fft there is no need to use a length that is a power of two, I would just use the entire signal, then your delta-f is just the sampling rate divided by the #samples. fax_Hz = [0:fs/N:fs/2]

2. What is the implication of an Odd or Even N value and how should it be accounted for? I forget exactly but I think the result has a duplicate value for either 0Hz or Nyquist with an odd NFFT. I would just avoid using an odd NFFT, there is really no need to do it.

3. Should an alternative method be used to convert the time domain signals into the frequency domain in order to be able to correct for the transfer function?

In general, your approach should work.

1. Is it acceptable to interpolate the TF to the frequencies at which the FFT was calculated in order to do correction? And are the frequencies used correct?

I have done this to calculate a coherent power time history. like the TF, you need an averaged spectra to get coherence, but you need the complete FFT to get back to the time history. So you will need to interpolate down to your FFT frequencies (interp1). This is where your bin width selection comes in, a large bin width will give you a smoother TF, but will spread out any narrow peaks.

I would have to think about it more, but make sure you do not need to be multiplying by the square root of the TF. Coherence is based on the power-level result and so I needed to multiply by the square root.