Transfer function extraction from frequency response

Hello,

I am trying to get a step response to a system whose magnitude and impulse reponse data for different frequencies are known.

num= xlsread('C:\Users\PrasadNGanes\Desktop\Try4_TF.xlsx');

FrequencyVector= num(:,1);

mag = num(:,2);

phase = num(:,3);

phase = deg2rad(phase);

ResponseData = horzcat(mag, phase);

sys = frd(ResponseData,FrequencyVector);

[num,den] = bode(sys);

H = tf(num,den);

step(H);

I am getting an error stating :

Error using DynamicSystem/bode

The "bode" command operates on a single model when used with output arguments.

How to resolve this error.

Answer 1

[Star Strider]

You will need to use the System Identification Toolbox for this.

Then:

ResponseData = mag.*exp(1j*phase);

data = idfrd(ResponseData,FrequencyVector,SamplingInterval);

H = tfest(data)

And proceed from there.

See the documentation on idfrd for details.

.

[Star Strider]

The frequencies all have to be below the Nyquist frequency (on-half the sampling frequency). That is an essential requirement of sampled systems. Be certain that you supplied the correct sampling interval to idfrd.

[Ganesh Prasad]

Thanks for the help with the code.

I am getting an error "Cannot simulate the time response of FRD models."

and a warning stating-"Frequency points above the Nyquist frequency have been removed from the data. " The snapshot is attached.

How do i remove this ?

The code is below:

num= xlsread('C:\Users\PrasadNGanes\Desktop\Try4_TF.xlsx');

FrequencyVector= num(:,1);

mag = num(:,2);

phase = num(:,3);

phase = deg2rad(phase);

ResponseData = mag.*exp(1i* phase);

sys = idfrd(ResponseData,FrequencyVector,0.1e-6);

H = tfest(sys,20,20);

t=0:0.1:10;

u=sin(t);

y=lsim(sys,u,t);

[Ganesh Prasad]

Ok.Why isn't lsim() not generating response for Frd objects?

[Star Strider]

I replied to the error message in my previous Comment:

• The frequencies all have to be below the Nyquist frequency (one-half the sampling frequency). That is an essential requirement of sampled systems. Be certain that you supplied the correct sampling interval to idfrd.

I suggest using the compare function to see how well the estimated system fits the supplied data. Start with a simple system (for example 2 poles and let tfest choose the zeros), then increase it until you get a good fit. Tweak (reduce) the number of zeros after that to see if the fit improves. If it does not, use the original estimate.

Frequency response models are not generally able to produce time domain results. You might be able to reconstruct your identified system as a separate, new system and create time domain results with it. You can get the various components (numerator, denominator, sampling interval) from the tfest output. Also, I suggest using step first, then lsim with your own inputs. Be sure the sampling frequency of the time vector and ‘u’ match the sampling frequency of the system. The way I construct time vectors that meet that criterion is:

L = length of the time vector in units of the time vector;

Fs = sampling frequency;

t = linspace(0, L*Fs, L*Fs+1)/Fs;

Then use ‘t’ to create ‘u(t)’.

EDIT — (20 Feb 2024 at 15:56)

Corrected typographical errors, expanded discussion and explanation.

[Ganesh Prasad]

The sampling frequency of the freq response data I have is 100kHz.The sampling interval time is given as 1e-5.I am getting a tf of 100 poles and zeros of very long complex numbers.What can I do to reduce the number of poles and zeros and make them as short as possible?

[Star Strider]

Start with a simpler system (for example 2 or 3 poles and let tfest choose the zeros), then gradually increase the number of poles until you get a good fit. Then experiment by decreasing the number of zeros to see if the fit improves. If it does not, use the original or the best fit values. (This is interactive using the compare function and assessing the resullts, so it may take a few minutes.)

[Ganesh Prasad]

Thanks for the advice.

I have implemented a 2 pole and zero transfer function and seeing the step response to the system.I have attached the bode plot and the step response plot.The code is below:

num= xlsread('C:\Users\PrasadNGanes\Desktop\Try4_TF.xlsx');

FrequencyVector= num(:,1);

mag = num(:,2);

phase = num(:,3);

phase = deg2rad(phase);

ResponseData = mag.*exp(1i* phase);

sys = idfrd(ResponseData,FrequencyVector,0.1e-6);

H = tfest(sys,2);

[num]=[-1.26e+07+1.83e+04i -1.84e+13 - 4.98e+11i];%taken from H.Numerator

[den]=[1+ 0i -2.04e+06+8.94e+04i -2.9e+13-1.08e+11i];%taken from H.Denominator

