Time response from T parameter matrix

Hello,I need to pass a signal(CSV file having time and voltage columns) through a T parameter matrix of 5000 rows of frequency.Can somebody let me know which function to use to process a signal through a T parameter matrix?

[Mathieu NOE]

hello

I don'tunderstand the "T parameter matrix of 5000 rows of frequency" ?

are we speaking about a transfer function ?

[Ganesh Prasad]

Hi, I have converted S parameter matrix of 2*2*5000 to a T parameter multidimensional matrix of 2*2*5000.The T matrix is essentially a 2*2 matrix for 5000 different frequency values.I need to see the response of the T parameter matrix to a time domain data having time and voltage columns.Is there a way to do this?

[Mathieu NOE]

well, which channel of the T matrix would you use ?

IMHO you define which input / output port of T must be used, so we get a 1D array of magnitude and phase from which we make either a FIR or IIR filter using invfreqs or invfreqz

then you can use this filter model and make a time domain simulation using filter

[Ganesh Prasad]

Hello Mathew,

Thanks for the response

I see that

[b, a] = tfdata(H);

H_inf = infreqs(b, a)

where:

• "H" is the transfer function
• "b" and "a" are the numerator and denominator coefficient vectors of "H"
• "H_inf" is the infinite frequency response of "H"

I have just 1-D T parameter vector.How to convert it to transfer function H in the above example?

[Mathieu NOE]

hello again

yes , can you share a T vector ? zip it if it's big

[Ganesh Prasad]

Hello Mathieu,

I have converted using S2T() using the below code.

sparams = sparameters('2Port.s2p');

S = sparams.Parameters;

freq = sparams.Frequencies;

T1=s2t(S);

The S2P file is attached(the extension is .txt .Please change to .s2p while using)

Please let me know if you need any further data.Thanks.

[Mathieu NOE]

here I am just plotting the 4 transfer functions from your data

we recognize that there are in fact only 2 different TF (symetry)

now the next question is up to which frequency must the TF be identified ? this depends of your input signal content

the reason for my question is : do we need an accurate model up to the max plotted frequency or can we somehow simplify the problem by limiting the upper frequency range constraint

%% code if you have the following toolboxes
%   Antenna Toolbox
%   RF Toolbox

% % Load the .s2p file using the sparameters function
% s2pObject = sparameters('2Port.s2p');
%  
% % Extract frequency and S21 parameter
% frequency = s2pObject.Frequencies;
% S21 = s2pObject.Parameters(:,2,1); % S21 parameter
% % Find the magnitude of S21
% S21_magnitude = abs(S21);

%% my home made solution 
% data file is organized s follows : 
% !! freq	mag(S11)	phase(S11)	mag(S12)	phase(S12)	
% !     	mag(S21)	phase(S21)	mag(S22)	phase(S22)	
% !
% #  GHz  S  MA  R 50
data = readmatrix('2Port.s2p','FileType','text'); 
%freq/mag(S11)/phase(S11)/mag(S12)/phase(S12)/mag(S21)/phase(S21)/mag(S22)/phase(S22)
freq = data(:,1)*1e9; % Hz
FRF11 = data(:,2).*(cosd(data(:,3))+1i*sind(data(:,3))); % complex FRF
FRF12 = data(:,4).*(cosd(data(:,5))+1i*sind(data(:,5))); % complex FRF
FRF21 = data(:,6).*(cosd(data(:,7))+1i*sind(data(:,7))); % complex FRF
FRF22 = data(:,8).*(cosd(data(:,9))+1i*sind(data(:,9))); % complex FRF

figure(1)
subplot(2,1,1),semilogx(freq,20*log10(abs(FRF11)),freq,20*log10(abs(FRF22)));
xlabel('Frequency(Hz)');
ylabel('Modulus (dB)');
legend('FRF11','FRF22');
subplot(2,1,2),semilogx(freq,180/pi*(angle(FRF11)),freq,180/pi*(angle(FRF22)));
xlabel('Frequency(Hz)');
ylabel('angle (deg)');
legend('FRF11','FRF22');

figure(2)
subplot(2,1,1),semilogx(freq,20*log10(abs(FRF12)),freq,20*log10(abs(FRF21)));
xlabel('Frequency(Hz)');
ylabel('Modulus (dB)');
legend('FRF12','FRF21');
subplot(2,1,2),semilogx(freq,180/pi*(angle(FRF12)),freq,180/pi*(angle(FRF21)));
xlabel('Frequency(Hz)');
ylabel('angle (deg)');
legend('FRF12','FRF21');

[Ganesh Prasad]

Hello Mathieu,

I agree that the S-matrix is symmetrical.

Please find the code for the FFT of the input waveform.

data=xlsread('C:\Users\PrasadNGanes\Desktop\snp_data\probe.xlsx');

time=data(:,1);

output=data(:,2);

Ts=0.78*10^-12;%difference between samples in probe.xlsx

Fs=1/Ts;

N=length(output);

X=fft(output);

f=(0:N-1)*Fs;

figure;

plot(f,X);

xlabel('Frequency');

ylabel('FFT');

Please see the attached FFT waveform and the input waveform file attached.

[Mathieu NOE]

so your input signal looks like a low pass filtered square wave (frequency = 0.77 GHz + harmonics)

the signal is a bit short so I could only do a single buffer fft (no averaging) to have finest resolution (still limited to 0.2 GHz) , so if you pick two peaks from the time plot you can see the time delta (signal period) is around 1.3e-9 s , therefore my frequency = 0.77 GHz.

the first peak from the fft graph is matching this value (at 0.8 GHz +/- 0.2 GHz resolution)

now if we want to keep the general shape of your signal , that means we need to take care of the upper harmonics as well.

We can see from the fft graph that above 10 GHz the spectrum amplitude drops rapidely , so we can forget about what is above this limit .

so back to our FRF models , wee can say that we don't look for perfect match above 10 GHz

FYI, this is my code to generate the time and fft plots

data=readmatrix('probe.xlsx');
time=data(:,1);
signal=data(:,2);
dt = mean(diff(time));
Fs = 1/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = length(signal);    % for maximal resolution take NFFT = signal length (but then no averaging possible)
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
ind = freq>0;
freq = freq(ind);
spectrum_raw = spectrum_raw(ind);

figure(2),semilogx(freq,20*log10(spectrum_raw),'b');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');


function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[Mathieu NOE]

now let's try to put a low order model on your FRF

here I opted for a very low order (NA = NB = 1) and IMHO it should suffice

now , as we have the B and A coefficients of the IIR filters , we will use them latter on to apply them on your input signal

NB that the invfreqz operation is done with the sampling rate of your input signal , in order afterwards to make the filtering (with filter) : ( Fs = 1.2800e+12 )

% model for FRF11 / FRF22

% model for FRF12 / FRF21

%% code if you have the following toolboxes
%   Antenna Toolbox
%   RF Toolbox

% % Load the .s2p file using the sparameters function
% s2pObject = sparameters('2Port.s2p');
%  
% % Extract frequency and S21 parameter
% frequency = s2pObject.Frequencies;
% S21 = s2pObject.Parameters(:,2,1); % S21 parameter
% % Find the magnitude of S21
% S21_magnitude = abs(S21);

%% my home made solution 
% data file is organized s follows : 
% !! freq	mag(S11)	phase(S11)	mag(S12)	phase(S12)	
% !     	mag(S21)	phase(S21)	mag(S22)	phase(S22)	
% !
% #  GHz  S  MA  R 50
data = readmatrix('2Port.s2p','FileType','text'); 
%freq/mag(S11)/phase(S11)/mag(S12)/phase(S12)/mag(S21)/phase(S21)/mag(S22)/phase(S22)
freq = data(:,1)*1e9; % Hz
FRF11 = data(:,2).*(cosd(data(:,3))+1i*sind(data(:,3))); % complex FRF
FRF12 = data(:,4).*(cosd(data(:,5))+1i*sind(data(:,5))); % complex FRF
FRF21 = data(:,6).*(cosd(data(:,7))+1i*sind(data(:,7))); % complex FRF
FRF22 = data(:,8).*(cosd(data(:,9))+1i*sind(data(:,9))); % complex FRF


%% create model for FRF11 / FRF22
% IIR model
iter = 100;
Fs =   1.2800e+12; % sampling rate of your time signal probe
W = 2*pi*freq/Fs;
ind = freq <20e9; % frequency weighting ;
WT = zeros(size(W));
WT(ind) = 1;
NB = 1;
NA = 1;
[B,A] = invfreqz(FRF11,W,NB,NA,WT,iter);
% my trick to force both values of B numerator to have same value but opposite
% sign, in order to have a true differentiator at the numerator
Bmean = mean(abs(B));
B = Bmean*[1 -1];

% check bode plots
H = freqz(B,A,2*pi*freq/Fs);

figure(1)
subplot(2,1,1),semilogx(freq,20*log10(abs(FRF11)),freq,20*log10(abs(FRF22)),freq,20*log10(abs(H)),'k--');
xlabel('Frequency(Hz)');
ylabel('Modulus (dB)');
legend('FRF11','FRF22','identified model');
subplot(2,1,2),semilogx(freq,180/pi*(angle(FRF11)),freq,180/pi*(angle(FRF22)),freq,180/pi*(angle(H)),'k--');
xlabel('Frequency(Hz)');
ylabel('angle (deg)');
legend('FRF11','FRF22','identified model');


%% create model for FRF12 / FRF21

% IIR model
iter = 100;
Fs =   1.2800e+12; % sampling rate of your time signal probe
W = 2*pi*freq/Fs;
ind = freq <3e9; % frequency weighting ;
WT = zeros(size(W));
WT(ind) = 1;
NB = 1;
NA = 1;
[B,A] = invfreqz(FRF12,W,NB,NA,WT,iter);

% check bode plots
H = freqz(B,A,2*pi*freq/Fs);

figure(2)
subplot(2,1,1),semilogx(freq,20*log10(abs(FRF12)),freq,20*log10(abs(FRF21)),freq,20*log10(abs(H)),'k--');
xlabel('Frequency(Hz)');
ylabel('Modulus (dB)');
legend('FRF12','FRF21','identified model');
subplot(2,1,2),semilogx(freq,180/pi*(angle(FRF12)),freq,180/pi*(angle(FRF21)),freq,180/pi*(angle(H)),'k--');
xlabel('Frequency(Hz)');
ylabel('angle (deg)');
legend('FRF12','FRF21','identified model');