signal processing notes

APR 2003
symbols
going from continuous to discrete fourier transform
sin(x)/(x) and its properties
polyphase filter banks
gaussian width after convolving 2 gaussians
Converting spectrum back to rms volts
Abbreviations/Symbols used
ts sample  time (time between samples)
fs
fn		 sampling frequency (1/ts)
nyquist frequency  fs/2=1/(2*ts)
N total number of time samples used.
fres

 

 frequency resolution: 1/(ts*N)
# of independent samples: N/2 . If complex samples are taken then each complex number is two independent numbers so there will be N independent frequency channels.
FT Fourier Transform
x(t),X(t) x(t) is the time series, X(f) is the FT of x(t)
i sqrt(-1)
m,M
j,J etc..		 if m,j are indices of sets (say to sum over), then M,J will be the number of elements in the respective sets.

Going from a continuous to discrete fourier transform:  (top)

    The continuous FT of the function x(t) is: X(f)= integral(x(t)*exp(-2*pi*i*f*t)*dt). (f=freqHz, t=timeSecs)
When going to the discrete, finite case of N samples we get:
t=>ts* n      n=0..N-1 since the samples are spaced by ts.
f=> k/(ts*N)     k=0..N-1.  N/2 samples are duplicated for real input samples.
The FT becomes:
X(k/(ts*N))= sumn(x(n*ts)*exp(2*pi*i*(k/(ts*N))*n*ts) * dt) n=0..N-1 or
X(k/(ts*N))= sumn(x(n*ts)*exp(2*pi*i*k*n/N) * dt) n=0..N-1  .
If we agree that Xk is in steps of our frequency resolution: 1/(tsN),  that xn are the time points sampled at n*ts , and that dt is our unit of time then we can simplify the equation to:
Xk = sumn(xn*exp(2*pi*i*kn/N))  n=0..N-1, k=0,N-1

sin(x)/x and its properties:  (top)

    When you transform a data set of N points (with no explicit windowing) then you are weighting the data time points with 1 and all other time points with 0. This is a boxcar window of duration N.  sin(x)/x is the transform of a boxcar window.
    To show that sin(x)/x is the FT of a boxcar window,  we  integrate the window from -T/2 to T/2 (to  make it symmetric  so the anit-symmetric imaginary part goes to zero):
X(f)=integral(exp(-2*pi*i*f*t)dt)  -T/2 <=t<=T/2 ..
The imaginary part is odd so it integrates to 0. The real part integrate to:
X(f)= 1/(2*pi*f) * [ 2*sin(2*pi*f*T/2)] = sin(pi*f*T)/(pi*f)= T* ( sin(pi*f*T)/(pi*f*T)
If we change to the discrete format:
Xk=N * sin(pi*(k/N)*N)/(pi*k/N*N) = N*sin(pi*k)/(pi*k)
The zeros of sin(x)/x are at:
zeros: sin(x) = 0  at x=(n+1)*pi (n=0..)
The extrema of sin(x)/x are found by differentiating and setting it to zero:
d/dx (sin(x)/x) = (cos(x)/x - sin(x)/x^2)=0
x^2*cos(x)=x*sin(x) or x=tan(x)
Tan(x) goes to infinity at x=(2n+1)*pi/2. As x goes to infinity the extrema must approach (2n+1)pi/2 so for large x this is a pretty good approximation for the peaks. Looking at (sin(x)/x)^2 (the power) the difference between (2n+1)*pi/2 and the correct value is: (1st: .2 db, 2nd: .07db, 3rd: .03db,4th: .02db...)

The sin(x)/x plot for (x>0) has  a linear plot of the voltage on top and a db plot of the power (sinx/x)^2 below. Since both sin(x) and (x) are anti symmetric in (x), sin(x)/x is symmetric and x<0 is the same as x>0.
    The first sidelobe is -13.26 db down in power from the peak. The values for the first 14 power sidelobes are printed on the bottom plot.

processing:x101/polyphase/sinx.pro

polyphase filter banks  (top)

    A filter bank will take a single bandwidth BW and channelize it into K narrower channels BN. If K is 1 then the filter bank is a normal lowpass/bandpass filter.
    The simplest multi channel filter bank is created by fourier transforming (FT) the data. We take K input points and perform a FT. This uses a boxcar window in the time domain. The convolution theorem tells us that this will convolve the channels in the frequency domain with sin(x)/x (since sin(x)/x is the FT of a boxcar window).  After computing the FT we need to square the data to compute the power in each channel BNk. The zeros and maxima of sin(x)/x and (sin(x)/x) ^2 are shown above.
 

    The individual filters BNk (k=0,K-1) will have a width that depends on how the filter bank was constructed. The filter BNk will have some overlap with adjacent filters.  When a strong signal falls in filter BNj, it will also appear (almost always) at a reduced level in BNk. This overlap in response is called the filter leakage. The polyphase filter technique allows you to decrease this leakage at the expense of taking larger time sequences.

straightforward way to decrease filter leakage.

    The easiest way to decrease the filter leakage is to compute a spectra with more frequency resolution than needed, and then smooth and decimate this spectrum to the desired resolution. To be more specific lets assume that we want 1024 frequency channels on output. To get better frequency resolution you need to take a longer run of data. Instead of 1024 points lets take 4096 points (4 times the needed points). We do a 4096 point FT to go to the frequency domain. We then define a 4 point filter that will be used to smooth the 4096 points spectra and reduce the filter leakage. After smoothing we can then decimate every 4th point to get back to the desired resolution. Writing this out using * as the convolution operator, x for mulitiplication, C4 for the 4 filter coefficients, y(t) as the 4096 time points, and Dn for decimate by n we get:
    D4(FT4096(y(t)) * C4 )

A little harder but quicker:

     Since the FT is a linear operation we can think of doing the smoothing/decimating operation before the FT. The smoothing (in the frequency domain) is a convolution so it becomes at multiply by the FT of C4 in the time domain.  This is just a slowly varying function of time (since there are no high freq components with only 4 coefs). The decimation by 4 is a multiply by a set of delta functions spaced by 4 points. In the time domain this becomes a convolution by delta functions spaced by 1024 points. Writing this out and using c4096 as the 4096 transform of C4 , and Delta1024 as a set of delta functions spaced by 1024 points we get:
    (y(t) x c4096 )  * Delta1024  . After doing this we have our 1024 points (still in the time domain) with the requested bin shape. We then do a FThh1024  transform to go the the frequency domeain.

   processing: x101/polyphase/polyphase.pro

gaussian width after convolving 2 gaussians  (top)

    Convolving two gaussians together gives a gaussian of a larger width. The plot shows how the gaussian width changes as a function of the input gaussians relative widths.
processing: x101/misc/gsconvolve.pro

Converting spectrum to rms volts.  (top)

    Given a spectral density function (spectrum) compute what the rms input voltage was. This was done for a sine wave and noise (generated on the computer). This is handy when you want to know the rms counts at a digitizer given the average value of the spectrum. The spectrum was computed from the fft using N=4096 pnts and the "forward transform" formula:
spc(f) = (abs( 1/N*St(v(t)*exp(-2pi*f*t/N)))^2
   the forward fft is getting divided by 1/fftlen in the idl routine.
 The data was set so that the rms value was 10 counts. Below i refer to I and Q. These are the inphase (Real) and Quadrature (imaginary) inputs.
The plots show going from spectrum back to rms volts (.ps) (.pdf):

Conclusions:

    processing: x101/070926/spcToRms.pro

page_up
home_~phil