Ts=0.1e-6;

sys1=tf(num,den,Ts);

bode(sys1);

t=0:0.1e-6:10;

[y,t]=step(sys1,t);

plot(t,y);

What do i need to change in the code to get a valid step response?

[Aquatris]

your bode diagram indicates that you have mostly a constant gain as your system. Neither your magnitude nor your phase change that much. Is that actual data from a system?

You can try to use the 'EnforceStability' option when calling tfest to for a stable system estimate, so you can do step. However accuracy will probably be worse.

opt = tfestOptions('EnforceStability',true);
H = tfest(sys,2,opt);

I think you do not have full grasp of the basics. So I recommend you first revisit what is a transfer function, what is a bode diagram, how the frequency domain relates to time domain and then try to do this.

[Ganesh Prasad]

Yes,I see that the magnitude and the phase of the system does not change much over a wide freq range.

However, my question is why the output is so different from the input.For a sine wave input , the output should not be undergoing a drastic change in magnitude and phase right?

For a sine wave input,I am getting a output which is drastically different.I need an explanation as to how to correct this.

(I am new to this field of work but I need to analyse this issue to go ahead with my simulation work).

Solution to the issue would be greatly appreciated.

[Star Strider]

This is not perfect, however it is likely the best that can be done.

There are severe problems with this system, not the least of which being that the frequency vector is in steps of Hz. The Nyquist frequency is the maximum frequency, so the sampling interval isthe inverst of two times that frequency (that being the sampling frequency). With those changes, I could ‘sort of’ get this to work (and it turns out that a 20-order transfer function actually works best, according to the compare plot).

Unfortunately, the best fit produces 3 right-half-plane poles (one real and two complex), so the system is unstable.

format shortE
num = readmatrix('Try4_TF.xlsx')

num = 218×3

1.0e+00 * 0 9.9973e-01 4.4975e+01 1.0000e+05 9.9973e-01 4.4975e+01 2.0000e+05 9.9973e-01 4.4975e+01 3.0000e+05 9.9973e-01 4.4975e+01 4.0000e+05 9.9973e-01 4.4975e+01 5.0000e+05 9.9973e-01 4.4975e+01 6.0000e+05 9.9973e-01 4.4975e+01 7.0000e+05 9.9973e-01 4.4975e+01 8.0000e+05 9.9973e-01 4.4975e+01 9.0000e+05 9.9973e-01 4.4975e+01

FrequencyVector= num(:,1); 
mag = num(:,2); 
phase = num(:,3);
Fs = 2*max(FrequencyVector);
Ts = 1/Fs

Ts = 
   3.5714e-11

figure
tiledlayout(2,1)
nexttile
plot(FrequencyVector, mag2db(mag/max(mag)))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
nexttile
plot(FrequencyVector, phase)
grid
xlabel('Frequency (Hz)')
ylabel('Phase (°)')

phase = deg2rad(phase);
ResponseData = mag.*exp(1i* phase);
data = idfrd(ResponseData,FrequencyVector,Ts);
H = tfest(data,21)

Warning: The frequency response of the data is not symmetric. The estimated model may have complex coefficients.

