Determining the input from the output using IFFT

Hello,

I have an output waveform y(t) after having passed through an S-matrix.I need to get the input signal x(t) and the method followed is:

1.Extract S21 from S-matrix and take IFFT of S21 to get the impulse response h(t)

2.convolve output waveform y(t) with 1/h(t) to get back input signal x(t).

The code is written below:

%Read the S parameter matrix and seperate the freq and s matrix data
sp=sparameters('C:\Users\PrasadNGanes\Desktop\snp_data\2Port.s2p');
freq=sp.Frequencies;
s=sp.Parameters;
%Seperate the S21 column from S-matrix and plot S21
S21=rfparam(sp,2,1);
rfwrite(S21,freq,'C:\Users\PrasadNGanes\Desktop\snp\newS21.s2p');
sp_inv=sparameters('C:\Users\PrasadNGanes\Desktop\snp\newS21.s2p');
figure;
rfplot(sp_inv);
title('S21 plot');

%Take IIFT of S21 to get impulse response h(t) and plot Impulse response
impulse_response=ifft(S21);
figure;
plot(impulse_response);
Fs = 2*max(freq);
Ts=1/Fs;
N=numel(impulse_response);
tvec=(0:(N-1))*Ts;
figure;
plot(tvec,impulse_response);
title('impulse response');

%convolve output y(t) with 1/h(t) to get back the input x(t)
inverse_impulse=1./impulse_response;

%read the output y(t)
data=xlsread('C:\Users\PrasadNGanes\Desktop\snp_data\probe_output.xlsx');
out_time=data(:,1);
output=data(:,2);

convolution_result=conv(output,inverse_impulse);
convolution_time=linspace(min(out_time)+min(tvec),max(out_time)+max(tvec),length(convolution_result));
figure;
plot(convolution_time,convolution_result);
title('input plot after convolution');

%Plot the real input and output
data=xlsread('C:\Users\PrasadNGanes\Desktop\probe.xlsx');
inp_time=data(:,1);
input=data(:,2);
figure;
plot(inp_time,input,'b');
hold on
plot(out_time,output,'r');
title('Actual output in red and input in blue');

However, I see that the input waveform got after running the code and the actual waveform do not match.I have attached the images for reference.

Can somebody help me what needs to be corrected to get the correct result.

[Ganesh Prasad]

I have tried taking the conjugate of the S21 row vector, flipping it and adding to the original S21 and then taking IFFT.But the output is not as expected.

%Read the S parameter matrix and seperate the freq and s matrix data
sp=sparameters('C:\Users\PrasadNGanes\Desktop\snp_data\QVR_MX0100A_InfMode_SldrInFlat_Zin_2Port.s2p');
freq=sp.Frequencies;
s=sp.Parameters;
%Seperate the S21 column from S-matrix and plot S21
S21=rfparam(sp,2,1);
rfwrite(S21,freq,'C:\Users\PrasadNGanes\Desktop\snp\newS21.s2p');
sp_inv=sparameters('C:\Users\PrasadNGanes\Desktop\snp\newS21.s2p');
figure;
rfplot(sp_inv);
title('S21 plot');

fs=26999*10^6;
N=length(freq);
dt=1/fs;
t=(0:N-1)*dt;
S21_t=transpose(S21);
reconstructed_fft = [S21_t, conj(fliplr(S21_t(2:end-1)))];


%Take IIFT of S21 to get impulse response h(t) and plot Impulse response
impulse_response=ifft(reconstructed_fft);
figure;
plot(impulse_response);
title('impulse response');

%convolve output y(t) with 1/h(t) to get back the input x(t)
inverse_impulse=1./impulse_response;

%read the output y(t)
data=xlsread('C:\Users\PrasadNGanes\Desktop\snp_data\probe_output.xlsx');
out_time=data(:,1);
output=data(:,2);

convolution_result=conv(output,inverse_impulse);
convolution_time=linspace(min(out_time)+min(t),max(out_time)+max(t),length(convolution_result));
figure;
plot(convolution_time,convolution_result);
title('input plot after convolution');

