Send q VRf hpvxd_xc@hpislsup.lvld.hp.com

from: rondes@rostra.com

to: VRf hpvxd_xc@hpislsup.lvld.hp.com

date: Tuesday, 29 April 1997

Hello,

Other than the Spectrum object, how do I get a numeric frequency reading from

a waveform? I need to obtain the number so I can save it with other data to

a file.

Thanks for your time and attention,

Ron DesGroseilliers Jr.

rondes@rostra.com

Rostra Precision Controls

910-276-4853

from: rondes@rostra.com

to: VRf hpvxd_xc@hpislsup.lvld.hp.com

date: Tuesday, 29 April 1997

Hello,

Other than the Spectrum object, how do I get a numeric frequency reading from

a waveform? I need to obtain the number so I can save it with other data to

a file.

Thanks for your time and attention,

Ron DesGroseilliers Jr.

rondes@rostra.com

Rostra Precision Controls

910-276-4853

gvg@lvld.hp.com / 970-679-3030 / FAX 970-679-5971

to: Ron DesGroseillers JR

date: Tuesday, 29 April 1997 0929 MDT

> from: rondes@rostra.com

> to: VRf hpvxd_xc@hpislsup.lvld.hp.com

> date: Tuesday, 29 April 1997

>

> Hello,

>

> Other than the Spectrum object, how do I get a numeric frequency reading from

> a waveform? I need to obtain the number so I can save it with other data to

> a file.

>

> Thanks for your time and attention,

>

> Ron DesGroseilliers Jr.

> rondes@rostra.com

> Rostra Precision Controls

> 910-276-4853

Sir:

This is more complicated than you might think ... see below:

--------------------------------- cut here ----------------------------------

[16.8] SIGNAL PROCESSING, FFTS, & VEE

* An Australian customer contacted us with a problem: he was sampling a

sine-type signal with a data-acquisition box and wanted to know how to use

VEE to determine the signal amplitude, DC offset, phase, and frequency from

the data returned.

The first temptation was to use a few simple Formula boxes to do the job.

Given an array of data named "D", the amplitude could be given by:

(max(D) - min(D))/2

-- and the DC offset (average value, really) could be given by:

(max(D) + min(D))/2

However, figuring out the frequency leads to a nightmare consideration of

zero-crossing detectors, and figuring out the phase is even worse. Such

techniques are also very sensitive to noise.

This was clearly a job for VEE's Fast Fourier Transform (FFT) function ...

but the details on exactly how to use the FFT to get the job done are tricky.

This article explains the issues involved.

* The first problem is understanding exactly what an FFT is and what it does.

Most engineers have a simple understanding of the Fourier Transform in

general: that it converts a signal represented as, say, a voltage over a

period of time (a "time domain" signal) into a "spectrum" of sine functions

of different amplitude and phase (a "frequency domain" signal).

This is perfectly accurate as far as it goes, but when you get into dealing

with "discrete" data samples on a computer (instead of simply manipulating

math equations on the chalkboard) then the details get a little murky.

In this case, you must use a Discrete (or "Digital", if you like) Fourier

Transform (DFT) ... the elementary DFT algorithm is slow, so in practice we

generally use a more complicated but much faster FFT, but as far as the user

is concerned (with one catch that will be discussed at the end of this

article) they are exactly the same thing, and so this discussion will refer

to the DFT.

A DFT has some interesting properties. Let's consider a DFT on a sampled

time domain data set of, say, 256 points ... though the following comments

apply to any other number of points just as well:

% If you perform a DFT on this set of 256 points, you get back the same

number of frequency-domain values -- 256. The values returned by the DFT

are in the form of a set of complex numbers, in order to represent both

the phase and the magnitude (or amplitude, if you prefer) of sine curves

that make up the spectrum ... note that the math works out much more

easily if these values are indexed starting from 0 (making the last complex

value 255).

% The frequencies of the sine curves in the spectrum returned by the DFT are

determined by the rate at which the original time-domain data was sampled

... *not* the frequency of the original signal, mind you -- the

frequency at which it was *sampled*. If your sampling rate was 50 KHz,