H =
                                                                                                                                                                                                                                                               
    (6.409e10 - 0.25i) s^20 + (3.043e20 - 1.4e10i) s^19 + (1.778e31 + 0.12i) s^18 + (9.83e40 + 1.1e-10i) s^17 + (1.888e51 - 1.5e-21i) s^16 + (9.109e60 - 8.6e-31i) s^15                                                                                        
                                                                                                                                                                                                                                                               
            + (5.606e70 + 1.3e-41i) s^14 + (2.129e80 + 7.1e-51i) s^13 + (4.829e89 - 8e-62i) s^12 + (1.75e99 - 6.1e-71i) s^11 + (-3821001066788808618254674132373508191857529213743881419397435822626074505318229650308583957089548312499978240                 
                                                                                                                                                                                                                                                               
            + 3e-82i) s^10 + (-669308385220897775822354873232736644105116927840756096139407434288190428346768678834045661822189748810905864196063232 + 3.7e-91i) s^9                                                                                           
                                                                                                                                                                                                                                                               
            + (3.746e125 - 7.7e-100i) s^8 + (2.283e132 + 2.4e-108i) s^7 + (-897522051716783999776694737879386738178028932801116838438975297801635568837388526411373842729542020636605371739273253938739191261223651901440                                      
                                                                                                                                                                                                                                                               
            + 4.3e-116i) s^6 + (2.094e147 - 1.4e-123i) s^5 + (5.364e154 - 8.2e-131i) s^4 + (-49783119306763802095808224367901859260341395168623567100591006380047610427499421780990029208444881630627634403016937243655815973582150348026484740435906178056192 
                                                                                                                                                                                                                                                               
            + 2.3e-137i) s^3 + (-75857188497597857641298616451815620992988776266525793680825342057338430735307997872752682463363096649544153256692252106381378394501049906334132963078207179534383448064                                                       
                                                                                                                                                                                                                                                               
                                                                                                       + 4e-144i) s^2 + (1.543e172 - 6.6e-150i) s + (1.227e177 - 1.4e-157i)                                                                                    
                                                                                                                                                                                                                                                               
  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                                                                                                                                 
  s^21 + (1.223e09 + 0.17i) s^20 + (8.384e20 - 1.8e08i) s^19 + (6.898e29 + 8.8i) s^18 + (1.775e41 + 3.1e-11i) s^17 + (1.234e50 - 7e-20i) s^16 + (1.389e61 + 6.5e-31i) s^15                                                                                       
                                                                                                                                                                                                                                                                 
          + (8.477e69 + 5.7e-40i) s^14 + (3.407e80 - 6.4e-51i) s^13 + (8.877e88 - 4.8e-60i) s^12 + (3.249e99 + 5.2e-71i) s^11 + (-1388136594865208938608928947659574094111170936969782633967986657404914106863655121395937700789486069101887488                  
                                                                                                                                                                                                                                                                 
          + 4.1e-80i) s^10 + (-1824121275203458022646298360623905608391186090971718326907382319479084925883460448103034817384106076041825180595519488 - 4.9e-91i) s^9                                                                                            
                                                                                                                                                                                                                                                                 
          + (2.087e125 + 4.7e-100i) s^8 + (2.219e133 + 8.8e-109i) s^7 + (-670081589406671093038378258373967414658263607962424697683727645451983707849333914896921690967280874069435439309545058239155399812566401155072                                          
                                                                                                                                                                                                                                                                 
          - 5.6e-116i) s^6 + (-6297844119119924774018922148959963918398164054359256773015625618332403386627824188484545592819537632977248323057506999589316494855220126351389360128                                                                              
                                                                                                                                                                                                                                                                 
          + 1e-124i) s^5 + (5.053e154 + 1.4e-130i) s^4 + (3.775e160 - 7.7e-138i) s^3 + (-88364346780427700392668889710016214689873526899655100135632884632174793641741930732221871526807952030967512922595958488429902089151968750600222417171450644558805729280 
                                                                                                                                                                                                                                                                 
          - 9.3e-144i) s^2 + (-2346318975757809678853854917348473017882638240017756637765077347074317436141829431663256417368419933588065873742644185291413073370784336481693259484277180561654488449941504                                                      
                                                                                                                                                                                                                                                                 
                                                                                                                                   + 3e-150i) s + (1.718e177 + 7.4e-156i)                                                                                        
                                                                                                                                                                                                                                                                 
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 21   Number of zeros: 20
   Number of free coefficients: 42
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                           
Estimated using TFEST on frequency response data (complex) "data".
Fit to estimation data: 92.72%                                    
FPE: 0.009037, MSE: 0.006117                                      
 

N = H.Numerator

N = 
  Columns 1 through 7

   6.4092e+10 - 2.5158e-01i   3.0434e+20 - 1.3991e+10i   1.7782e+31 + 1.1955e-01i   9.8296e+40 + 1.0816e-10i   1.8877e+51 - 1.5144e-21i   9.1094e+60 - 8.6308e-31i   5.6064e+70 + 1.3336e-41i

  Columns 8 through 14

   2.1286e+80 + 7.1227e-51i   4.8291e+89 - 8.0144e-62i   1.7503e+99 - 6.0768e-71i -3.8210e+108 + 3.0240e-82i -6.6931e+116 + 3.6968e-91i  3.7464e+125 -7.6658e-100i  2.2835e+132 +2.3999e-108i

  Columns 15 through 21

 -8.9752e+140 +4.2889e-116i  2.0944e+147 -1.3612e-123i  5.3643e+154 -8.2252e-131i -4.9783e+160 +2.3487e-137i -7.5857e+166 +4.0439e-144i  1.5431e+172 -6.6041e-150i  1.2267e+177 -1.3671e-157i

D = H.Denominator

