## Archive for the ‘**DSP**’ Category

## Draft Paper for Citizen Scientist

## Auto-Correlation of Rectified Signal re DEMON spectral Analysis

(Slightly modified version to be published on Citizen Scientist in Feb of March 2010)

It has been pointed out (reference 1) that the speed of a domino wave can be estimated from peaks in the auto-correlation of the rectified recording of the noise of the domino array collapse (note here I am ignoring filtering etc which is useful with the auto-correlation method as it is also with the DEMON method (references 2 and 3). There can be no argument with this in principle since for periodic phenomenon the auto-correlation and spectrum of the rectified signal both encode the frequency and/or period of the phenomenon. The auto-correlation will contain a series of peaks at intervals corresponding to the period of the phenomenon (and the time delay of the for such peak will be at a value corresponding to the period), while the DEMON spectrum will have a peak at a frequency corresponding to the frequency (1 over the period) of the phenomenon.

### Time Domain Auto-Correlation Complexity

One option for computing the auto-correlation of a data vector is to do the calculations in the time domain. The pseudo-code below shows how this could be done.

##### ————————————————————————————————–

##### Psuedo-Code for Time Domain Auto-Correlation:

// Assume input data (real) is in array(1..Npoints) and // output data is in array ss(1..Ndelays) with the first // point corresponding to zero delay ss=zeros(1,Ndelays) for delay=0 to Ndelays-1 for idx=delay+1 to Npoints ss(delay+1)=ss(delay+1)+data(idx)*data(idx+delay); endfor endfor -----------------------------------------------------------

This allows us to compute the number of floating point add and multiply operations needed for this time domain auto-correlation, which is:

Now if as is often the case , then we can write:

Note: loop control variable increments and array index arithmetic operations are not counted as these are generally much faster than the (usually) double precision floating point operations.

### Frequency Domain Auto-Correlation Complexity

The main alternative to the above time domain method for computing the auto-correlation is to use the convolution theorem to convert the auto correlation calculation from what is a discrete convolution in the time domain to a point wise product in the frequency domain, and then transforming back into the time domain. To do the conversion from time to frequency domain we will be using the Fast Fourier Transform (FFT) (see reference 4 for some general background for the FFT), which is an efficient implementation of the Discrete Fourier Transform (DFT). This is the Fourier transform for a discrete periodic signal, so to avoid wrap around effects we need to pad the data with sufficient zeros so the wrap around does not effect the part of the result we are interested in (that is the FD auto-correlation we will compute using the DFT will be a circular auto-correlation). This means that if we want points of the auto-correlation we will need to pad our data vector with zeros giving an effective data vector length of .

Now the equivalent of the convolution theorem for correlations tells us that:

Where denotes the (circular) correlation operator (not to be confused with the convolution operator ), the inverse DFT operator, the DFT operator, “.” the pointwise product operator, is the correlation delay which starts at and index arithmetic is modulo the length of . For convenience it is best to regard the indicies to follow the C/C++ convention and start at and run up to , where is the length of rather than form to .

That this works can be demonstrated on zero mean random complex simulated data ( points sample delays), the absolute value of the difference between the time and frequency domain results are shown in figure 1 below, where we can see that the differences are of the order of (all calculations done in double precision), which is pretty good agreement.

where is the length of the transform, and is a constant of the order of . for a complex transform and half that for a real transform.

So if we wish to compute points of the auto-correlation of a signal of length we will need two FFT’s of length (one forward transform and one inverse transform) and an additional complex multiplies, which is equivalent to three times as many real multiplies and one real add but half of these can be eliminated if the original signal is real because of the symmetries in the FFT of a real signal.

So we have (to a handwaving approximation for the coefficient and a conservative allowance for the number of multiples) for a frequency domain auto correlation the floating point operations count is:

### Comparison Of Operations Counts for Time and Frequency Domain Auto-Correlation

We may plot the operations counts corresponding to equations 1 and 2, when we get the plot in figure 2, where we see a crossover where the frequency domain auto-correlation becomes more efficient than the time domain algorithm somewhere in the region of 50 delays (this is for a data sample of 30000 points, a value typical of a recording of a domino wave of duration slightly over 1 second). As for domino wave we know that we are looking for impact frequencies down to something between 15 and 20 Hz for a sampling frequency of the order of 20 kHz we will need to look at auto-correlation delays of up to something over 1000 sample intervals. In which case the frequency domain auto-correlation is significantly more efficient.

Also of note in Figure 3 is the variation in the timing of the FD Auto-Correlation, this is because a mixed radix FFT is in use and its efficiency varies somewhat with the length of the FFT.

### Complexity of DEMON analysis

In the previous section we have seen that if we wish to analyse a signal using the auto-correlation of the rectified signal we will use a frequency domain implementation of the auto-correlation. This is more or less equivalent in operations count to two FFT’s of length equal to the signal length plus the number of points in the auto correlation we require. However when analysing a recording of a domino collapse wave we can just as well look for the frequency of the peak in the DEMON spectrum (the spectrum of the absolute value of the signal) (references 2 and 3). This involves just one FFT so it requires no more than half the operations that looking at the auto-correlation requires.

This difference is largely academic if we are doing a one-off analysis but will become significant when a large number of analyses are to be performed (and in an embeded system analysing signals in real time even more so).

### Refences

1. Hannon J., “Domino Speed”, Citizen Scientist, May 2009, http://www.sas.org/tcs/weeklyIssues_2009/2009-05-01/project1/index.html

2. Larham R., Measuring the Speed of the Domino Effect Using the Windows Sound Recorder, Part 1, Citizen Scientist, Nov 2007, http://www.sas.org/tcs/weeklyIssues_2007/2007-11-02/project2/index.html

3. Larham R., Measuring the Speed of the Domino Effect Using the Windows Sound Recorder, Part 2, Citizen Scientist, Dec 2007, http://www.sas.org/tcs/weeklyIssues_2007/2007-12-07/project2/index.html

4. Press W. H, Teukolsky S. A. , et al. “Numerical Recipes in C” (or any other of the series), CUP, 2-nd ed., 1992.

## More Domino Wave Speed Stuff

Discussions in Backscatter Column of Citizen Scientist (ref 4) of my paper and experiments on measuring the speed of the Domino Effect (ref 1 and 2) and attempts to repeat the experiments also reported in CS seem to miss the point of using DEMON spectra (ref 3).

I may be partially responsible for this since I have in the past commented that the spikes of noise which I presume to be the clicks of domino impacts on closer examination appear to be tonals modulated by a narrow envelope. I now believe that this was partially a misinterpretation of the data and partially a result of over fierce band pass filtering on my part.

Whatever the truth of the above we can observe directly from a plot (fig 1) of the auto-correlation of a recording of the domino wave, that there is no clear signal at the expected delay (of about 0.04 s for this data file):

Fig 2 shows the DEMON spectrum for this data file and here we see a “Sore Thumb” signature at about 24.5 Hz, which would correspond to a delay of about 0.04 seconds in the auto-correlation for this data set, where there is no clear feature. (the 24.5 Hz frequency corresponds to a wave speed close to 1 m/s on an 0.04 m pitch domino array, or a normalised wave speed of ~1.4 for dominoes of height ~0.052 m and thickness ~-.008 m (the width does not affect the wave speed but here its is ~0.025m)

Note: The datafile d9p2.wav is the recording of one of the repeat runs, not reported in refs 1 and 2, that I have done to look at the variability in speed measurements. In this case the speed measured is within about 2% of that previously reported.

References:

1. Larham R., *Measuring the Speed of the Domino Effect Using the Windows Sound Recorder, Part 1*, Citizen Scientist, Nov 2007, http://www.sas.org/tcs/weeklyIssues_2007/2007-11-02/project2/index.html

2. Larham R., *Measuring the Speed of the Domino Effect Using the Windows Sound Recorder, Part 2*, Citizen Scientist, Dec 2007, http://www.sas.org/tcs/weeklyIssues_2007/2007-12-07/project2/index.html

3. Smith R., *Attempted Recreation of Ron Larham’s Domino Wave Experiment*, Citizen Scientist, March 2009, http://www.sas.org/tcs/weeklyIssues_2009/2009-03-06/project1/index.html

4. Hannon J., *Still More About Domino Waves*, Backscatter Column, Citizen Scientist, April 2009, http://www.sas.org/tcs/weeklyIssues_2009/2009-04-03/backscatter/index.html

## Blocking Matrices

**Proof that the Blocking Matrices do Indeed Block Signals from Specified Directions.**

For some/certain signal processing tasks we wish to generate beam patterns with nulls in specified directions (or equivalently; filters with narrow notches at specified frequencies). Typically these can be used to null out stationary interference sources, directions to multipath images of a desired source, …

We assume that we have a number of sensors ( )which may be array elements, subarrays or formed beams and we wish to null signals from directions , and , and responses for a source in direction (being a column vector of complex signals one element for each sensor). In a narrowband system we may assume there are simply complex numbers representing the sensor amplitude and phase response. We let denote the normalised versions of these signals.

It is obvious (when pointed out anyway) that the matrix:

is a blocking matrix for directions (where is the matrix with for its columns). The space of possible sensor outputs constitutes a (complex) vector space of dimension and the subspace is the null space of and behaves like the identity transformation on the orthogonal complement of the null space.

To see that this is a blocking matrix for we need just consider:

so:

Now the columns of are the vectors and so these are zeros vectors as expected.

To see that behaves like the identity on the orthogonal complement consider in the orthogonal complement of the space spanned by the columns of . Then since each component of this is the dot product of a column of and and so zero. Hence:

(I hope that is all clear, WordPress LaTeX seems to have started deleting all the \ characters from the LaTeX strings again so there may still be some omissions from when I have tried to restore them)