View Raw SPL
/*****************************************************************************
* *
* TFESTIMATE.SPL Copyright (C) 2024 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Welch transfer function estimate *
* *
* Revisions: 11 Jul 2024 RRR Creation *
* *
*****************************************************************************/
#include
#if @HELP_TFESTIMATE
TFESTIMATE
Purpose: Estimates the transfer function using the Welch method.
Syntax: TFESTIMATE(x, y, win, olag, nfft, fs, detrend, zeropad, range)
x - A series or array, the input.
y - A series or array, the response.
win - An integer or series, the segment size or window.
If a series, each input segment is multiplied by
WIN with the segment size set to LENGTH(win).
Defaults to a Hamming window of where the length
is 1/8 the length of the input series.
olap - Optional. An integer, the number of points to
overlap. Defaults to length(WIN)/2.
nfft - Optional. An integer or series. If NFFT is an
integer, NFFT is the FFT length. If NFFT is a
series, the series contains the frequencies used
to evaluate the transfer function estimate.
Defaults to the integer bestpow2(WIN).
fs - Optional. A real, the input series sample rate.
Defaults to RATE(s).
detrend - Optional. A string, the segment detrend mode:
"none" : no detrending (default)
"constant" : remove the mean of the segment
"linear" : remove the linear trend of the
segment
zeropad - Optional. A string, the short segment zero pad mode:
"none" : no zero padding, ignore the last column
if it is shorter than the other columns
(default)
"zeropad" : add zeros to the last column if it is
shorter than the other columns to force
the same length for all columns
range - Optional. A string, the frequency range of the
transfer function estimate.
"onesided" : one sided transfer function estimate
from 0 to Fs/2 (default)
"twosided" : full transfer function estimate
from 0 to Fs
"center" : centered transfer function estimate
from -Fs/2 to Fs/2
Returns: A complex series or array, the 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 contains the response of the system in W1 to the random noise
in W2.
W4 computes the onesided 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))
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.location=1
Same as above except W6 and W7 display the magnitude of the actual
and estimated transfer functions. W8 visually compares the results.
Remarks:
TFESTIMATE estimates the complex transfer function between
an input series X and a corresponding output series
Y using the Welch method (see PWELCH). The input and output
series are divided into overlapping segments and the FFTs
of the windowed segments are averaged.
The output of TFESTIMATE is complex.
By convention, the complex conjugate of X is used to compute
the transfer function estimate.
The magnitude squared of the transfer function is determined by:
mag(TFESTIMATE(x, y)^2)
The phase of the transfer function:
phase(TFESTIMATE(x, y))
If NFFT is an array of one or more target frequencies, TFESTIMATE
uses the Goertzel algorithm to compute the TFESTIMATE magnitudes
at each input frequency. The input frequencies can be any real
values. Frequency values less than zero or greater than rate(S)
are wrapped to the range between 0 and rate(S).
If DETREND is "constant" or "linear", the mean or linear
trend is computed and removed from each segment before FFT
processing.
By default, the RANGE value is "onesided" and transfer function
estimate is displayed from 0 to rate(X)/2.
If the RANGE value is "twosided", TFESTIMATE produces the full
transfer function from 0 to rate(X) in wrap around order. For
example, a 10 point full transfer function estimate 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 RANGE is "center", the full transfer function estimate 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
estimate 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}
See Also:
Cpsd
Hamming
Hanning
Kaiser
Mscohere
Psd
Pwelch
References:
[1] Oppenheim & Shafer
Digital Signal Processing
Prentice Hall, 1975
pp. 548-556
[2] Oppenheim & Shafer
Discrete-Time Signal Processing
Prentice Hall, 1989
pp. 730-742
[3] Programs for Digital Signal Processing
IEEE Press, 1979
Section 2.1-1 - 2.1-10
[4] Press, Flannery, Teukolsky, Vetterling
Numerical Recipes in C
Cambridge Press, 1988
pp. 407-552
[5] S. Lawrence Marple, Jr.
Digital Spectral Analysis with Applications
Prentice Hall, 1987
pp. 152-158
#endif
/* welch transfer function */
ITERATE tfestimate(x, y, win, olap, nfft, fs, detrend, zeropad, range, output = "complex", est = "")
{
local pxy = {}, pxx = {}, pyx = {}, pyy = {}, tf;
if (argc < 2 || not(isarray(x)) || not(isarray(y)))
{
error(sprintf("%s - two input series required", __FUNC__));
}
/* parse arguments */
(win, olap, nfft, fs, dflag, range, output, est) = welch_parse_args(x, y, win, olap, nfft, fs, detrend, zeropad, range, output, est, __FUNC__);
if (est == "h2")
{
/* cross PSD */
pyx = cpsd(y, x, win, olap, nfft, fs, detrend, zeropad, range, output);
/* input PSD */
pyy = pwelch(y, win, olap, nfft, fs, detrend, zeropad, range);
/* transfer function */
tf = pyy / pyx;
}
else
{
/* cross PSD */
pxy = cpsd(x, y, win, olap, nfft, fs, detrend, zeropad, range, output);
/* input PSD */
pxx = pwelch(x, win, olap, nfft, fs, detrend, zeropad, range);
/* transfer function */
tf = pxy / pxx;
}
if (isxy(pxy) || isxy(pxx) || isxy(pyx) || isxy(pyy))
{
/* xy */
tf = xy(nfft, tf);
sethunits(tf, gethunits(pxy));
}
return(tf);
}