D = 
  Columns 1 through 7

   1.0000e+00 + 0.0000e+00i   1.2233e+09 + 1.7424e-01i   8.3835e+20 - 1.7649e+08i   6.8977e+29 + 8.8494e+00i   1.7746e+41 + 3.0840e-11i   1.2338e+50 - 7.0297e-20i   1.3887e+61 + 6.4778e-31i

  Columns 8 through 14

   8.4771e+69 + 5.7149e-40i   3.4066e+80 - 6.3707e-51i   8.8769e+88 - 4.7648e-60i   3.2486e+99 + 5.1790e-71i -1.3881e+108 + 4.0517e-80i -1.8241e+117 - 4.9119e-91i  2.0871e+125 +4.7019e-100i

  Columns 15 through 21

  2.2192e+133 +8.8322e-109i -6.7008e+140 -5.6408e-116i -6.2978e+147 +1.0444e-124i  5.0531e+154 +1.4264e-130i  3.7748e+160 -7.6706e-138i -8.8364e+166 -9.2685e-144i -2.3463e+171 +2.9519e-150i

  Column 22

  1.7183e+177 +7.3614e-156i

figure
compare(data, H)
grid

% [num]=[-1.26e+07+1.83e+04i -1.84e+13 - 4.98e+11i];%taken from H.Numerator
% [den]=[1+ 0i -2.04e+06+8.94e+04i -2.9e+13-1.08e+11i];%taken from H.Denominator
% Ts=0.1e-6;
% sys1=tf(num,den,Ts);

figure
bode(H)
grid

% t=0:0.1e-6:10;

figure
pzmap(H)
grid

figure
nyquist(H)
grid

[y,t]=step(H);
SI = stepinfo(H)

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

SI = struct with fields:
         RiseTime: NaN
    TransientTime: NaN
     SettlingTime: NaN
      SettlingMin: NaN
      SettlingMax: NaN
        Overshoot: NaN
       Undershoot: NaN
             Peak: Inf
         PeakTime: Inf

figure
plot(t,real(y), 'DisplayName','Re')
hold on
plot(t, imag(y), 'DisplayName','Im')
plot(t, abs(y), 'DisplayName','Abs')
hold off
grid
legend('Location','bestoutside')
ylim([-1 1]*100)

.

[Ganesh Prasad]

@Star Strider:

Thanks a lot for taking your time out and explaining in detail the best possible solution to the issue.

Your solution has helped me with my work immensely.

[Star Strider]

As always, my pleasure!

I am happy to have been able to help you!

[Ganesh Prasad]

@Star Strider

Hello,I am facing an issue related to IFFT.The problem is in the ink:

https://in.mathworks.com/matlabcentral/answers/2085218-determining-the-input-from-the-output-using-ifft?s_tid=srchtitle

Can you please help me in correcting the code.Thanks.

[Star Strider]

I looked at the link in your Comment when it was posted. I cannot help with it.

The inverse Fourier transform of a complex double-sided frequency response (such as that produced by the fft function without using the fftshift function) should reproduce the time-domain response exactly without needing any sort of other filtering or post-processing. If you only have the single-sided complex Fourier transform (or can synthesize it from the magnitude and phase information, denoted as ‘Fourier1’ here), you can approximate the two-sided Fourier transform using:

Fourier2 = [Fourier1, flip(conj(Fourier1))];

time_domain = ifft(Fourier2, 'symmetric');

This assumes row-vectors. Column vectors require the semicolon concatenation operator in the ‘Fourier2’ calculation instead of the comma operator.

As for determining the input given only the output, that requires also knowing the system that created the output. It is necessary to have any two items of information (input signal, plant or filter, output) in order to determine the third. Time domain information is best, although it is possible to use frequency domain information to get some information about the plant or filter.

[Ganesh Prasad]

@Star Strider Ok.Thanks for replying to the question.

[Star Strider]

My pleasure.

Apparently that has some system defined in it, however I was unable to determine what the system was, or how it was defined. If it is a plant or filter, you may be able to get the information you want from it.

[Ganesh Prasad]

The system(A probe in my case) is defined in the form of Scattering parameter(S2P) matrix.

The theory i read in the internet is that impulse response is the IFFT of S21 vector(transmission coefficients) and the input is determined by the convolution of output with derived inverse impulse response.

However, that theory is not working in my case.

Answer 2

[Aquatris]

First of all, bode command does not give you numerator and denominator. It gives you magnitude and phase of the system, or FRD in your case.

In order to get a step response with the 'step' function, you need to provide a transfer function or state space representation of a system to the step function.

That is why you first need to fit a model to your FRD. You can use the system identification toolbox if you have access to it. Alternatively, if you know the structure of your model, you can also formulate a simple optimization problem or manually try to fit a model by adjusting the parameters.

Once you fit a model to your FRD, then you can use the step function to obtain the step response.