1 / 16

Spectrum Estimation W. Rose 2013-04-06

Spectrum Estimation W. Rose 2013-04-06. Signal x(t) t=0 to T Δ T=sampling interval=1/ f samp 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)

hanne
Download Presentation

Spectrum Estimation W. Rose 2013-04-06

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Spectrum Estimation • W. Rose 2013-04-06 Department of Kinesiology and Applied Physiology

  2. 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

  3. 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

  4. Example of the “simple” spectrum estimate Suppose where f=1 , d1=4/, d3=4/(3).424, d5=4/(5).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

  5. 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

  6. 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

  7. 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

  8. 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

  9. 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

  10. A Better Spectrum Estimate: Details 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 non-overlapping 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

  11. 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

  12. A Better Spectrum Estimate: Details 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

  13. Department of Kinesiology and Applied Physiology

  14. 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

  15. “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. Reason: The average power of a signal is defined as 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

  16. 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

More Related