How to turn S11 to time domain by Matlab

Hi,everyone.I run HFSS to get S11 for my circuit,HFSS also provide time domain for S11.I want to obtain the same time domain response by Matlab.However the result is not even close.Can anyone help~Thx!

The input signal is an impulse and I perform the frequency domain simulation from 0~7.5GHz.

And both attached files are the result from HFSS.

Here's my code:

load S11_re.txt -ascii

load S11_im.txt -ascii

freq=1e+9*S11_re(:,1);

S11=S11_re(:,2)+1i*S11_im(:,2);

TDR=ifft(S11);

Fs=2*max(freq);

Ts=1/Fs;

N=numel(TDR);

tvec=(0:(N-1))*Ts;

plot(tvec,abs(TDR))

[Mathieu NOE]

the frequency domain data must contain the phase also (we need to compute the ifft of a complex valued transfer function, here you provide only the modulus)

[Guan Hao]

@Mathieu NOE Sorry,I'm not that familiar with the time domain signal.

Thanks for your help ! The phase data of frequency domain has been uploaded.

I've shared the complex value of S11 and also modified my code.

Answer 1

[Mathieu NOE]

hello again

I am mot sure to understand the shape of the impulse response from your S11_time.txt file :

why only two positive peaks ? from the Bode plot (first figure) I was expecting some kind of damped oscillations - alike what we get from the ifft (nb the symmetrical computation, including negative and positive frequency data)

also I was unsure what is the time unit in the S11_time.txt file ?

% freq domain (transfer function)
load S11_re.txt -ascii % freq + S11_re
load S11_im.txt -ascii % freq + S11_im
freq=1e+9*S11_re(:,1);
frf=S11_re(:,2)+1i*S11_im(:,2);

figure(1),
subplot(2,1,1),plot(freq,abs(frf))
xlabel('freq (Hz)')
ylabel('amplitude')
title('Impulse Response (IR)');
subplot(2,1,2),plot(freq,180/pi*angle(frf))
xlabel('freq (Hz)')
ylabel('phase(°)')


% IR (impulse response) obtained with ifft method
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
    
TDR = real(ifft([frf; frf_sym])); % NB we need the negative and positive frequency complex FRF
TDR = TDR(1:101); % truncation is possible if TDR decays fast enough

Fs=2*max(freq);
Ts=1/Fs;
N=numel(TDR);
tvec=(0:(N-1))*Ts;

% compare to impulse response generated by your software
load S11_time.txt -ascii % time + IR
time = S11_time(:,1);
IR = S11_time(:,2);
time = time/1e9; % assuming time data was given in nanoseconds

figure(2),
plot(tvec,TDR./max(TDR),'b',time,IR./max(IR),'r')
xlabel('time (s)')
ylabel('amplitude')
title('Impulse Response (IR)');
legend('IR from re/im data','IR from file');

[Guan Hao]

@Mathieu NOE Thanks for your help!The S11_time file is the magnitude of the signal,so it returns two positive peak.I'll send another file.

By the way,the time unit is nanosecond.

[Mathieu NOE]

this time it's ok - how can you obtain either two positive peaks or pos / neg ?

results with provided file are matching , I simply removed the header lines so load will work fine (see attached txt file)

all the best

% freq domain (transfer function)
load S11_re.txt -ascii % freq + S11_re
load S11_im.txt -ascii % freq + S11_im
freq=1e+9*S11_re(:,1);
frf=S11_re(:,2)+1i*S11_im(:,2);

figure(1),
subplot(2,1,1),plot(freq,abs(frf))
xlabel('freq (Hz)')
ylabel('amplitude')
title('Impulse Response (IR)');
subplot(2,1,2),plot(freq,180/pi*angle(frf))
xlabel('freq (Hz)')
ylabel('phase(°)')


% IR (impulse response) obtained with ifft method
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
    
TDR = real(ifft([frf; frf_sym])); % NB we need the negative and positive frequency complex FRF
TDR = TDR(1:101); % truncation is possible if TDR decays fast enough