then the highest frequency represented in the spectrum (by the complex

value with index 255) is 50 KHz.

All the lower values have periods one sample interval longer than this

last value -- that is, the element with index 254 has a period of two

sample intervals, that with index 253 has a period of three sample

intervals, all the way down to the element with index 1, which has period

255 times the sample interval.

In reverse terms, the highest frequency is the highest "harmonic" --

integer multiple -- of the sine curve represented by the value in the

spectrum with index 1 ... that means this "base harmonic" has a frequency

1/255th of the highest harmonic.

Note, for curiosity, that the only harmonic values in the spectrum are

those that have periods some multiple of the sample rate. This might seem

odd, but it is inherent in the fact that you are dealing with a sampled

signal: the sampling process only obtained values of the input signal at

the sample intervals, and so the reflection of the input signal into the

frequency domain is governed by the sample rate as well.

Note also that this discussion "talked around" the complex value with

index 0. This is the "zeroth harmonic" ... which, in plain language, is

just the "DC component" of the signal. It will always have an imaginary

part of 0.

% The magnitude of the complex values returned by the DFT algorithm

provides information that is a little hard to explain. What it actually

gives back is the sum of all the samples of data that represent a given

harmonic.

For example, suppose we have our 256 points and they are sampling exactly

four cycles of a sine wave with a peak amplitude of 3. This means that

after we take the DFT of this data, we end up with a set of complex values

that are zero (or effectively zero) everywhere but in element 4 ... which

has a magnitude of 768! This seems crazy, but the input has 256 samples

of a harmonic with amplitude 3, so 3 * 256 = 768 is the right value.

This is result is actually the way the DFT is defined, and for

applications like image processing it's what you want ... but people who

are used to dealing with spectrum analyzers and the like get rather

upset with it -- the sine signal has a value of 3, by Bob, and they expect

to get a value of 3 out of the DFT!

This can actually be obtained in practice by "scaling" the output of the

DFT ... all you have to do is divide it by the number of original data

points, which is 256 in this case.

% Now, a useful digression: there's no information in the DFT itself that

tells you what any of the frequencies really are ... the complex values in

the spectrum are really only defined as multiples of the base harmonic and

do not have any specific frequency assigned to them. This makes sense

because the sampled data set could represent anything from picoseconds to

centuries -- the spectrum that results would look exactly the same. You

have to keep track of the actual frequencies yourself and do some

additional computations to factor them in to the results.

% OK, another twist: according to a little rule of signal processing called

the Nyquist criterion, a sampling process can only detect sine harmonics

that are no more than half the sampling rate ... think about it: a sine

curve varies from positive to negative -- and if you sample one over an

interval any longer than half a cycle, you can end up getting two positive

samples and completely miss the negative half of the cycle (or the

reverse). This critical frequency is known as the "Nyquist frequency".

Attempting to sample a signal that is over half the sample rate can result

in bogus harmonic values known as "aliases", and you can't compensate for

them in software -- you must have hardware filtering on your

data-acquisition system to ensure that the input signal doesn't

incorporate any frequencies higher than this.

This leads to an awkward question: if the highest harmonic in the complex

value set returned by the DFT has the same frequency as the sampling

frequency -- but you can't really detect any frequencies that are any

higher than half that frequency ... then what sense does this highest

harmonic make? In fact, what sense do any of the harmonics that

correspond to sampling any faster than half the sampling rate make?

