Getting an equation from a signal transfer function

60 visualizzazioni (ultimi 30 giorni)
Hi guys, I dont know if this is possible or not, but I have two audio signals, an input and an output, I then got the transfer function of those two signals using fft, but now I would like to get a mathematical expression for that transfer function, do you guys know of anyway I can achieve this ?
  6 Commenti
Sifiso Mzobe
Sifiso Mzobe il 14 Ott 2023
how would i go about this? my signal is continious time signal.
Walter Roberson
Walter Roberson il 15 Ott 2023
audio signals are always discrete time, never continuous time.

Accedi per commentare.

Risposte (4)

Walter Roberson
Walter Roberson il 14 Ott 2023
Modificato: Walter Roberson il 15 Ott 2023
Your formula for calculating the transfer function is not correct. Observe this simple example:
syms t
f1 = sin(t)
f1 = 
f2 = sin(2*t)
f2 = 
L1 = laplace(f1)
L1 = 
L2 = laplace(f2)
L2 = 
L12 = simplify(L1./L2)
L12 = 
That is the transfer function
f12 = ilaplace(L12)
f12 = 
That is the time-domain of the transfer function.
Now let us try numerically:
T = linspace(0,2*pi,128);
y1 = double(subs(f1, t, T));
y2 = double(subs(f2, t, T));
Y1 = fft(y1);
Y2 = fft(y2);
Y12 = Y1./Y2;
n12 = ifft(Y12);
plot(T, y1, T, y2); legend("numeric " + string([f1 f2]))
figure()
fplot(L12, [0 2*pi]); legend("symbolic transform")
figure()
subplot(2,1,1); plot(T, real(n12)); legend("numeric transform real");
nnz(isfinite(real(n12)))
ans = 0
subplot(2,1,2); plot(T, imag(n12)); legend("numeric transform imaginary")
The real plot is blank because none of the coefficients are finite.
Consider the discrete fourier transform of a pure sine wave: provided that it is sampled finely enough, it will be non-zero only at the bin corresponding to the primary frequency. Consider then a pure sine wave of a different frequency: the fft for it will be non-zero only at the bin corresponding to the primary frequency of the second wave. But they are different primary frequencies, so the non-zero bin in the second one is in a different location than the non-zero bin for the other. Divide the fft's and you will get a lot of 0/0 entries at the places that both are 0, and you will get infinite entries at the bin that the first frequency is non-zero, and you will get 0 entries at the bin that the second frequency is non-zero (because it would be 0 for the first fft divided by something finite.) Lots of nan, two infinities, two zeros. ifft() that and the nan "pollute" the entire calculation and the result is all nan. This is clearly not a meaningful transfer function.
Is there a meaningful transfer function in this case? There certainly is, and it is easy to find using laplace transform.
  3 Commenti
Walter Roberson
Walter Roberson il 15 Ott 2023
and then how do i get n12 as an equation?
In the case of this example,
n12 = @(t) nan(size(t))
This would also apply for the case where the fft of the second signal is zero anywhere
Walter Roberson
Walter Roberson il 15 Ott 2023
In this below example we can see that the ratio of fft gives a completely different time domain function than if you work using laplace transform.
I am pretty sure your ratio-of-fft is wrong.
syms t
f1 = sin(2*t)
f1 = 
f2 = t^2 - 1
f2 = 
fplot([f1, f2], [0 pi])
L1 = laplace(f1)
L1 = 
L2 = laplace(f2)
L2 = 
L12 = simplify(L1/L2)
L12 = 
f12 = ilaplace(L12)
f12 = 
char(f12)
ans = '- (4*cos(2*t))/3 - (2*cosh(2^(1/2)*t))/3'
f12s = simplify(f12)
f12s = 
fplot(f12s, [0 2*pi])
T = linspace(0,2*pi,128);
y1 = double(subs(f1, t, T));
y2 = double(subs(f2, t, T));
Y1 = fft(y1);
Y2 = fft(y2);
Y12 = Y1./Y2;
n12 = ifft(Y12);
plot(T, y1, T, y2); legend("numeric " + string([f1 f2]))
figure()
subplot(2,1,1)
fplot(L12, [0 2*pi]); legend("symbolic transform")
subplot(2,1,2)
fplot(f12, [0 2*pi]); legend("symbolic transform in time domain")
figure()
subplot(2,1,1); plot(T, real(n12)); legend("numeric transform real");
subplot(2,1,2); plot(T, imag(n12)); legend("numeric transform imaginary")