Fs=2*max(freq);
Ts=1/Fs;
N=numel(TDR);
tvec=(0:(N-1))*Ts;

% compare to impulse response generated by your software
load S11_time.txt -ascii % time + IR
% remove the header lines otherwise load will throw an error message
% ============================================================================================
% ANSOFT                                                                              04/18/24
% Terminal S Parameter Plot 2                                                         00:46:05
% 
% --------------------------------------------------------------------------------------------
%         Time [ns]                   St(1,1) []
                                              
       
time = S11_time(:,1);
IR = S11_time(:,2);
time = time/1e9; % assuming time data was given in nanoseconds

figure(2),
plot(tvec,TDR./max(TDR),'b',time,IR./max(IR),'r')
xlabel('time (s)')
ylabel('amplitude')
title('Impulse Response (IR)');
legend('IR from re/im data','IR from file');

[Guan Hao]

@Mathieu NOE The software can display S11 in many ways,I sent the absolute value in former file,so it shows two positive peak.If I display the original data,it shows one postive peak and one negative peak.

Thank you for your help,I do learn a lot about how to do ifft.

By the way,I think the software inputs a rectangular pulse with the width of Ts(Time step) and the magnitude is 1/Ts.So maybe I should create a rectangular pulse and convolution with the impulse response?

But I got a question,how do I know where the rectangular pulse locates on the axis of time?

% create input rectangular pulse
t = 0 : Ts : 2e-9; 
Amp=1/Ts; % amplitude of the pulse
u=[zeros(1,numel(0:Ts:Ts)) Amp*ones(1,numel(Ts:Ts:2*Ts)) zeros(1,numel((2*Ts)+(Ts):Ts:2e-9))]; 
plot(u)
TDR=real(ifft([frf; frf_sym])); % NB we need the negative and positive frequency complex FRF
TDR=TDR(1:numel(u)); % truncation is possible if TDR decays fast enough
TDR1=conv(TDR,u)

[Mathieu NOE]

yes, you need to display the data with positive and negative peaks , otherwise we are comparing oranges and apples.

for the convolution with an impulse , an impulse will not modify your IR as a "true" impulse has falt spectrum (infinite)

as you create a rectangular implse with some samples delay, the effect is to delay your IR with a bit or smoothing (less ringing) as the spectrum of a rectangualr impulse has decaying energy as we go into higher frequencies.

% freq domain (transfer function)
load S11_re.txt -ascii % freq + S11_re
load S11_im.txt -ascii % freq + S11_im
freq=1e+9*S11_re(:,1);
frf=S11_re(:,2)+1i*S11_im(:,2);

figure(1),
subplot(2,1,1),plot(freq,abs(frf))
xlabel('freq (Hz)')
ylabel('amplitude')
title('Impulse Response (IR)');
subplot(2,1,2),plot(freq,180/pi*angle(frf))
xlabel('freq (Hz)')
ylabel('phase(°)')


% IR (impulse response) obtained with ifft method
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
    
TDR = real(ifft([frf; frf_sym])); % NB we need the negative and positive frequency complex FRF
TDR = TDR(1:101); % truncation is possible if TDR decays fast enough

Fs=2*max(freq);
Ts=1/Fs;
N=numel(TDR);
tvec=(0:(N-1))*Ts;

% compare to impulse response generated by your software
load S11_time.txt -ascii % time + IR
% remove the header lines otherwise load will throw an error message
% ============================================================================================
% ANSOFT                                                                              04/18/24
% Terminal S Parameter Plot 2                                                         00:46:05
% 
% --------------------------------------------------------------------------------------------
%         Time [ns]                   St(1,1) []
                                              
       
time = S11_time(:,1);
IR = S11_time(:,2);
time = time/1e9; % assuming time data was given in nanoseconds

