The use of IFFT

Can any kind soul please help or provide some kickstart?

I have a transfer function/frequency response, h(f). Where h(f) is calculated by applying FFT to the csv data and then V2(f) is divided with V1(f).

Answer 1

[Star Strider]

Your data are frequency and only the corresponding amplitudes of the input and output signal at each frequency. You need to measure the phase data at each frequency as well, and since this is missing, half of your data are missing.

You could calculate the transfer function from the complex fft of both the input and output time domain data, but lacking the phase data, reconstructing the impulse response is not possible.

If you have complex frequency data, you can use the Signal Processing Toolbox invfreqz function to estimate the transfer function, and then use impz to calculate the impulse response.

[Alycia]

How do I measure the phase data at each frequency? Right now all I have is csv data obtained from oscilloscope.

[Star Strider]

You can’t measure the phase from only the data you have. The oscilloscope has to be able to give you the phase information. If the oscilloscope is doing the Fourier transform, there may be a way to get the phase information from it. Otherwise, you’re missing half the information it’s creating in calculating the Fourier transform, and that information is necessary for you to do what you want.

If the original data are stored in the oscilloscope, you might be able to have it give you the complex results. That way, you can calculate the complex transfer function, and from that, estimate the impulse response. It would actually be much better if you can get the original time-domain signals. Then you could import them into MATLAB and do the entire fft and transfer function analysis there.

[Alycia]

The attached is is my original time-domain signals data. By using the code that you have shared months ago (as below), will I be able to calculate the entire fft? or will it be easier to just use the Fourier transform from the oscilloscope?

