### spectrum_estimation

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