**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Transfer function extraction from frequency response

27 views (last 30 days)

Show older comments

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.

### Accepted Answer

Star Strider
on 20 Feb 2024

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.

.

##### 17 Comments

Star Strider
on 20 Feb 2024

Ganesh Prasad
on 20 Feb 2024

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);

Star Strider
on 20 Feb 2024

Edited: Star Strider
on 20 Feb 2024

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
on 20 Feb 2024

Star Strider
on 20 Feb 2024

Ganesh Prasad
on 21 Feb 2024

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
on 21 Feb 2024

Edited: Aquatris
on 21 Feb 2024

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
on 21 Feb 2024

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
on 21 Feb 2024

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
on 22 Feb 2024

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.

Ganesh Prasad
on 26 Feb 2024

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

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

Star Strider
on 26 Feb 2024

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.

Star Strider
on 27 Feb 2024

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
on 27 Feb 2024

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.

### More Answers (1)

Aquatris
on 20 Feb 2024

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.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)