[d,s,r] = xlsread('Alycia Tay scope_0.csv');    % Load File
t = d(3:end,3);                                 % Time Vector
sig = d(3:end,4:5);                             % Signal Vectors
L = size(t,1);                                  % Record Length
Ts = mean(diff(t));                             % Sampling Interval
Fs = 1/Ts;                                      % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency
Fsig = fft(sig)/L;                              % Normalised FFT
Fv = linspace(0, 1, fix(L/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector
figure(1)
semilogy(Fv, abs(Fsig(Iv,:)))                   % Use ‘semilogy’ For Plot
grid
legend('Signal 1', 'Signal 2')

I'm still quite new to DSP concepts and my only reference is from dspguide.com. Do you have any useful websites or textbooks to share so that I could grasp the concept within weeks?

[Star Strider]

Apologise for the delay. I’m GMT-7 here and it was the end of my day (and energy).

When I open your ‘Trace 3’ file in Excel, it has one time column and one voltage column. What does it represent? I thought you wanted to compare two signals (input and output) and get the transfer function.

[Alycia]

Sorry! This is the correct file. Thanks so much for your help! Very appreciated!

[Star Strider]

As always, my pleasure!

To clarify, ‘1 (VOLT)’ and ‘2 (VOLT)’ are the signals we’re interested in. (I get the impression you’re in Oz or NZ, so Happy Valentine’s Day!)

Not so happy with the results, becasue both signals seem to be the same. I need yet another file to calculate a transfer function.

My code:

[D,S,R] = xlsread('Alycia 1MHz.xls');
Time = D(:,2);
V1   = D(:,3);
V2   = D(:,3);

Ts = mean(diff(Time));                                              % Sampling Interval (Mean Interval Difference)
Td = std(diff(Time));                                               % Test For Consistency (Standard Deviation)
Fs = 1/Ts;                                                          % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency

L = length(Time);                                                   % Signal Length
FT1 = 2*fft(V1)/L;                                                  % Fourier Transform: V1 (One-Sided, Normalised)
FT2 = 2*fft(V2)/L;                                                  % Fourier Transform: V2 (One-Sided, Normalised)
HT  = FT2./FT1;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                 % Frequency Vector
Iv = 1:length(Fv);                                                  % Index Vector

figure(1)
subplot(3,1,1)
plot(Fv, abs(FT1(Iv)))
grid
axis([0  1E+8    ylim])
subplot(3,1,2)
plot(Fv, abs(FT2(Iv)))
grid
axis([0  1E+8    ylim])
subplot(3,1,3)
plot(Fv, abs(HT(Iv)))
grid
axis([0  1E+8    ylim])

Check = [mean(V1-V2); std(V1-V2)];                                  % Check To See If The Signals Are Different

The ‘Check’ assignment returned a vector of zeros, showing the signals are exactly the same.

No worries or rush on my part — this seems to be a relatively straightforward problem, so I’ll stay with it until we’ve got it sorted and get good results.

Meanwhile, I’m eagerly and happily awaiting a new file...

[Alycia]

Happy Valentine's Day to you as well.

awwww... I don't have any new set of data with me right now and the lab will be closed till Tuesday. Instead, I made amendments to the previous data (1MHz.xls). I have shifted the Time and ‘1 (VOLT)’ and ‘2 (VOLT)’ to positive values in this current attached file.

However, the TF seems rather weird and I can't proceed on to IFFT with such disappointing result.

[Star Strider]

When I plotted it, it looks like just broadband noise.

I intend to be here until we get this sorted, so I’ll wait til Tuesday and then see what the new data are.

[Alycia]

Hello! Hope you are doing fine. Please refer to the code as below. Correct me if I am wrong. Thank you!

array = xlsread('AC_Width100nsEdge20ns.xls');
time = array(:,2);
v1 = array(:,3);
v2 = array(:,4);
N = size(time,1);
Ts = mean(diff(time));     
Fs = 1/Ts;           
Fn = Fs/2;           
FT1 = fft(v1)/N;
FT2 = fft(v2)/N;
HF = FT2./FT1;

FreqV = linspace(0,1,fix(N/2)+1)*Fn;
Iv = 1:length(FreqV);

%Magnitude only
mag1 =abs(FT1(Iv));
mag2 =abs(FT2(Iv));
TF = mag2./mag1;

%Magnitude in dB
magdB1 = 20*log10(abs(FT1(Iv)));
magdB2 = 20*log10(abs(FT2(Iv)));

%Phase in degree
phdg1 = (180/pi)*unwrap(angle(FT1));
phdg2 = (180/pi)*unwrap(angle(FT2));

I have a question regarding to your code. Why do you multiple 2 to your FT? In addition, how can I utilize both the mag and phase results to calculate ifft and convolution?

[Star Strider]

Fine here, I hope you are as well.

You’re lucky in having both your input and output signals in the time domain, so you can do all your calculations in the complex domain. You don’t have to bother with magnitude and phase conversions. I did the convolution and ifft to get the impulse response in my code here.

I added a few lines of code to calculate the transfer function and impulse response correctly (it’s necessary to do all these in the complex frequency domain), and added three plots, but otherwise your code is unchanged. The lines I added are commented, and the plots are self-explanatory.

Realising your transfer function in system terms (that is, creating a filter that reproduces it) is probably impossible. I nevertheless gave that a go with the invfreqz call. That isn’t necessary, but I was curious, so I did it. You can delete that calculation in the associated freqz plot call without affecting the rest of the code. I left it in so you can experiment with it if you want to.

The multiplication by 2 in the plot of the fft is necessary to normalise the amplitudes. The fft is actually symmetric about the origin, so half the energy is in each half. Multiplying the single-sided fft by 2 corrects for this. The 'symmetric' argument tells ifft that it is the two-sided fft. This improves the accuracy.

In the ifft call, I multiplied the transfer function created from the original fft data by ‘N’ to ‘denormalise’ it so the amplitudes would be correct.

Here you go:

array = xlsread('Alycia AC_Width100nsEdge20ns.xls');
time = array(:,2);
v1 = array(:,3);
v2 = array(:,4);
N = size(time,1);
Ts = mean(diff(time));     
Fs = 1/Ts;           
Fn = Fs/2;           
FT1 = fft(v1)/N;
FT2 = fft(v2)/N;
HF = FT2./FT1;

FreqV = linspace(0,1,fix(N/2)+1)*Fn;
Iv = 1:length(FreqV);

%Magnitude only
mag1 =abs(FT1(Iv));
mag2 =abs(FT2(Iv));

% Complex Calculations
TFc = FT2./FT1;                                         % Complex Transfer Funciton
TF = abs(TFc(Iv));                                      % Magnitude
Ph = (angle(TFc(Iv)));                                  % Phase (Radian)
ImpRsp = ifft(TFc*N, 'symmetric');                      % Impulse Response

% Call ‘invfreqz’
Wn = FreqV/Fn;                                          % Normalise Frequencies To (0,1)
[b,a] = invfreqz(TFc(Iv), Wn, 10, 11);                  % Estimate Transfer Funciton Coefficients

%Magnitude in dB
magdB1 = 20*log10(abs(FT1(Iv)));
magdB2 = 20*log10(abs(FT2(Iv)));

%Phase in degree
phdg1 = (180/pi)*unwrap(angle(FT1));
phdg2 = (180/pi)*unwrap(angle(FT2));

figure(1)
subplot(2,1,1)
plot(FreqV, 2*TF)
grid
axis([0  1E+8    ylim])
title('Transfer Function')
ylabel('Amplitude')
subplot(2,1,2)
plot(FreqV, Ph*180/pi)
grid
axis([0  1E+8    ylim])
xlabel('Frequency')
ylabel('Phase (°)')

figure(2)
plot(time, ImpRsp)
grid
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')

figure(3)
freqz(b, a, 2048, Fs)

[Alycia]

I was sick since Monday but I'm feeling great now. Thank you so much for the very well explained codes! I am a newbie in matlab, really appreciate that you extended your help for almost a week. Will give it a try later! May I approach you again if I have encountered doubts to the code above?

[Star Strider]

Sorry you were ill. I picked up something last week as well. It was accompanied by the unusual complication of viral labyrinthitis, so that provided a bit of fun.

As always, my pleasure!

Of course you can! My Accepted Answers come with implied support as long as the topic is reasonably close to the original.