SpectralLeakage

Download Report

Transcript SpectralLeakage

Spectral Leakage
Pp289 jdsp
Freq of kth sample, No centering
kf s
k
fk 

N t N
Centering and the PSD
• PSDs are centered if the DC component is
at N/2, where N= number of samples.
• Before and after centering
How do I center the FFT?
UlawCodec ulc = new UlawCodec();
•
•
final double[] doubleArray = ulc.getDoubleArray();
double d[] =
Mat1.truncateToIntegralPowerOfTwo(doubleArray);
•
•
•
•
•
•
FFT1dDouble f = new FFT1dDouble();
double in_im[] = new double[d.length];
Mat1.centering(d);
f.computeForwardFFT(d, in_im);
f.normalize();
final double[] psd = f.getPSDNotNormalized();
•
Graph.graph(psd, "", "", "psd");
•
•
•
•
f.computeBackwardFFT(d, in_im);
Mat1.centering(d);
UlawCodec ulc2 = new UlawCodec(f.getRealData());
ulc2.play();
How do I implement centering?
• public static void centering(double r[]) {
•
int s = 1;
•
for (int i = 0; i < r.length; i++) {
•
s = -s;
•
r[i] *= s;
•
}
•
}
How do I compute freqs?
fsk fs
fk 

N
2
f s  8000, N  8000, k  4000
8000* 4000 8000
fk 

0
8000
2
Computing Freqs
• The center freq
should be zero.
• But what about 4001
and 3999?
8000* 4001 8000
fk 

 1hz
8000
2
8000*3999 8000
fk 

 1hz
8000
2
Thus every bin gives one freq
•
•
•
•
Bins give a unique freq.
K=n/2=0 hz
K=n = 4000 hz
K=0= -4000 hz
Real signals
• Real signals have zero complex
components
• Real signals have symmetrical PSDs
• Real signals are positive on the right
• Real signals are negative on the left
• One bin maps to one frequency
Spectral leakage
• frequency analysis effect
• Applies to:
– finite-length signals
– finite-length segments
• energy has "leaked" out of the original
signal spectrum into other buckets.
Windowing
• window function (also known as an
apodization function or tapering
function)
• a function that is zero-valued outside of
some chosen interval
Windowing
• windowing is the root cause of spectral
leakage
Does Windowing Really Help?
• Non-rectangular window functions actually
increase the total leakage
• they can also redistribute it to places
where it does the least harm, depending
on the application
• windowing predistorts the input samples
so that the spectral leakage is evened out
For example
k
fk 
N t
k
fk 
8000  10k
800
k  Nfk / f s
Frequency Bin
• Often the kth array element in the
frequency domain is referred to as a
frequency bin
The 400 Hz spectrum
• we expect the maximum amplitude to
occur at.
k  Nfk / f s  2048*400 / 8000 =102.4
k=102, and 103<- energy spreads around two bins
Spectral Leakage
• When a single freq leaks into more than
one bucket in the fft.
• This occurs because of rectangular
windowing function use in a STFT
• It smearing the PSD across multiple
buckets.
To limit spectral leakage
• Use a non-rectangular windowing function
• A function that is limited in time and zero
everywhere else.
• Finite support.
• An a-periodic function
Select a windowing function
•
•
•
•
Use either speed
Or accuracy
Or ease of implementation
Do a straight multiply of the window and
your input signal.
• You may need to do this several times, if
you don’t want to loose data…
But I need speed!
• If you don’t overlap, you will loose data,
but you will gain speed.
• You need to know the nature of the input
signal!
There are many windows
•
•
•
•
•
Hanning
Blackman
Bartlet
Hamming
Etc.
Windows
Code for Hanning
• public static double[] makeHanning(int n) {
•
double window[] = new double[n];
•
double arg = 2.0 * Math.PI / (n - 1);
•
for (int i = 0; i < n; i++)
•
window[i] = 0.5 - 0.5 *
Math.cos(arg * i);
•
return window;
•
}
Lyon Window
Lyon Window Code
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
public static double[] makeLyon(int n) {
double window[] = new double[n];
double u = 0.0;
double du = 2.0 / n;
for (int i = 0; i < n / 2; i++) {
window[i] = y5(0, 1.0, u);
u += du;
}
u = 0;
for (int i = n / 2; i < n; i++) {
window[i] = y5(1.0, 0.0, u);
u += du;
}
return window;
}
y5
• public static double y5(double y0, double
y1, double u) {
•
double t2 = u * u;
•
double t3 = t2 * t2;
•
return (6 * y1 - 6 * y0) * t3 * u +
•
(-15 * y1 + 15 * y0) * t3 +
•
(10 * y1 - 10 * y0) * t2 * u + y0;
•
}
Lyon vs Hanning Window