Real spectrum analysis with Octave and MATLAB
I use Fourier Transforms on time domain data all the time now, whether it's measuring OPAMP noise or looking for spurious signals with my digital oscilloscope acting as the digitizer. Being able to get a calibrated spectrum display is very useful when verifying and troubleshooting nearly any design.
Most modern oscilloscopes now have a DFT/FFT1 display mode built in and that's fine, but you are stuck using the built-in definitions and DFT implementation and I have yet to see one that will handle noise measurements properly. Also, you may just have data from some other source that you would like to quickly analyze.
While there are many examples on the Web of using DFT functions that show how to make a spectrum for signal analysis or power spectral density for noise analysis, many are not fully accurate, hard to use, or even worse, not compatible with each other. All of this confusion leads to a lot of wasted time. These problems lead me to create a set of consistent analysis functions for Octave/MATLAB2 that I use all the time for calibrated DFT signal and noise analysis in a unified, easy to remember, and consistent form.
With Octave, data from any source may be analyzed and generally Octave is quicker than traditional programming languages when a solution needs to be developed fast.
Every DFT analysis starts with an input signal
When testing DFTs a method of quickly generating test signals is needed. There are literally millions of function snippets on the Web for generating signals but they all use a different approach and that makes them hard to remember and reuse. I have found that I need to generate deterministic signals in two forms:
- Generate a deterministic signal with a known number of sine wave cycles. That way I can test windowing functions with whole and fractional cycles to easily determine amplitude accuracy, etc.
- Generate a deterministic signal based on a sampling approach. This form is useful when I want to make a 10,000 point vector of a 10 kHz sine wave sample at a sampling rate of 100 kHz. This is useful when thinking in terms of an actual digitized signal.
Then for white noise signals:
To generate Gaussian white noise with a given Power Spectral Density (PSD) in real units that we analog folks use all the time, like: Vrms/√Hz.
Note that in all these examples, and for the rest of this discussion, the units for the input signal are normally assumed to be in Vrms. If the data is in different units then it should be scaled appropriately back into Vrms units before continuing the signal processing.
Starting with the generation of sine wave signals. For scenario #1 above the function: “tonecycles.m” is used as shown in Figure 1.
% Make an input signal (Sine Wave Tone) based on number of cycles
cycles = 3.95; % Number of cycles to make
points = 1000; % Number of points to generate
amp = 1; % Amplitude Vrms
time_data = tonecycles(amp, cycles, points);
If a sinewave based on sampling is needed, use the tonesampling() function as shown in Figure 2.
% Make an input signal (Sine Wave Tone) based on sampling
fs = 100000; % Sampling Rate Hz
fsig = 500; % Tone Frequency Hz
points = 1000; % Points
amp = 1; % Amplitude Vrms
time_data = tonesampling(amp, fsig, fs, points);
Figure 2 To make a sine wave signal based on a sampling approach, the function “tonesampling.m” is used. The resulting output plot is shown (bottom).
See the example file: “generate_sinewaves_example.m” on the last page of the article for a complete Octave example of Figures 1 and 2 with plots.
To make white noise of a specified power spectral density, the function: “noisepsd.m” is used as shown in Figure 3. A complete working Octave example of the noisepsd() function is provided in the file: “generate_noise_example.m”.
% Noise PSD Setup (Random Noise)
fs = 10000; % Sampling Rate Hz
points = 10000; % Points
amp = 1.0; % Amplitude Vrms / rt-Hz
% Calculate the noise
noise_data = noisepsd(amp, fs, points);