figure(2),
plot(tvec,TDR./max(TDR),'b',time,IR./max(IR),'r')
xlabel('time (s)')
ylabel('amplitude')
title('Impulse Response (IR)');
legend('IR from re/im data','IR from file');



% create input rectangular pulse
t = 0 : Ts : 2e-9; 
Amp=1/Ts; % amplitude of the pulse
u=[zeros(1,numel(0:Ts:Ts)) Amp*ones(1,numel(Ts:Ts:2*Ts)) zeros(1,numel((2*Ts)+(Ts):Ts:2e-9))]; 
figure(3),
% plot(u)
TDR=TDR(1:numel(u)); % truncation is possible if TDR decays fast enough
TDR1=conv(TDR,u)
tt = (0:numel(TDR1)-1)*Ts;
plot(tvec(1:numel(u)),TDR./max(TDR),'b',tt,TDR1./max(TDR1),'r')
xlabel('time (s)')
ylabel('amplitude')
title('Impulse Response (IR)');
legend('TDR','TDR1');

[Guan Hao]

@Mathieu NOE Was it nearly impossible to reconstruct the data produced by the software?

I saw some rectangular shape on the peak which was produced by the software.

[Mathieu NOE]

you notice that the IR lacks a bit of time resolution

we would need to acquire data at a faster rate (and your frequency data must also contain more high frequency data)

the rectangular shape is simply due to a "relatively" low sampling rate - so as you would obtain after downsampling a sharp peak , result is a coarse peak with rectangular shaped peak

[Guan Hao]

@Mathieu NOE I see.Thank you for your help sir!

[Mathieu NOE]

my pleasure !

[Guan Hao]

@Mathieu NOE Hi,sir.Sorry to bother you again.I'm thinking that maybe there's a possibility to obtain the third little peak (Multi-reflection) by the former two large peak (Two huge reflection in my circuit) by doing correlation.I've tried for few days but I can't get my result correct.The data of S11 was attached as an txt file.

Thx a lot !!!

%%% IFFT
% Reconstruct HFSS signal
% freq domain (transfer function)
load onepeakhighcircuit_re_im.txt -ascii % freq + S11_re
freq=1e+9*onepeakhighcircuit_re_im(:,1);
frf=onepeakhighcircuit_re_im(:,2)+1i*onepeakhighcircuit_re_im(:,3);
% IR (impulse response) obtained with ifft method
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
 
Fs=2*max(freq);
Ts=1/Fs;
TDR=real(ifft([frf; frf_sym])); % NB we need the negative and positive frequency complex FRF
TDR=TDR(1:101); % truncation is possible if TDR decays fast enough
N=numel(TDR);
tvec=(0:(N-1))*Ts;
plot(tvec,TDR/max(TDR),'LineWidth',2)
hold on
grid on
% Cross Correlation
v=3e8;
t1=0.1/v;  % First reflection ending time,0.1 is the physical length of my circuit
l=t1/Ts;   % Number of points
R1=TDR(1:1:l+1); % First peak TDR
R2=TDR(l+1:1:(2*l)+1); % Second peak TDR
[c,lags]=xcorr(R1,R2); % Obtain multi-reflection from correlation
c1=[TDR(1:1:15);c;TDR(29:40)]; % Except the correlation part,the rest should remain the same
Tvec=(0:numel(c1)-1)*Ts;
plot(Tvec,c1/max(TDR),'--','LineWidth',3)
xlabel('time (ns)')
ylabel('normalized amplitude')
set(gca,'XLim',[0 2e-9])
set(gca,'XTick',0:2e-9/4:2e-9)
set(gca,'XTickLabel',0:0.5:2)
legend('IFFT From S11','Correlation Model')

[Mathieu NOE]

hello

no problem

I had a look at your code and played a bit with it , but I am not sure what method / principle you are using here to "create" that 3rd peak

so I may not really help you unless there is a paper or publication that coule help me understand what you are trying to achieve

[Guan Hao]

@Mathieu NOE Thanks for your reply.