Accedi per commentare.


Star Strider
Star Strider il 15 Ott 2023
Modificato: Star Strider il 16 Ott 2023
I would like to get a mathematical expression for that transfer function
If you have access to the System Identification Toolbox, start with the iddata function and then either tfest or ssest, depending on what you want. Use the compare function to view the results of the estimation. Tweak the estimation function call order argument until you get the appropriate order that fits best.
The System Identification Toolbox functions are the mosst robust, however the Signal Processing Toolbox has the invfreqz and freqz functions (or their continuous equivalents) that can produce similarl results.
EDIT — (16 Oct 2023 at 02:03)
One option that I had not mentioned is that you can estimate the poles and zeros and their locations from the imaginary part of the transfer function, plotted as a function of frequency. The poles will be at the frequencies where the imaginary part of the transfer fucntion has singularities, and the zeros will be the zero-crossings between the poles. There may be a pole or zero at the origin and infinity as well.

Paul
Paul il 15 Ott 2023
Modificato: Paul il 15 Ott 2023
Hi Sifiso,
Here is an example that might be helpful.
Generate an input signal
inputAudio = sweeptone(6,0.1);
Fs = 44.1e3;
N = numel(inputAudio)
N = 269010
t = (0:N-1)/Fs;
DFT of the input
U = fft(inputAudio);
Frequency vector (rad/sec)
w = (0:N-1)/N*2*pi*Fs;
Transfer function for example that we'll want to reconstruct. I'm assuming that the audio is going through an analog device.
w0 = 2.5e3;
h = tf(w0^2,[1 2*.3*w0 w0^2])
h = 6.25e06 ---------------------- s^2 + 1500 s + 6.25e06 Continuous-time transfer function.
Generate the outputAudio as the ouput of the transfer function excited by the input.
outputAudio = lsim(h,inputAudio,t);
DFT of output
Y = fft(outputAudio);
Ratio of DFT's output/input
H = Y./U;
Bode plot of H
bode(frd(H(1:N/2),w(1:N/2)))
The Bode plot looks reasonable at all frequencies (up to the Nyquist frequency), so we can use all of the data to estimate the transfer function. Note that the inputs is designed to excite the system at all frequencies. Other inputs may have different characteristics. Anyway ...
Estimate the transfer function using tfest, assume two poles and don't specify the number of zeros. Use only the data up to the Nyquist frequency.
sys1 = tfest(frd(H(1:N/2),w(1:N/2)),2)
sys1 = -69.19 s + 6.243e06 ----------------------- s^2 + 1499 s + 6.246e06 Continuous-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 1 Number of free coefficients: 4 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on frequency response data. Fit to estimation data: 99.75% FPE: 2.877e-07, MSE: 2.876e-07
Same thing, but specify no zeros.
sys2 = tfest(frd(H(1:N/2),w(1:N/2)),2,0)
sys2 = 6.143e06 ----------------------- s^2 + 1475 s + 6.145e06 Continuous-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 0 Number of free coefficients: 3 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on frequency response data. Fit to estimation data: 98.02% FPE: 1.83e-05, MSE: 1.83e-05
tfest can also be used directly on the time domain data
sys3 = tfest(iddata(outputAudio,inputAudio,1/Fs),2,0)
sys3 = From input "u1" to output "y1": 6.25e06 ---------------------- s^2 + 1500 s + 6.25e06 Continuous-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 0 Number of free coefficients: 3 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data. Fit to estimation data: 100% FPE: 4.421e-24, MSE: 4.421e-24
Compare the DFT ratio, the true transfer function, and the estimated transfer functions.
hplot = bodeplot(frd(H(1:N/2),w(1:N/2)),h,sys1,sys2,sys3,w(1:N/2));
setoptions(hplot,'PhaseWrapping','on')

