1
$\begingroup$

I have been trying to replicate the results from the tfestimate function in Matlab by performing, from my understanding, the corresponding signal processing steps in a separate Matlab script file. The results simply differ.

The ultimate goal is to be able to run a tfestimate-like replica in C/C++.

Anyhow, my input signal $x(n)$ being injected into a certain plant is a PRBS of order $\mathcal{O}=15$, length $N=2^{\mathcal{O}}-1$, shifted and scaled onto $\\{-a,a\\}$ where $0 < a < 1$ rad/s. The output signal $y(n)$, unit rad/s, contain readings from a rategyro. In- and output signals are sampled at $f_s=128$ Hz.

Algorithm

My processing steps go like this,

  1. Calculate auto- and cross-correlations, \begin{equation} r_{xx}(l)=\sum_{n=0}^{N-1}x(n)x^{\ast}(n-l)\textrm{ , }l=0,\pm 1,\ldots,\pm (N-1) \end{equation} \begin{equation} r_{yx}(l)=\sum_{n=0}^{N-1}y(n)x^{\ast}(n-l) \textrm{ , }l=0,\pm 1,\ldots,\pm (N-1) \end{equation} where $\ast$ denotes complex conjugate.
  2. Divide $r_{xx}(l)$ and $r_{yx}(l)$ into data segments of length $M=1024$ with $50\\%$ overlap between successive data segments. Doing so yield $L=128$ segments of length $M$ for $r_{xx}(l)$ and $r_{yx}(l)$ respectively.
  3. Let us denote the segments from 2. $r_{xx}^{(i)}(n)$ and $r_{yx}^{(i)}(n)$ where $n=0,1,\ldots,M-1$ and $i=0,1,\ldots,L-1$.
  4. Select a windowing function $w(n)$ e.g. Hamming where $n=0,1,\ldots,M-1$.
  5. Calculate modified periodograms, \begin{equation} P_{xx}^{(i)}(f)=\dfrac{1}{MU}\sum_{n=0}^{M-1}r_{xx}^{(i)}(n)w(n) e^{-j 2\pi f n}=\dfrac{1}{MU}X^{(i)}(f)\textrm{ , }i=0,1,\ldots,L-1 \end{equation} \begin{equation} P_{yx}^{(i)}(f)=\dfrac{1}{MU}\sum_{n=0}^{M-1}r_{yx}^{(i)}(n)w(n) e^{-j 2\pi f n}=\dfrac{1}{MU}Y^{(i)}(f)\textrm{ , }i=0,1,\ldots,L-1 \end{equation} using FFT of sequence $r_{xx}^{(i)}(n)w(n)$ and FFT of sequence $r_{yx}^{(i)}(n)w(n)$ and where, \begin{equation} U=\dfrac{1}{M}\sum_{n=0}^{M-1}w(n)^2 \end{equation}
  6. Calculate the Welch power spectrum estimates, \begin{equation} P_{xx}(f)=\dfrac{1}{L}\sum_{n=0}^{L-1}\tilde{P}^{(i)}_{xx}(f) \end{equation}

\begin{equation} P_{yx}(f)=\dfrac{1}{L}\sum_{n=0}^{L-1}\tilde{P}^{(i)}_{yx}(f) \end{equation}

  1. Calculate the transfer function estimate, \begin{equation} G(f)=\dfrac{P_{yx}(f)}{P_{xx}(f)} \end{equation}

Note 1 I am using (cross-)correlation in step 1. because I am unaware of other ways of calculating cross-spectrum $P_{yx}(f)$.

Note 2 The following identity can be used in step 5. \begin{equation} P_{xx}(f) = \mathcal{F}|r_{xx}(m)| = \sum_{m=-(N-1)}^{N-1}r_{xx}(m)e^{-j 2\pi f m} = \dfrac{1}{N}\left\lvert\sum_{n=0}^{N-1}x(n) e^{-j 2\pi f n}\right\rvert^2 =\dfrac{1}{N}\left\lvert X(f)\right\rvert^2 \end{equation} which is more in line with the Welch paper,Welchs method: Averaging modified periodograms.

Note 3 I am calculating correlation according to this in step 1.,

rxx = conv(x,fliplr(x))';

ryx = conv(y,fliplr(x))';

where x and y are row vectors. Actually I have written a separate function for this in order to not use the conv function of Matlab. What I am trying to say is that I got strange results, at least to me, when using the xcorr function in Matlab and hence I stayed away from it.

Results and comparison

It is clear from comparing result 1. and 2. below my algorithm is not really doing what Matlab does using its tfestimate function. It is also clear from result 3. the effect of squaring the DFTs in step 5. greatly enhances the estimate which is quite close to result 1.

  1. Using Matlab command,
    [G,f] = tfestimate(x,y,1024,[],[],fs);

Result of Matlab command [G,f] = tfestimate(x,y,1024,[],[],fs).

  1. Using steps 1.-7. above,

Result of step 1.-7.

  1. Using steps 1.-7. above but squaring DFTs in step 5. meaning, \begin{equation} P_{xx}^{(i)}(f)=\dfrac{1}{MU}\left(X^{(i)}(f)\right)^2 \end{equation} \begin{equation} P_{yx}^{(i)}(f)=\dfrac{1}{MU}\left(Y^{(i)}(f)\right)^2 \end{equation} Result of step 1.-7. but squaring DFTs in step 5.

Questions

  1. What is missing to the algorithm in step 1.-7. for it to replicate the tfestimate result better?
  2. How come squaring the DFTs in step 5. yield much better results? Squaring them breaks the identify in Note 2 at least to my understanding.
  3. Is there an other way of calculating cross-spectrum $P_{yx}(f)$ instead of taking the Fourier transform of the cross-correlation? I know Matlab offers such a function too, cpsd but that is a dead end to me right now. From what I have understood from the help description of the cpsd function it uses cross-correlation and Welch’s averaged, modified periodogram method of spectral estimation.

Any help is highly appreciated!

$\endgroup$
9
  • $\begingroup$ this might help. You can calculate the cpsd in the frequency domain exactly the way you would calculate the PSD. Instead of $X^*X$, it’s just $X^*Y$ where these are Fourier transforms of $x$ and $y$, not $r_{xx}$ and $r_{yy}$ $\endgroup$
    – Jdip
    Commented Apr 3 at 21:33
  • $\begingroup$ If you include the Matlab script I’ll take a look and can help better. You might be averaging things wrong… no way to know without seeing the code $\endgroup$
    – Jdip
    Commented Apr 3 at 21:35
  • $\begingroup$ Thanks. Yes I just realized the same thing. This should work $P_{yx} = Y^{*}X$. I suppose you would segment $x$ and $y$ first, window them and calculate the periodograms? Will give it a try tomorrow. Once again, thanks for responding. I’ll upload Matlab code should it be needed 😀 $\endgroup$
    – wolfiesax
    Commented Apr 3 at 21:43
  • $\begingroup$ that's correct! except I believe it would be $YX^*$ $\endgroup$
    – Jdip
    Commented Apr 3 at 23:33
  • 1
    $\begingroup$ Jdip: I can confirm you were right about the conjugation. I will compose a proper answer/work around to my issue later today. $\endgroup$
    – wolfiesax
    Commented Apr 4 at 11:28

1 Answer 1

1
$\begingroup$

I can confirm the following algorithm successfully replicates the tfestimate function for my signals $x(n)$ and $y(n)$,

Algorithm

  1. Divide $x(k)$ and $y(k)$, where $k=0,1,\ldots,N-1$, into data segments of length $M=1024$ with 50 % overlap between successive data segments. Doing so yield $L=64$ segments of length $M$ for $x(n)$ and $y(n)$ respectively.
  2. Let us denote the segments from 1. $x^{(i)}(n)$ and $y^{(i)}(n)$ where $n=0,1,\ldots,M-1$ and $i=0,1,\ldots,L-1$.
  3. Window $x^{(i)}(n)$ and $y^{(i)}(n)$ by a $M$-length Hamming window $w(n)$ and denote the result accordingly, \begin{equation} x^{(i)}_{w}(n) = \sum_{n=0}^{M-1}x^{(i)}(n)w(n)\textrm{, }i=0,1,\ldots,L-1 \end{equation} \begin{equation} y^{(i)}_{w}(n) = \sum_{n=0}^{M-1}y^{(i)}(n)w(n)\textrm{, }i=0,1,\ldots,L-1 \end{equation}
  4. Calculate $M$-point DFTs of $x^{(i)}_{w}(n)$ and $y^{(i)}_{w}(n)$, \begin{equation} X^{(i)}(f) = \mathcal{F}\left(x^{(i)}_{w}(n)\right)\textrm{, }i=0,1,\ldots,L-1 \end{equation} \begin{equation} Y^{(i)}(f) = \mathcal{F}\left(y^{(i)}_{w}(n)\right)\textrm{, }i=0,1,\ldots,L-1 \end{equation}
  5. Calculate modified periodograms, \begin{equation} P_{xx}^{(i)}(f) = \dfrac{1}{MU}X^{(i)}(f)\left(X^{(i)}(f)\right)^{*}\textrm{, }i=0,1,\ldots,L-1 \end{equation} \begin{equation} P_{yx}^{(i)}(f) = \dfrac{1}{MU}Y^{(i)}(f)\left(X^{(i)}(f)\right)^{*}\textrm{, }i=0,1,\ldots,L-1 \end{equation} where $*$ denotes complex conjugate and, \begin{equation} U = \dfrac{1}{M}\sum_{n=0}^{M-1}w^2(n) \end{equation}
  6. Calculate the Welch power spectrum estimates, \begin{equation} P_{xx}(f) = \dfrac{1}{L}\sum_{i=0}^{L-1}P_{xx}^{(i)}(f) \end{equation} \begin{equation} P_{yx}(f) = \dfrac{1}{L}\sum_{i=0}^{L-1}P_{yx}^{(i)}(f) \end{equation}
  7. Calculate the transfer function estimate, \begin{equation} G(f)=\dfrac{P_{yx}(f)}{P_{xx}(f)} \end{equation}

Note 1: In practice we do not have do divide by $U$ and $M$ in step 5. due to the calculation in step 7. The same goes for division by $L$ in step 6.

Results

The following figures illustrate a comparison between the results from the algorithm above and the corresponding from the tfestimate function in Matlab.

Magnitude comparison

Phase comparison

Note 2: The tfestimate function was run according to,

W         = M         ;
NO        = M/2       ;
NFFT      = Nfft      ;
Fs        = fs        ;
EST       = 'H1'      ;
FREQRANGE = 'twosided';
[Ghat1024,f1024] = tfestimate(x,y,W,NO,NFFT,Fs,'Estimator',EST,FREQRANGE);

Summary

The results show the algorithm replicates the tfestimate function quite well. There are some minor differences for magnitude, perhaps due to precision differences? As for the phase, it differs from 48 Hz and upwards but this ok for me realizing the magnitude is well damped for this frequency region. It would have been interesting though to get the correlation variant of the algorithm working but since this modified algorithm solves my problem I am simply happy sticking to it.

$\endgroup$
1
  • $\begingroup$ Cool! You can accept your own answer fyi ;) $\endgroup$
    – Jdip
    Commented Apr 5 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.