When a signal is traveling in the circuit , it will generate the reflection when the impedance varies. There're two greatest impedance deficit in the circuit that produce two peak of reflection.When I send a signal from input,first reflection occurs at the first impedance deficit and the second reflection occurs at the second impedance deficit,both of the reflections travel back to input port.However,the reflection produced by the second peak produced another reflection when it met the first impedance deficit and it bounced between these two impedance deficit.That's why the third peak was seen in the IR. (multi-reflection)

And the third peak was related to the first and second peak.That's the reason why I'm trying to do correlation to try to predict the third peak.

Hope this is helpful !!!

[Mathieu NOE]

tx for the explanation

before coding it, simply , how do you predict when the third peak should be time localised from the two first peaks ? what is the math behind this ?

[Guan Hao]

@Mathieu NOE I can't predict the amplitude and the time precisely by simple math,that's why I'm doing IFFT from frequency domain.

But I can roughly predict when the third peak should be time localised.

t3 will be where the third reflection was located on time axis.

[Mathieu NOE]

does t1 and t2 correspond to the two first positive peaks

then if t3 = t1 +2*t2 , the third peak should appear around t = 1.45e-9 s

what if the reflection is so attenuated that you don't see the effect anymore ?

[Guan Hao]

@Mathieu NOE

The second reflection peak should be the negative peak which is located at 6e-10s.

And t1=2d1/v, t2=2(d1+d2)/v, t3=(2(d1+d2)+2d2)/v=t2+(t2-t1)=2t2-t1=1.15e-9s.

The third peak should be at 1.15ns as predicted.However,the third peak occurs earlier than 1.15ns.

TBH,I'm designing a taper,which means the impedance is function of x(length).So I think that maybe the reflection path of third peak was wrong (since it's just my hypothesis).

I've checked that the third peak will cause a significant error in frequency domain for my model,so it's non-neglegible.

I'm curious about how did you do the correlation,it seem that you were more closer to the answer than me!!!

Thx for your reply again.

[Mathieu NOE]

hello again

yes you're correct , I made a mistake in the t3 formula

I'm still stratching my head with your problem - not sure to find out where is the issue and how to solve it

[Guan Hao]

@Mathieu NOE Hello again ! I got some trouble doing fft this time...

If there's two pulse 1 and -1 at the beginning and the end ,the rest of the points are all zero (as my code) .I was expecting to see a periodic signal in frequency domain(As the attached figure shows).However, the fft result isn't periodic.Did I miss something? Thx a lot!

Fs = 50e9;          % Sampling frequency
T = 1/Fs;           % Sampling period
Fn = Fs/2;          % Nyquist frequency
t = 0:T:0.667e-9;   % Time vector
L = length(t);      % Signal length
m = zeros(1,L-2);
n = [1 m -1];        % TDR starts and end with impulse of value 1 

stem(t*1e9,n,'LineWidth',1);
xlabel('time (ns)')
title('Time Domain Signal')

m = zeros(1,L-1);
n = [1 m -1];        % TDR starts and end with impulse of value 1 
X = fft(n);
f = Fn*(0:(L))/L;

subplot(2,1,2)
plot(f*1e-9,abs(X./max(X)),'LineWidth',1.5)
xlabel('GHz')
title('Frequency Domain Signal')

[Mathieu NOE]

hello

I changed a bit your code so that I coud add some zeroes before and after the pulses

nz = 2; % choose number of zeroes before and after pulses

see the effect of changing nz on the result

to get the frequency domain signal you expected , you needed to have a least one zero (before and after the pulses)

t = linspace(0,0.667e-9,200);   % Time vector (please choose even number of points)
T = mean(diff(t));           % Sampling period
Fs = 1/T;
L = length(t);      % Signal length
nz = 2; % choose number of zeroes before and after pulses
m = zeros(1,L-2-2*nz);
n = [zeros(1,nz) 1 m -1 zeros(1,nz)];        % TDR starts and end with impulse of value 1 