William Rose
William Rose il 16 Ott 2023
There are two distinct threads in this discussion.
One thread is an answer to your original quesiton: 'I measured a transfer function; how do I find a mathematical function that fits it?" @Star Strider gave an excellent answer to that quesiton, and I gave an answer in my comment.
The other thread is "How does one measure (or, to be more accurate, estimate) a transfer funciton?" which was not part of your original question. This thread arose because you said you measured the TF by fft(y)/fft(x). @Walter Roberson gave good examples showing that this method of estimation can give completely wrong results, and he said "I am pretty sure your ratio-of-fft is wrong." (More on that later.) @Paul gave an example where estimating H(f) with , where Y(f)=fft(y(t)) and X(f)=fft(x(t)), gives a reasonable answer. @Paul also showed how you can use tfest() to get a mathematical function for H(f), which works well for the second order lowpass system he used.
@Walter Roberson is correct that the ratio-of-fft method is a bad idea. The ratio-of-fft method worked in the example of @Paul, because the signal x(t) was a swept sine wave, so the denominator fft was not zero, or close to zero, at any frequency. But for many signals, that will not be the case. And in the example of @Paul, there was no noise on x(t) or y(t).
I tried the same thing as you - the ratio-of-fft method - when I first tried to measure transfer functions (between blood flow and blood pressure) experimentally. The results were unsatisfactory, even when I used an external pump to add white noise to the blood flow. (This was to assure the denominator fft was never zero, or close to zero.) I did some reading, and then implemented an estimation method that is statistically superior. Then I got good results, using the experimental data I had already collected.
@Star Strider and @Paul recommend tfest(). It will implement a good method of estimation. It has many options which you will have to read about in the help.
I am attaching some notes on the topic.
  4 Commenti
Paul
Paul il 24 Ott 2023
Hi William,
I'm not sure what's all that different between your example and the example below, other than that this example uses a much higher sampling frequency and many more points. Maybe that makes a lot of difference. For this example, there really isn't that much difference between using a simple FFT ratio and tfestimate, which I think is similar, at least in spirit, with your approach, even with a large amount of noise, or only noise, at the input. I didn't look at what happens when noise is added at the output. In the original question, the goal was to derive a model, and the differences in the derived models are not visible to my eye (using the stock Bode plotting functions, not zooming in).
I'm not advocating that the FFT ratio is always the right approach to generate the data to estimate the model, just showing that in some cases it might not be so bad, even if tfestimate may be better.
I do have onne question about the cross PSD approach that you showed and is used in tfestimate. It's my understanding that the psd and cpsd estimates assume that the input signal is stationary. But the sweeptone signal is not stationary, though tfestimate seems to work fine anyway in this example. Are there ever any concerns using cpsd/psd approach for non-stationary signals?
inputAudio = sweeptone(6,0.1);
Fs = 44.1e3;
rng(100);
noise = 1*randn(size(inputAudio)); % one sigma value is 2x larger than the input amplitude
analyze(inputAudio,Fs);
analyze(inputAudio + noise,Fs);
analyze(noise,Fs);
function analyze(u,Fs)
N = numel(u);
alpha = 0.8;
h = tf(1-alpha,[1 , -alpha],1/Fs,'Variable','z^-1');
y = lsim(h,u);
U = fft(u);
Y = fft(y);
Hfft = Y./U;
w = (0:N-1)/N*2*pi*Fs;
w = w(1:N/2);
Hfft = frd(Hfft(1:N/2),w,1/Fs);
[Htfe,wtfe] = tfestimate(u,y);
Htfe = frd(Htfe,wtfe*Fs,1/Fs);
figure
hplt = bodeplot(Hfft,Htfe);
title('Estimated Transfer Functions')
setoptions(hplt,'PhaseWrapping','on')
sys1 = tfest(Htfe,1,'Ts',1/Fs);
sys2 = tfest(Hfft,1,'Ts',1/Fs);
figure
bode(sys2,sys1,h,{1 1e6});
title('Estimated Models v True Model')
end
William Rose
William Rose il 24 Ott 2023
Great points. I don't have time to reply appropriately at the moment. I hope to do so next week.

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!

Translated by