1
$\begingroup$

Before moving to the actual question, I would like to emphasize on the following points (maybe they are obvious to some of you, but I still would like to list them, since they make the difference):

  • I don't have any reliable information on the input, thus I have only the measurements of the output $y(t)$.

  • I have only 1 set of recorded data, thus no repetitive measurements are in this case available.

Objective: Identifying a reliable transfer function model for this data in the frequency domain (e.g. rational fraction polynomial with numerator $B$ and denominator $A$)

Question: After obtaining a transfer function model, how does one validate this model and thus can be assured that it does not contain any systematic/modelling errors? Please, keep in mind the above 2 points!

To my limited experience, calculating the uncertainty of the transfer function residuals will be not possible, since the missing information on the input (and thus also not having the variance of the input). A second approach that I thought of is to do a 'simple' whiteness test of the residuals (sum of squares due to errors (SSE) or the R-square both mentioned in the Matlab manual), but is this approach even significant? I mean, sure... One obtains a result, but it's complex valued due to the data ($Y(k)$) which is also complex valued. Consequently, there is always an imaginary part which makes it fail the criterion (i.e. model is ok if SSE is close to 0). Any thoughts/solutions on this? Or am I mixing things up?

Any help, advise/guidance (tips) or references that I can consult in order to solve this issue is very much appreciated! Because I'm really lost and totally confused.

$\endgroup$
  • $\begingroup$ Welcome to DSP.SE! This is an interesting problem. A comment: not having the variance of the input will mean you can't know the gain, so you'll have to assume some value for the input. Question: What are you actually going to do with the model once you have it? That may have a bearing on how to approach the problem. $\endgroup$ – Peter K. Apr 25 '16 at 14:00
  • 1
    $\begingroup$ @Peter - Thank you for the warm welcome. Regarding your question: once I've the validated model, then I would like to use the poles (and the corresponding zeros) that stems from that model in a subsequent step. However, it doesn't involve any feedback-loop, so it's all within the context of an 'open loop' system where the outcome (i.e poles of the transfer function) is deduced/taken and substituted in another separate step, which will assume that the information is correct. Therefore, the need of some kind of tangible form of model validation, in order to substantiate this made assumption. $\endgroup$ – WaffleTeX Apr 25 '16 at 15:29
  • $\begingroup$ Just another question: you mention y(t) and Y(k) are they the same? And your data, y(t) is complex? $\endgroup$ – Peter K. Apr 26 '16 at 19:34
1
$\begingroup$

I don't know of a standard way to do this, but I think taking a leaf out of the training of artificial intelligence algorithms might help:

  • divide your output data into two disjoint sets,
  • perform whatever identification procedure you want to obtain the ARMA parameters, and
  • then compare the two transfer functions using whatever distance metric makes sense for your application.

In the figure below, I've decided that my distance metric is that the poles and zeros of each ARMA model are "not too far" from each other (or the ARMA poles and zeros found using all the data).

  • The $\color{black}{\bf black} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from first half of the data ($G_1$).
  • The $\color{red}{\bf red} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from second half of the data ($G_2$).
  • The $\color{grey}{\bf grey} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from all of the data ($G_{\rm tot}$).

enter image description here

Another way to do it might be to look at the transfer function $$ E = G_1 - G_2 $$ That is plotted below using the freqz command in R.

enter image description here


R Code Below

#30332
library(signal)
library(tseries)

# Create a time series that is not just white noise
T <- 1000

do_plot <- 1

#for (k in seq(1,100))
{
  x <- rnorm(T, 0, 1)
  bf <- butter(4,0.41)
  y <- filter(bf$b, bf$a, x)

  # Split data in two and compare ARIMA models
  L2 <- length(y)/2
  y1 <- y[100:L2] # Start away from 1 to avoid startup transients
  y2 <- y[(L2+1):length(y)]

  orders <- c(3,0,3)

  # First half
  tryCatch(
  {arima1 <- arima(y1, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_1 <- coefficients(arima1)
  num_1 <- all_1[4:6]
  den_1 <- all_1[1:3]

  # Second half
  tryCatch(
  {arima2 <- arima(y2, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_2 <- coefficients(arima2)
  num_2 <- all_2[4:6]
  den_2 <- all_2[1:3]

  # All data
  arima_tot <- arima(y, order = orders)
  all_tot <- coefficients(arima_tot)
  num_tot <- all_tot[4:6]
  den_tot <- all_tot[1:3]

  # Plot  boundaries
  xs <- c(-1,1)
  ys <- xs

  if (do_plot == 1)
  {
    par(mfrow=c(1,1), pty='s')
    plot(roots(den_1), xlim=xs, ylim=xs, lwd = 6, pch=4, xlab="Re(z)", ylab="Im(z)")
    do_plot <- 0
  }
  else
  {
    points(roots(den_1), lwd = 6, pch=4)
  }
  points(roots(den_2), lwd = 4, col="red", pch=4)
  points(roots(den_tot), lwd = 2, col="grey", pch=4)

  points(roots(num_1), lwd = 6, pch=1)
  points(roots(num_2), lwd = 4, col="red", pch=1)
  points(roots(num_tot), lwd = 2, col="grey", pch=1)
}


par(pty="m")
subtract_num <- conv(num_1, den_2) - conv(num_2, den_1)
subtract_den <- conv(den_1, den_2)

freqz(subtract_num, subtract_den)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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