TFESTIMATE

Purpose:

Estimates the transfer function using the Welch method.

Syntax:

TFESTIMATE(x, y, win, overlap, nfft, fs, detrend, zeropad, range, output, estimator)

Returns:

A complex series or array, the complex transfer function estimate.

Example:

W1: {1, -3, 4, 6, 2}

W2: gnorm(10000, 1)

W3: conv(w1, w2)

W4: tfestimate(w2, w3)

W1 contains the impulse response of the system.

W2 contains 10000 samples of normally distributed random noise.

W3 computes the response of the system in W1 to the random noise in W2.

W4 computes the one-sided complex transfer function estimate by dividing W2 and W3 into 8 overlapping segments multiplied by a Hamming window and averaging the FFTs of the segments.

Example:

W1: {1, -3, 4, 6, 2}

W2: gnorm(10000, 1)

W3: conv(w1, w2)

W4: tfestimate(w2, w3, "twosided")

W5: real(ifft(w4));extract(w0, 1, length(w1))

Same as above, except the twosided complex transfer function estimate is computed in W4.

W5 computes the impulse response from the transfer function estimate. The estimated values compare with the actual impulse response values in W1.

Example:

W1: {1, -3, 4, 6, 2}

W2: gnorm(10000, 1)

W3: conv(w1, w2)

W4: tfestimate(w2, w3, hamming(128), 64, 1024, "twosided")

W5: real(ifft(w4))

Same as above, except the input and response series are divided into segments of 128 samples overlapped by 64 samples. Each segment is multiplied by a Hamming window and the 1024 point FFTs are averaged.

W5 computes the impulse response from the transfer function estimate. The estimated values compare with the actual impulse response values in W1.

Example:

W1: {1, -3, 4, 6, 2}

W2: gnorm(10000, 1)

W3: conv(w1, w2)

W4: tfestimate(w2, w3, hamming(128), 64, 1024, "twosided")

W5: real(ifft(w4))

W6: mag(fft(w1, 10000))

W7: mag(w4)

W8: W6;overp(W7);h=legend(0, 0, "Actual","Estimated");h.target=PAPER;h.location=5;h.orientation=1

Same as above, except W6 and W7 display the magnitude of the actual and estimated transfer functions. W8 visually compares the results.

Remarks:

Given a perfect, noise-free system where the Fourier transform of the input x[n] is X( f ) and the Fourier transform of the output y[n] is Y( f ), the transfer function of the system, H( f ) is:

However, real world systems are almost always affected by noise. Simply dividing the Fourier transforms can lead to significant errors particularly where the signal power is low.

Two widely used estimators that account for noise in the input and output series are derived from the cross-power spectral density (CPSD) and the power spectral density (PWELCH) functions.

The H₁ Estimator is considered optimal when the output noise is uncorrelated with the input signal. It is the most common choice in practice, as it assumes the noise primarily affects the measurement of the output. It is calculated as the ratio of the cross-power spectral density of the output and input, S_yx( f ), to the power spectral density of the input, S_xx( f ):

The H₂ Estimator is more suitable when the input noise is uncorrelated with the output signal. It is calculated as the ratio of the power spectral density of the output S_yy( f ) to the cross-power spectral density of the input and output S_xy( f ):

TFESTIMATE estimates the complex transfer function between an input series x and a corresponding output series y using the Welch method (PWELCH) to estimate the cross-power spectral density and the power spectral density. Both signals are segmented into overlapping windows, each windowed and transformed using the FFT. The averaged spectral estimates are then used to calculate the transfer function according to the H₁ or H₂ estimation method, as specified by estimator option. By convention, the complex conjugate of x is used to compute the H₁ transfer function.

If nfft is an array of one or more target frequencies, TFESTIMATE uses the GOERTZEL algorithm to compute the transfer function at each input frequency. The input frequencies can be any real values. Frequency values less than zero or greater than rate(x) are wrapped to the range between 0 and rate(x).

If detrend is "constant" or "linear", the mean or linear trend is computed and removed from each segment before FFT processing.

If zeropad is "zeropad", the last column is padded with zeros to make it the same length as the other columns. If zeropad is "none" or empty, the last column is ignored if it is shorter than the others.

If range is "twosided", the TFESITMATE produces the full transfer function from 0 to rate(x) in wrap-around order. For example, a 10 point full transfer function with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5}

If range is "center", the full transfer function is displayed from -rate(x)/2 to rate(x)/2 with the zero frequency in the middle. For example, a 10 point full transfer function with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{0, 1, 2, 3, 4, 5, -4, -3, -2, -1}

If both inputs are real, range defaults to "onesided", else if any input is complex, range defaults to "twosided".

By default, the output of TFESTIMATE is complex. Use MAGNITUDE and PHASE to compute the magnitude and phase components from the complex transfer function. For example, the magnitude-squared of the transfer function is computed with:

mag(TFESTIMATE(x, y))^2

and phase of the transfer function is computed with:

phase(TFESTIMATE(x, y))

If output is "mag2" or "mag", TFESTIMATE returns the magnitude-squared or the magnitude of the transfer function directly.

See CPSD to estimate the complex cross-power spectral density.

See MSCOHERE to estimate the magnitude-squared coherence function.

See PWELCH to estimate the power spectral density.

See DADiSP/FFTXL to optimize the underlying FFT computation.

See Also:

CPSD

FFT

GOERTZELDFT

HAMMING

HANNING

KAISER

MSCOHERE

PSD

PWELCH

SPECTRUM

References: