VF_spectrumVD_spectrumVE_spectrum
FunctionPower-density spectrum
Syntax C/C++#include <VFstd.h>
float VF_spectrum( fVector Spc, ui specsiz, fVector X, ui xsiz, fVector Win );
C++ VecObj#include <OptiVec.h>
T vector<T>::spectrum( const vector<T>& X, const vector<T>& Win );
Pascal/Delphiuses VFstd;
function VF_spectrum( Spc:fVector; specsiz:UIntSize; X:fVector; xsiz:UIntSize; Win:fVector ): Single;
DescriptionThe 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_Hann. 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. *)

Error handlingIf size is not a power of 2, VF_FFT (on which VF_spectrum is based) complains "Size must be an integer power of 2" and the program is aborted.
Return valuePSD at the Nyquist frequency
See alsoVF_FFT,   VF_convolve,   VF_autocorr,   VF_xcorr,   VF_filter

VectorLib Table of Contents  OptiVec home