The answer is that, for purely real time-domain signals (which is all we

will be dealing with here, avoiding the nasty question of what a "complex

signal" is), the spectrum is split in two: the first half consists of

the DC term and all the harmonics to that just below the Nyquist

frequency; the upper half begins with the harmonic with a frequency

corresponding to the Nyquist frequency, while all the following complex

values "mirror" those of the lower half ... ... they have the same

magnitudes as the lower half, though their phase is the negative of their

counterparts -- that is, they are "complex conjugates" of the lower half

of the spectrum.

The last value in the set corresponds to the complex conjugate of the

first harmonic; the value before that corresponds to the complex conjugate

of the second harmonic; and so on.

The fact that the DFT algorithm generates two sets of values for each

harmonic component means that it splits the signal power in half; so if

you are scale the complex values, not only do you divide it all by the

number of samples ... but you only take the first half of the harmonics,

and multiply them by two.

This means that if you have a DFT of a 256-point time-domain signal, then

what you actually get back is a DC component (element 0 of the complex

value set) plus 128 harmonic values (elements 1:128 of the complex

value set) ... you can throw away the other 127 values.

In fact, this is exactly what VEE does for its FFT -- it only returns 129

complex values for a 256-value time-domain signal -- a DC value and 128

harmonic values. To scale them, you would have to divide them by the

original number of time-domain samples -- 256 -- and then, because you

only have half the harmonics, multiply all the complex values (except the

DC value) by 2.

Having waded through all this headache, let me summarize briefly:

% If you put 256 time-domain signals into a DFT, you get 256 complex values

back out.

% The complex values consist of a DC component (the first element in the

set) and a series of harmonics of increasing frequency -- starting at the

base frequency, then the second harmonic, then the third, and working up

to a harmonic that matches the sampling frequency (the last element in the

set). The periods of all the harmonics are integer multiples of the

sampling interval.

% To get actual signal magnitudes from a DFT, you have to divide its output

by the original number of time-domain samples.

% The DFT does not have any frequency information in it ... all it gets is a

set of sampled data and does not factor in the sampling interval. If you

want actual frequencies out of DFT complex values, you have to figure that

in yourself.

% The only frequencies valid in the output of a DFT are those less than half

of the sampling frequency (due to the Nyquist criterion). For real

time-domain data, the upper half (except the first element) of the output

of the DFT contains a redundant set of complex-conjugate "aliases" of the

lower half. This means that the signal magnitudes have been cut in half,

so you not only need to divide all the data by the number of samples to

scale it ... you also must multiply all the harmonics by 2.

* So ... to extract information from an input signal, you must first

scale it ... you can easily do this with the following UserFunction:

+----------------------------------------------------------------------+

| UserFunction: ScaleFFT |

+----+----------------------------------------------------------+------+

| | | |

| | +--------------------------+ | |

| | | Formula | | |

| | +----+---------------------+ | |

| | +-->| Td | fft(Td)/totsize(Td) +--+ | |

| | | +----+---------------------+ | | |

| | | | | |

| Td +--+ +---------------------------+ +-->| CpxD |

| | | | | |

| | | +--------------------------------------+ | | |

| | | | Formula | | | |

| | | +---+----------------------------------+ | | |

| | +-->| A | concat(A[0],2*A[1:totsize(A)-1]) +--+ | |

| | +---+----------------------------------+ | |

| | | |

+----+----------------------------------------------------------+------+

The first Formula box divides the FFT of the input data by the number of

points in the original signal; the second multiplies all the harmonic

components by 2 to generate the proper result.

* OK, given this, what sort of data can be extracted with it? For the

problem we're supposed to address, we're trying to sort out a sine wave from

sampled data and determine its amplitude, DC offset, phase, and frequency.

The following VEE objects can extract this data from the complex data output

by the ScaleFFT() function above;

+-------------+

---+-->| ScaleFFT(x) +--+ +-------------+

| +-------------+ | | Formula |

| | +---+---------+

| +----------------+------------->| C | re(C[0] +--> DC offset

| | +---+---------+

| |

| | +---------------------------------+

| | | Formula |

| | +--------+ +---------------------------------+

| +-->| mag(x) +--+-->| 1 + maxIndex(X[1:totsize(A)-1]) +--+

| | +--------+ | +---------------------------------+ |

| | | |

| | | +--------------------------------+

| | | |

| | | | +----------+

| | | | | Formula |

| | | | +---+------+

| | +-------|--------->| M | M[I] +--> Amplitude

| | +--------->| I | |

| | | +---+------+

| | |

| | | +-----------------+

| | | | Formula |

| | | +---+-------------+

| +-----------------------|-->| C | phase(C[I]) +--> Phase

| +-->| I | |

| | +---+-------------+

| |

| | +-----------------+

| | | Formula |

| +------------+ | +-----------------+

+-->| totsize(x) +----------|-->| S | R * (I/S) +--> Frequency

+------------+ +-->| I | |

+----->| R | |

| +---+-------------+

sample frequency ---+

Some notes to explain how this stuff works:

% You obtain the DC offset just by taking the real component of the first

value in the complex value output from the ScaleFFT() function.

% You obtain the amplitude of the main harmonic by taking the magnitude

of the ScaleFFT() output data and searching for the index at which the

value is maximum (while avoiding the DC value in element 0); once this

index is obtained, you can then get the amplitude for that index.

% You want to get the index because it allows you to determine phase and

frequency as well. To get the phase, you simply get the element of

the complex value output for that max-value index and take its phase.

% To obtain the signal frequency, you take the ratio of the maximum-value

index to the total number of input sample points and multiply it times

the sample frequency (which has to be specified previously). Please

make sure you set the data type on the input pins to Real, since if you

get integer inputs you may end up with an integer divide (and always

mysteriously get "0" as a result).

Experiments with these tools show them to return highly accurate results --

as long as you have a set of samples that match an integer number of cycles.

If the input contains a fractional number of cycles, it creates an apparent

discontinuity in the signal that adds high-frequency components to the FFT

spectrum, reducing the validity of the conclusions.

The DC offset value tends to remain stable under such variations, as does the

frequency -- but the signal amplitude varies considerable, and the phase

varies widely. Further investigation did not resolve these problems and I

can only suggest that the input waveform be trimmed to an integer number of

cycles (using displays with marker inputs) to get the best results.

In principle, use of VEE windowing functions to filter the input data can

help reduce the effects of such discontinuities, but my experiments with them

showed them to be ineffectual in this case. I also discovered that if you

want to get scaling to come out right, you have to multiply the output of the

windowing function by 2 ... this didn't seem right, but I didn't have the

background to investigate further.

* While I said earlier that an FFT was really just a version of the DFT that

was much faster, it does have the peculiarity that it wants to perform

computations on data sets whose number of elements is a integer power of two

(2, 4, 8 , 16, 32, 64 ... ). VEE will use an FFT algorithm if the input data

is this size; if it isn't, it really uses a slower DFT algorithm.

For N values, the time required to perform a DFT is proportional to N^2; for

an FFT, it is proportional to N*LOG(N). For large values of N, the

difference is one of days versus minutes, so you want to make sure your data

is a power of two. All you have to do is add zero elements to the end of

your array ... it doesn't really affect the calculation.

The following function, PadZero(), does this for you:

+--------------------------------------------------------------------+

| UserFunction: PadZero |

+-+----------------------------------------------------------------+-+

| | | |

| | +----------------------------------------+ | |

| | | Formula | | |

| | +---+------------------------------------+ | |

| | +-->| D | intPart(log(totsize(D)+0.1)/log(2) +--+ | |

| | | +---+------------------------------------+ | | |

| | | | | |

| | | +--------------------------------------------+ | |

| | | | | |

| | | | +--------------------------------------------------+ | |

| | | | | Formula | | |

| | | | +---+----------------------------------------------+ | |

| | | +-->| E |(S=2^E ? A : concat(D,ramp((2^(E+1)-S),0,0))) | | |

|D+--+----->| D | +-->|R|

| | | +-->| S | | | |

| | | | +---+----------------------------------------------+ | |

| | | | | |

| | | +---------------+ | |

| | | | | |

| | | +-----------+ | | |

| | +-->| totsize() +--+ | |

| | +-----------+ | |

| | | |

| | | |

+-+----------------------------------------------------------------+-+

The formula at top takes the number of elements of the input data, determines

the power of two that gives that number, and then truncates it to an integer

(note that I added a fudge-factor of 0.1 to ensure that roundoff error will

not cause problems).

The formula on the bottom then determines if the number of elements in the

input data is exactly the same as 2 raised to the integer power defined in

the first formula. If so, it just passes on the array; if not, it

concatenates the array with an array of zeroes long enough to make the number

of elements the value of 2 raised to the next highest integer power. (Note

how the ramp() function is used to generate this array.)

* A test function (xfftfunc.vee) showing off ScaleFFT() and PadZero() is

available.

[%%]