Report

Spectrum Estimation W. Rose 2013-04-06 Department of Kinesiology and Applied Physiology Signal x(t) t=0 to T ΔT=sampling interval=1/fsamp N=number of samples T=N ΔT “Simple” spectrum estimate S(f) = power spectrum = k||X(f)||2 X(f)=discrete Fourier transform of x(t) (X(f) is complex) where f = 0 to Δf·N/2 = 0 to fsamp/2 (“single-sided”) f = 0 to Δf·(N-1) (“two-sided”) and where Δf=1/T Department of Kinesiology and Applied Physiology Spectrum Details S(f) = power spectrum = k||X(f)||2 where X(f) = discrete Fourier transform of x(t) X(f) is complex; S(f) is real If x(t) has N points Then two-sided spectrum has N points one-sided spectrum has N/2+1 points k=normalizing factor; depends on particular routine used to calculate spectrum: single or double sided, “peak” or “RMS” units, etc. (1) Department of Kinesiology and Applied Physiology Example of the “simple” spectrum estimate Suppose = 1 sin 2 · · + 3 sin 2 · 3 · + 5 sin 2 · 5 · where f=1 , d1=4/ = 1.273, d3=4/(3) = 0.424, d5=4/(5) = 0.255. This means x(t) is comprised of sinusoidal components with frequencies of 1, 3, and 5 Hz, with amplitudes given by d1, d3, and d5.* Suppose also that the sampling rate is fsamp=100 Hz, i.e. dT=1/fsamp=0.01 s, and suppose that we collect data for 5 seconds (Ttotal=5 s). *dk=4/(kπ)=amplitudes for ±1V square wave (k odd). Department of Kinesiology and Applied Physiology Example of the “simple” spectrum estimate Matlab code to make t and x(t) and plot: dt=0.01 t=(0:499)*dt; f=1; d1=4/pi; d3=4/(3*pi); d5=4/(5*pi); x=d1*sin(2*pi*f*t)… +d3*sin(2*pi*3*f*t)… +d5*sin(2*pi*5*f*t); plot(t,x,’.r-’); Department of Kinesiology and Applied Physiology Example of the “simple” spectrum estimate Matlab code to compute amplitude spectrum using fft(x): X=fft(x); Xamp_2s=abs(X); N=length(x); T=N*dt; df=1/T; freqs=(0:N-1)*df; plot(freqs,Xamp_2s,’.r-’); Note: fft(x) returns a 2-sided spectrum, as is evident from the graph below. Amplitude spectrum (two-sided) Department of Kinesiology and Applied Physiology Example of the “simple” spectrum estimate The non-zero values |X(f)| = dk·N/2, where dk=the amplitude of the sinusoid used to contruct x(t), as shown on an earlier slide. This shows that values returned by fft(x) scale like N/2, except at f=0 and f=fNyquist=fsamp/2, at which frequencies the values returned by fft(x) scale like N (not shown in this example). The “divide by two” scale factor is due to the fact that the energy that started out at 1 Hz is split between 1 Hz and 99 Hz in the two-sided spectrum. Amplitude spectrum (two-sided) Department of Kinesiology and Applied Physiology Example of the “simple” spectrum estimate Compute the single-sided power spectrum. Divide by N at f=0 and f=fNyquist; divide by N/2 at all other frequencies. S=([Xamp_2s(1) Xamp_2s(2:N/2)*2 Xamp_2s(N/2+1)]/N).^2; freq1s=(0:N/2)*df; plot(freq1s,S,'.r-'); xlabel('Frequency (Hz)'); ylabel('Power'); This gives S(f)=d12,d32,d52 at the appropriate frequencies. Single-sided power spectrum Department of Kinesiology and Applied Physiology A Better Spectrum Estimate Overview 1. Divide the time domain record into blocks. 2. Find power spectrum of each block. 3. Average the power spectra, frequency by frequency. Department of Kinesiology and Applied Physiology A Better Spectrum Estimate: Details 1. Divide time domain record into segments of equal length. If there are q non-overlapping segments, then also include q-1 half-overlapped segments. Example: Total data record length=Ntot=4000 points. Investigator chooses q=4. Then each segment has length Nseg=Ntot/q=1000 points. The four nonoverlapping blocks start at points 0, 1000, 2000, 3000. Three half-overlapped blocks start at 500, 1500, 2500. Total number of segments=2q-1=7. Department of Kinesiology and Applied Physiology A Better Spectrum Estimate: Details 2. Find power spectrum of each segment. Before computing spectrum of each segment, remove linear trend, resulting in block with zero mean value and zero slope. Some prefer to remove only the mean value and not the linear trend (if any). Window the data in the segment with Hann or Hamming window. Compute power spectrum of windowed data. Correct power spectrum for the loss of power caused by the window: Hann window: multiply power estimates by 8/3. Hamming: multiply power estimates by 2.516. Department of Kinesiology and Applied Physiology A Better Spectrum Estimate: Details 3. Average the power spectra, frequency by frequency. Savg(f)=(1/(2q-1))*Sum(i=1 to 2q-1){Si(f)} where Si(f)=power at frequency f of the ith segment Department of Kinesiology and Applied Physiology Department of Kinesiology and Applied Physiology Scale Factors for Power Spectrum Estimates “single-sided” spectrum: At each non-zero frequency (from Δf to highest frequency below the “halfway point”, which is f=fsamp/2), multiply the two-sided estimate at that frequency by 2 to get the single-sided power estimate. The single sided spectrum only goes to fsamp/2. At f=0 and at f=fsamp/2, the single sided estimate is the same as the two-sided estimate. See Labview help for Power Spectrum.vi (returns two-sided spectrum) and Auto Power Spectrum.vi. (returns one-sided spectrum). Department of Kinesiology and Applied Physiology “RMS” units Some routines, such as Labview’s AutoPowerSpectrum.vi, give S(f) in “root mean square” units by default. The net effect is to divide the power spectrum estimate by 2, if f>0, compared to the power spectrum estimate given in “peak” units. = 0 = > 0 2 Reason: The average power of a signal is defined as 1 2 = 0 This definition was chosen because it equals the power dissipated in a 1 ohm resistor by a current I=x(t). If x(t)=Acos(2πft), and if T=a whole number of cycles, then we can solve the integral to show that P=A2/2 if f>0, and P=A2 if f=0. It is called “root mean square”, because a time-varying signal has much power as a steady signal whose amplitude equals the square root of the mean of the squared amplitude of the time varying signal. Department of Kinesiology and Applied Physiology Partial front panel of FFT_and_Power_Spectrum_Units.vi Example VI in LV2012 Output from FFT.vi scales like record length n. Other VIs return output whichh is independent of record length. FFT and Power_Spectrum return two-sided spectra, which is why their non DC scaling is divided by 2. The other VIs return one-sided spectra. FFT_Power_Spectral_Density divides each power spectrum estimate by Δf, the frequency spacing between estimates (Δf =spectral width of each estimate). Department of Kinesiology and Applied Physiology