μblog: engineering from the trenches

When we are faced with the problem of system identification, it is convenient to determine the transfer function, or impulse response, of a system in order to generate a working model. Since at this stage of analysis, the system is almost a black-box, we cannot make too many assumptions about the transfer function beyond what can be determined from the input and output time series. All hope is not lost, however, since there exist formidable methods of transfer function estimation.

The transfer function is formally the Laplace transform of the unit impulse response of a system, or the frequency domain impulse response representation. Where the impulse response is a time series, the transfer function defines the amplitude modulation and phase lag of the system frequency spectra. The classical method to estimate this transfer function is to determine the cross spectral density of the input and output signals and then normalize it to the input spectral density.

To determine the cross spectral density of the two signals, we first convolve the two signals to yield the cross correlation. In the computer, the two time series are either wrapped around or zero padded so that one time series slides across the other time series and the overlap area is measured at a given offset (offset of one signal compared to the other.) If both of the signals have, for example, a 1Hz component, then the cross correlation will also have a 15Hz periodicity. That is to say, if our sampling rate is 10Hz, so each step in the time series represents 100ms, then there will be a peak in the cross correlation every ten data points. The spectral components of the cross correlation can then be examined through the Fourier transform to get the cross spectral density. What this will give, is a measure of which frequency components both signals share, there will be a peak at 1Hz only if both signals contain 1Hz components.

At this point, in an ideal world, we could take the cross spectral density and normalize its magnitude to the input signal spectrum and we would have the transfer function of our system. In a real world situation, we have a problem with system noise that often corrupts our transfer function.

In the real world, step functions are seldom applied for system identification as they are broad band for a brief period of time, instead, noise generators are employed to generate the input signal of the system and the resultant transfer function is processed in the same method as mentioned above. These transfer functions are averaged over long data sets to give a smooth function which is resistant to noise corruption. This trick only works well where the noise is flat (white) in the frequency band of interest.

The next step that can be taken is to try to identify the spectral distribution of the noise and later normalize the cross spectral density to both this noise power and the input signal power. A good example is a system where there is 60Hz noise leaking in and the bandwidth of interest is DC to 100Hz. If we know that there is a peak at 60Hz in the noise spectra, and we see that there is a peak at 60Hz in an otherwise smooth transfer function, then that peak is somewhat suspect and should be investigated to see if it is real or it is due to the noise input of the system. This technique works well when the noise is stable in the frequency domain and is time invariant, it does not work well when the noise frequency distribution moves around in a less-than-deterministic manner, such as in neural recordings.

Based on which state the brain is believed to be in, there are different spectral components to the recorded EEG. Trying to determine the functionality of the brain from these passive recording is a field on its own, however, I am interested in determining the transfer function of neural tissue for stimulation artifact removal. At this point I am still developing a strong method to determine the transfer function in the face of these adversities with some luck in the narrow-band signal domain. By supplying sinusoidal stimulation at known frequencies and averaging over time, I am able to overcome some of these obstacles. The downside is that there are less data-points to use and the transfer function looks to be less-smooth, such as the image at the top.

del.icio.us | digg