[Ganesh Prasad]

I have also tried using rational fit.The code is below:

S = sparameters('2Port.s2p');
freq = S.Frequencies;
data = rfparam(S,2,1);
rat_fit = rational(freq,data);
% Now rat_fit is the rational function object that models the S21 response
% Compute frequency response
[H, f] = freqresp(rat_fit, freq);
% Obtain impulse response by inverse Fourier transform
impulseResponse = ifft(H);
fs=26999*10^6;
N=length(freq);
dt=1/fs;
t=(0:N-1)*dt;

inverse_impulse=1./impulseResponse;
%read the output y(t)
data=xlsread('C:\Users\PrasadNGanes\Desktop\snp_data\probe_output.xlsx');
out_time=data(:,1);
output=data(:,2);

convolution_result=conv(output,inverse_impulse);
convolution_time=linspace(min(out_time)+min(t),max(out_time)+max(t),length(convolution_result));
figure;
plot(convolution_time,convolution_result);
title('input plot after convolution');

The output of the convolution is still not matching the required waveform.

[William Rose]

@Ganesh Prasad,

Perhaps the output is not as expected because

b2 = s21a1 + S22a2

and your code does not account for the S22a2 contribution to the output. Perhaps that term is significant.

[Ganesh Prasad]

Thanks for responding to the query.

I have an output waveform b2 and S21 data from S-matrix but there is no a2 waveform(no input from the load end/matched load) and hence using b2=S21 a1 given a2=0 for analysis.

Is there anyway the input can be determined from the output?

[William Rose]

@Ganesh Prasad, I'm afraid I do not have time to assist on this question at this time.

[Paul]

What is the justification for this?

% 2.convolve output waveform y(t) with 1/h(t) to get back input signal x(t).
inverse_impulse=1./impulse_response;
...
convolution_result=conv(output,inverse_impulse);

That doesn't seem right.

Sometimes it's helpful to verify code for a simple problem with a known solution.

[Ganesh Prasad]

Hello,

This was done considering the below condition:

If x(t)*h(t)=y(t) then x(t)=1/h(t) * y(t) where * is convolution and h(t) is the impulse response.

I have an output waveform y(t) after having passed through an S-matrix.I need to get the input signal x(t) and the method followed is:

1.Extract S21 from S-matrix and take IFFT of S21 to get the impulse response h(t)

2.convolve output waveform y(t) with 1/h(t) to get back input signal x(t).

[William Rose]

@Paul is right, as usual. It is not true that

"If x(t)*h(t)=y(t) then x(t)=1/h(t) * y(t) where * is convolution and h(t) is the impulse response."

When I say DFT below, I mean the discrete Fourier transform, which you can compute with fft().

I would do the deconvolution in the frequency domain. To do this, you will need the frequencies of the S21(f) response to be the same as the frequencies of Y(f)=DFT of y(t). Then you can estimate X(f)=Y(f)./S21(f), and x(t)=invDFT(X(f)). But remember that the DFT assumes that your recording of y(t) represents exactly one period of a periodic signal. Your recording of y(t) has about 3.5 periods. So there will be a lot of high frequency "noise" in the DFT, due to the mismatch between the beginning and the end. This is why you use a window such as Hamming or Hann. If you apply a window to Y(f) before computing its DFT, then the X(f) and x(t) you get from deconvolution will be the input that gives the windowed signal. Another option is to clip y(t) so that it wraps around smoothly, i.e. so that you have exactly 3 cycles of the signal, not 3.5. You may also need to apply a linear adjustment, i.e. a kind of "tilt", so that the beginning and the end of the signal line up smoothly. Then you can apply the opposite tilt to the final x(t) which you get from deconvolution.

[Paul]

@Ganesh Prasad

I think it would be helpful if you shared you s21 data (including the frequencies) or a representative example and the associated y(t). You can use the paperclip icon on the Insert menu to attach a file. And I think we need to better understand exactly how y(t) is generated.

[Ganesh Prasad]

