| Description | The data set X is analyzed for its power spectral density (PSD), i.e. the mean square amplitude. The result is stored in Spc. xsiz must be at least 2*specsiz, and specsiz has to be an integer power of 2. Internally, X is divided into xsiz*specsiz/2 segments and the average over the spectra of the individual segments is calculated. Each segment of length 2*specsiz yields the PSD for specsiz+1 frequencies (see VF_FFT). In order to keep specsiz an integer power of 2, there are only specsiz points stored in Spc and the last one, the PSD at the Nyquist frequency fNyquist = 0.5 / sampling_interval, is given as the return value of the function. It may either be neglected (by calling the function like a void function) or stored as the last element in Spc by calling the function as
Spc[specsiz] = VF_spectrum( Spc, specsiz, X, xsiz, Win );
in this case, Spc must have a length of specsiz+1.
Win is a window that is applied to the data segments. The size of the Win vector must be 2*specsiz. Within the VectorLib library, three functions are available that give suitable Windows: VF_Welch, VF_Parzen, and VF_Hanning. A square window (i.e. no windowing at all) is achieved by setting all elements of Win to 1.0 using VF_equ1. Use of the square window is not recommended here, though.
You may wish to test the quality of the calculated spectrum by applying Parseval's theorem (provided you called VF_spectrum as in the above example and stored the PSD for the Nyquist frequency):
1.0/xsize * VF_ssq( X, xsize ) must be about equal to VF_sum( Spc, specsiz+1 ).
If the deviation between both results is large, the sampling interval in X probably is too large.
About special versions with the prefixes VFp_, VFs_ and VFl_, consult chapter 4.8. |
| Example C/C++ (for DOS): | #include <VFstd.h>
#include <VFmath.h>
#include <Vgraph.h>
#include <math.h>
#include <conio.h>
void main( void )
{
fVector X, Spc, Win, Time, Freq;
float deltat, fOscill, fNyquist;
ui xsize=4096, specsiz=256;
X = VF_vector( xsize );
Spc = VF_vector( specsiz+1 );
Win = VF_vector( 2*specsiz );
Time = VF_vector( xsize );
Freq = VF_vector( specsiz+1);
deltat = 0.001; /* sampling interval 1 millisecond */
fNyquist = 0.5 / deltat; /* Nyquist frequency = 500 Hz */
fOscill = 100; /* say we are sampling a 100 Hz oscillation */
VF_ramp( Time, xsize, 0, deltat );
VF_ramp( Freq, specsiz+1, 0, fNyquist / specsiz );
VFx_sin( X, Time, xsize, 2*M_PI*fOscill, 0, 1 );
/* sine wave (omega*t) */
VF_cmpC( X, X, xsize, 0.7 );
/* convert into asymmetric square wave */
VF_Welch( Win, 2*specsiz ); /* or another window */
Spc[specsiz] = VF_spectrum( Spc, specsiz, X, xsize, Win );
V_initGraph( "\\BorlandC\\BGI\\" );
/* give the correct BGI path! */
VF_xyAutoPlot( Freq, Spc, specsiz+1,
PS_SOLID | SY_CROSS, GREEN );
getch(); /* hit any key when you have seen enough */
closegraph();
V_nfree( 5, Freq, Time, Win, Spc, X );
} |
| Example Pascal/Delphi: |
uses VFstd, VFmath, Vgraph,...;
const
xsize=4096; specsiz=256;
var
X, Spc, Win, Time, Freq: fVector;
deltat, fOscill, fNyquist: Single;
begin
X := VF_vector( xsize );
Spc := VF_vector( specsiz+1 );
Win := VF_vector( 2*specsiz );
Time := VF_vector( xsize );
Freq := VF_vector( specsiz+1);
deltat := 0.001; (* sampling interval 1 millisecond *)
fNyquist := 0.5 / deltat; (* Nyquist frequency = 500 Hz *)
fOscill := 100; (* say we are sampling a 100 Hz oscillation *)
VF_ramp( Time, xsize, 0, deltat );
VF_ramp( Freq, specsiz+1, 0, fNyquist / specsiz );
VFx_sin( X, Time, xsize, 2*PI*fOscill, 0, 1 );
(* sine wave (omega*t) *)
VF_cmpC( X, X, xsize, 0.7 );
(* convert into an asymmetric square wave *)
VF_Welch( Win, 2*specsiz ); (* or another window *)
VF_Pelement( Spc,specsiz )^ :=
VF_spectrum( Spc, specsiz, X, xsize, Win );
VF_xyAutoPlot( Freq, Spc, specsiz+1,
PS_SOLID + SY_CROSS, GREEN );
V_free( Freq ); V_free( Time ); V_free( Win );
V_free( Spc ); V_free( X );
end;
(* It is assumed that the graphics routines have already been initialized (DOS), or that the above example stands within a Paint (Pascal for Windows) or a ShowView (Delphi) routine. *) |
|