VF_spectrum VD_spectrum VE_spectrum
 Function Power-density spectrum
 Syntax C/C++ #include float VF_spectrum( fVector Spc, ui specsiz, fVector X, ui xsiz, fVector Win ); C++ VecObj #include T vector::spectrum( const vector& X, const vector& Win ); Pascal/Delphi uses VFstd; function VF_spectrum( Spc:fVector; specsiz:UIntSize; X:fVector; xsiz:UIntSize; Win:fVector ): Single;
 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_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. Example C/C++ (for DOS): #include #include #include #include #include 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 handling If 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 value PSD at the Nyquist frequency