Getting an equation from a signal transfer function

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 ?

[Walter Roberson]

How are you calculating the transfer function by using fft() ?

[Sifiso Mzobe]

TF = fft(audio1)./fft(audio2);

[Sifiso Mzobe]

or better put, TF = fft(output_audio)./fft(input_audio);

[William Rose]

@Sifiso Mzobe,

I have never seen a direct answer to your question. Depending on what the transfer funciton looks like, you could try using one of the classic filters, or a cascaded set of them. Or you could fit a rational polynomial in s=iw (continuous time), or a rational polynomial in z=exp(iw) (discrete time).

[Sifiso Mzobe]

how would i go about this? my signal is continious time signal.

[Walter Roberson]

audio signals are always discrete time, never continuous time.

Answer 1

[Walter Roberson]

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.

[Sifiso Mzobe]

and then how do i get n12 as an equation?

[Walter Roberson]

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]

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

Answer 2

[Star Strider]

‘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.

Answer 3

[Paul]

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

Answer 4

[William Rose]

@Sifiso Mzobe,

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.

[Paul]

I think that the statement "that the ratio-of-fft method is a bad idea" is too broad. It's certainly a simple approach, and it might not be the best approach to blindly use for all inputs, but that doesn't mean it's a bad idea in general. As shown, it can work well for swept sine waves, or chirps when focusing on specific frequencies. Doesn't tfestimate essentially use this approach (though with more advanced proessing in the frequency domain)? To be sure, one has to understand the properties of the input signal and not do the processing blindly.

Perhaps it would be instructive to show an example where the simple ratio-of-fft does not work (for a stable system with impulse response that has a Fourier transform).

[William Rose]

@Paul,

Here is a script that does what you requested: show an example where the simple ratio-of-fft does not work (for a stable system with impulse response that has a Fourier transform).

The input x(k) is Gaussian and white. The output y(k) is a first order lowpass filtered version of x(k):

The Z transform of the equation above is

and therefore the transfer function is

or

This filter has unity gain at DC. alpha is in the range (0,+1). When alpha=0, there is no filtering, and y(k)=x(k). As alpha increases, the lowpass filtering gets stronger.

The attached script estimates the transfer function H(f) by four methods, and it displays the theoretical transfer function:

H0=fft(y)./fft(x)

H1=Sxy./Sxx

H2=Syy./Syx

Hv=(H1.*H2).^0.5

where

[Sxy,f]=cpsd(y,x,window,noverlap,nfft,fs)

[Syx,f]=cpsd(x,y,window,noverlap,nfft,fs)

[Sxx,f]=cpsd(x,x,window,noverlap,nfft,fs)

[Syy,f]=cpsd(y,y,window,noverlap,nfft,fs)

Matlab's tfestimate() estimates the transfer function by equation H1 above, by default.

The script produces output such as below, when there is zero measurement noise on x and y. Even in this idealized case, it is clear that the estimate H0=fft(y)/fft(x) is very noisy compared to the other estimates. When measurement noise is added, the inferioty of H0=fft(y)/fft(x) is even more apparent.

Zoom in on the range f=0 to 50 Hz:

It has been shown that H1 is the optimal solution (minimum squared error) when there is measurement noise on y(t) and no measurement noise on x(t).

It has been shown that H2 is the optimal solution (minimum squared error) when there is measurement noise on x(t) and no measurement noise on y(t).

Vold et al (1984) proposed Hv as a compromise when there is measurement noise on both x(t) and y(t).

You can adjust the amount of measurement noise on x and y by changing the values of th constants xNoiseAmp and yNoiseAmp in the script. This allows you to "experimentally" test the statements above about H1, H2, Hv.

By the way: The Matlab function Pxy=cpsd(x,y,...) returns the expected value of X.conj(Y). Authors 1-4 (and probably other authors) define Pxy or Sxy as the expected value of conj(X).Y, i.e. they have the opposite convention about which gets the complex conjugate. This could be a source of confusion. In my code, and above, I use Sxy to denote the expected value of conj(X).Y, which is consistent with authors 1-4.

1. Jenkins & Watts (1968), Eq.8.3.3, p.341.

2. Bendat & Piersol (1971), Eqs.3-90,3-100, pp.82,84.

3. Otnes & Enochson (1972), Eq.7.33, p.304.

4. Bendat & Piersol (1980), Eqs.3.45,3.46, pp.54-55.

[Paul]

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]

@Paul,

Great points. I don't have time to reply appropriately at the moment. I hope to do so next week.