Please see the S21(please change .txt to .s2p extension) and output waveform file(.xlsx) attached.Thanks.

[Paul]

Taking a look at the data

copyfile S21.txt S21.s2p
S = sparameters('S21.s2p')

S = 
  sparameters: S-parameters object

       NumPorts: 1
    Frequencies: [1774×1 double]
     Parameters: [1×1×1774 double]
      Impedance: 50

  rfparam(obj,i,j) returns S-parameter Sij

The frequency data does not start at zero. I think if you're going to use the IFFT approach, you'll have to figure out how to extrapolate S21 to zero.

format short e
S.Frequencies([1 end])

ans = 2×1

1.0e+00 * 1.0000e+06 2.7000e+10

figure
semilogx(S.Frequencies,abs(squeeze(S.Parameters)))

The frequencies are not equally spaced (maybe they're spaced logarithmically?). Before taking using IFFT I think you'll have to do an interpoplation to an equally spaced frequency vector.

figure
plot(diff(S.Frequencies))

After those two issues are sorted out, you'll have to also account for the "upper half" of S21 in the frequency domain. Assuming that the impulse response is supposed to be real-valued, check this answer for a discussion on that topic.

I believe that all of these issues need to be addressed first to get the "full picture" of S21 before proceeding with the frequency domain approach suggested by @William Rose or with a time domain approach as you are attempting.

I didn't look at y(t), I suspect that there will be more work to do on that as well depending on decisions made on how you modify the S21 data.

There may be other options to to address this problem as well.

Answer 1

[William Rose]

@Ganesh Prasad,

Here is a script that does the deconvolution in the frequency domain, using the ideas in my comment above. It also does what @Paul suggested: it interpolates and extrapolates S21(f) to the frequencies of the output signal's Fourier transform.

The script makes five figures, which are mainly to reassure me that everything is working as I expect. The final figure, copied below, shows the output and esitmated input. I am not saying it is perfect.

See comments in the script for details.

[William Rose]

As @Paul mentioned, the frequency range for which S21 is specified, in file S21.txt, is much less than the frequency range of the signal y(t). y(t) is sampled at abut 1280 GHz, so the Nyquist frequency is about 640 GHz. S21 is specified from 1 MHz to 27 GHz, so we must extrapolate to f=0, and from 27 GHz to 640 GHz, in order to do the deconvolution. I am not worried about the extrapolation to f=0, since th extrapolation in this range looks pretty obvious on the plot of S21(f). I am less confident abut the extrapolation from 27 GHz to 640 GHz. I assume in my script that is flat from 27 to 640 GHz, because the plot of looks like it has flattened out by f=27 GHz. The estimated input would probably look different if I made a different assumption for the extrapolation of S21.

[Paul]

William,

Did you by any chance try downsampling y(t) to 54 GHz, which looks like it would be more than sufficient based on the plot of y? That way one could avoid the need to extrapolate S21(f). Just curious how that compares to what's shown above.

Another option may be to simply zero-pad S21 out to the Nyquist frequency, which is effectively interpolating the IDFT of the non-padded S21 in the time domain.

I also recently discovered that when the 'symmetric' flag is used with ifft there isn't a need to explicitly fill out the upper portion of the DFT before taking the IDFT; zeros will work just fine. For example

x = 1:10;
X = fft(x);
X(7:10) = 0;
ifft(X,'symmetric') - ifft(fft(x))

ans = 1×10

0 0 0 0 0 0 0 0 0 0

[William Rose]

@Paul, Good idea about down-sampling y(t). I didn't think of that. Thank you for the tip about 'symmetric' for ifft().

[Paul]

I just realized that zero-padding S21(f) will be a problem using this approach ... and I just realized that downsampling y might lead to other complications. Interesting problem.

[Ganesh Prasad]

@William Rose Thank you for taking the time out and explaining in detail with the code.This will help me verify the results I got from the simulation.

@Paul Thank you for giving your valuable inputs to get the solution to the problem.

Thanks a lot and sorry for the delayed response.

[William Rose]

@Ganesh Prasad, you're welcome.