subplot(2,1,1)
stem(t*1e9,n,'LineWidth',1);
xlabel('time (ns)')
title('Time Domain Signal')

% single sided FFT 
Y = fft(n);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/L*(0:(L/2));

subplot(2,1,2)
plot(f*1e-9,P1,'LineWidth',1.5)
xlabel('GHz')
title('Frequency Domain Signal')

[Guan Hao]

@Mathieu NOE Thanks for your help! I'm thinking that the zero points in frequency domain are'nt accurate,it should be 1.5GHz,3GHz...etc.

By the way,why did the periodic show by adding extra zeros in front and behind the pulses?

[Mathieu NOE]

in order to see the frequency of repetition rate between pulses , you need to have a time signal wich is a train of pulses , so more than 1 or pulses

there is a good explanantion here about what is the fft output of a single rectangular pulses (with different duty ratio) and what it becomes when you do the same fft then on a train of pulses

python - FFT of a periodic pulse signal - Signal Processing Stack Exchange

your code again modified to treat a train of (positive only) pulses distant of 0.667e-9 s , so yes the fft fundamental frequency should be 1/0.667e-9 = 1.4993 GHz

if you do the same with alterning polarity pulses , your frequency is devided by 2

ti = 0.667e-9; % time delta between pulses 
Np = 25; % number of pulses 

t = linspace(0,ti*(Np+1),10*(Np+1));   % Time vector (please choose even number of points)
T = mean(diff(t));           % Sampling period
Fs = 1/T;
L = length(t);      % Signal length
nz = 2; % choose number of zeroes before and after pulses
indp = nz + (0:Np-1)*floor(L/Np); % indices of pulses
signal = zeros(1,L);        % init signal with zeroes
signal(indp) = 1;        % pulses have all same positive polarity 
% signal(indp) = -(-1).^(1:numel(indp));        % pulses have alterning pos / neg polarity

subplot(2,1,1)
stem(t*1e9,signal,'LineWidth',1);
xlabel('time (ns)')
title('Time Domain Signal')

% single sided FFT 
Y = fft(signal);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/L*(0:(L/2));

subplot(2,1,2)
plot(f*1e-9,P1,'LineWidth',1.5)
xlabel('GHz')
title('Frequency Domain Signal')

[Guan Hao]

@Mathieu NOE Thank you again.I think I understood to create train of pulses in order to see the frequency of repetition.But the result seems a little bizzare to me,there are several pulse-like signal in frequency domain.I'm expecting to see something like the attached figure which is chubby and repeats every 1.5GHz

[Mathieu NOE]

hello again

are you sure that the figure you post is the spectrum of a pulse train with time distance = 0.667e-9 s ?

if we say that a pulse is a very narrow rectangular pulse , and we consider a train of such narrow rect pulses then we should get the result as written in this article

python - FFT of a periodic pulse signal - Signal Processing Stack Exchange

[Guan Hao]

@Mathieu NOE I'm not sure that the spectrum is exactly the pulse train...It's just a guess.

What I want to do is try to construct the frequency response with two pulse.

Anyway, I've learned how to solve the periodic problem from you! Thanks a lot again!

Answer 2

[Guan Hao]

@Mathieu NOE Hi,sir! I got a really simple question here. I attempted to describe a phenomenon in frequency domain by building a math model. If the model worked, it should obtain the same frequency response as the simulation softwares. However, it's too complicated to achieve my goal on frequency domain.

Let's say there were two factors A & B that caused my model wrong in frequency domain. Was it possible to turn the frequency response F1 which contained A & B and F2 which only contained A into time domain, turning the result of F1-F2 back to frequency domain,so that I could know how much the factor B influence in frequency domain.Thanks!

[Mathieu NOE]

hello

I am not 100% to understand how you want to proceed

why do you have this problem with the first frequency domain model ? can you share some data and code ?

[Walter Roberson]

You should start your own Question for this.