Transfer function extraction from frequency response

Question

Ganesh Prasad 2024-2-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2084333-transfer-function-extraction-from-frequency-response

评论： Ganesh Prasad 2024-2-27

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2024-2-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2084333-transfer-function-extraction-from-frequency-response#answer_1412673

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.

.

17 个评论
显示 15更早的评论隐藏 15更早的评论

Star Strider 2024-2-20

编辑：Star Strider 2024-2-20

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.

Star Strider 2024-2-21

在 MATLAB Online 中打开

Try4_TF.xlsx

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)

.

Star Strider 2024-2-26

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.

Star Strider 2024-2-27

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 2024-2-27

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 2024-2-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2084333-transfer-function-extraction-from-frequency-response#answer_1412668

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.