GG 313 Lecture 25 Transform Pairs and Convolution 11/22/05

1. GG 313 Lecture 25 Transform Pairs and Convolution 11/22/05

2. Today we continue our look at some important transform pairs and introduce convolution. Recall from last class:

3. Lets look at the sine wave transform. It consists of two peaks of the same amplitude, one at positive and one at negative frequency, one positive and one negative. What would this look like in Matlab? % sintran sine transform sfreq=20; % Hz len=256; %signal length dt=.01; % sample period dfreq=(1/dt)/len; for ii=1:256 time(ii)=(ii-1)*dt; % time scale freq(ii)=(ii-1)*dfreq; % frequency scale if ii>len/2+2; % if above the Nyquist, start frequency back down freq(ii)=freq(len+2-ii); end y(ii)=sin(2*pi*sfreq*time(ii)); end

4. plot (time,y) This sine wave has a frequency of 20 Hz and an amplitude of 1. It is samples at 100 samples/sec (1/5th of the sample period, or 2/5th of the Nyquist frequency.

5. Y=fft(y); p=Y.*conj(Y); % power nfig=figure; plot(p) The raw power spectrum shows two spikes, one at 1/5 and one at 4/5 of the sampling frequency.

6. plotting against frequency: nfig=figure; plot (freq,p) Hz What is the amplitude of the peak? The value at the peak is 14,264. What’s it’s significance?

7. Recall the delta function transform. The amplitude of the transform across all frequencies is the square of the height of the delta function, and the peak height of the sync function (transform of the boxcar) is the height of the boxcar squared times the number of non-zero points. So, what might we expect for the transform of the sine wave? First of all, note that This is an important fact (an another way to check on the behavior of the transforms you are using). This says that the sum of the squares of all the points in the time domain times the length of the time series is equal to the sum of the power in the frequency domain. Check it: in our example, len*sum(y.^2)=sum(p)=32,640.

8. Well, that’s fine, but shouldn’t all the energy be at the sine wave frequency? Why aren’t all the values away from dfreq equal to zero, and the amplitude of the peak equal to 32,640/2 (the divide by two because there are 2 peaks)? Check the value of dfreq in terms of the spectral estimates. Note that 20 Hz occurs at sample 256/5=51.2. This means that the frequency of our sine wave is not a perfect multiple (harmonic) of the sampling period, and thus an even number of cycles don’t fit in the length of the time series. Look at the end points plotted so that we butt the end of the previous series with the start of the next:

9. Plot([y(250:256) y(1:7)] If n cycles of our sine wave fit exactly into the length of the sample, then this glitch in the red circle wouldn’t be there, and ALL the energy would be in the peak at the correct frequency.

10. Let’s try a sine wave frequency that fits exactly in the sample. Recall the fundamental frequency is the one where one cycle fits in the sample. =2*pi*1/T. Where T=len*dt. To about match the frequency in our example, we want about 50 times this frequency, or  =100*pi/(256*.01), and f=50/T19.53 Hz. When we run the m file using this frequency, we find that ALL of the energy is exactly at this frequency (split evenly between the two peaks). If we again plot the end points butted together, we see that the “glitch is gone:

11. The lesson here is that signals that contain energy at frequencies that are NOT at harmonics will show “leakage” of those frequencies into other frequencies. This happens all the time, but it isn’t a serious problem, and it is one we can pretty well fix. ________________ Note the big difference between the transform pairs we’ve looked at. The first two were single occurrences of a TRANSIENT, and not periodic signals at all. The sine wave pair is a periodic signal in the time domain. Transients in the time domain lead to periodic signal in the frequency domain and visa versa. If we want to know the size of a transient in the time domain from its transform, we can get it very easily, but to get the size of a periodic signal, we need different constants.

12. How do we get the amplitude of a sine wave. Certainly, we can just measure the height of the signal in the time domain, but what if there are several sine waves that interfere with each other in the time domain. We can get the exact amplitude of sine waves by applying the following factor: p=2*sqrt(p/(len^2)); where p=Y.•conv(Y); but only where the sine wave(s) are exact harmonics. If not, the amplitude will be les than the true amplitude since the energy in the sine wave is smeared out. Check out the examples on the next page:

13. Perfect harmonics, amplitudes 2 and 1 arbitrary frequencies, amplitudes 2 and 1

14. Power Spectral Density The power spectrum is fine for transient analysis and pure sine waves, but when we’re dealing with continuous signals and noise, the amplitude units of the power spectrum depend on the number of estimates in our sample, and thus don’t give a “universal” estimate of our signal amplitude. So we use the PSD. PSD, or power spectral density is defined as the power per unit frequency, often the power per Hz, where power has the units of the square of the amplitude in the time domain. For example if our units of amplitude in the time domain are meters, then the psd units are meters2/Hz. The PSD is calculated as a power spectrum, and then each spectral estimate is divided by the len/bandwidth of each estimate. The bandwidth of each sample is the Nyquist frequency divided by the number of spectral estimates, with a correction dividing the Y(0) and Y(Nyq) estimates by 2.

15. The power spectral density thus yields an estimate of the power in the signal vs frequency that is independent of the nyquist frequency and of the number of samples. We can easily calculate the PSD in Matlab using the “periodogram” function, which returns the PSD: [Pyy,w]=periodogram(y,[ ],len,1/dt) [psd estimate vector, frequency vector]= periodogram(input time series, window,length of time series, sampling frequency) What’s the PSD of our example sine wave? Interestingly, this power DOES depend on the length of the time series, and thus the bandwidth. This is because our sine wave has a bandwidth less than the bandwidth of the spectral estimates, thus the PSD smears the energy over the bandwidth when it shouldn’t.

16. Frequency The figure above tries to explain some of these concepts. An 8-point transform is shown. The blue area is the region from zero frequency to the Nyquist frequency. The first spectral estimate is at 0 Hz with a bandwidth that goes from -N/8 to N-8. Since only half of this bandwidth is in the blue area, the power in that bandwidth should be divided by 2. The same goes for the estimate at the Nyquist frequency.

17. A sine wave at an arbitrary frequency will show up as two peaks, one at positive and one at negative frequencies (you can move the right-hand-side of the spectrum to negative frequencies). the sine wave does not have a bandwidth, and thus the PSD, which calculates the energy in the bandwidth of a spectral estimate, will yield a too-small value for the amplitude of the sine wave.

18. Before continuing to the last transform pair, we need to study CONVOLUTION. Convolution is a very important process by which FILTERS are applied to signals. A filter is an operation that changes the shape of a spectrum and/or applies a delay to a time series (changes the phase spectrum). An AMPLIFIER (or ATTENUATOR, if the amplitude is decreased) changes only the amplitude of a signal without changing the shape of the spectrum. Convolution operates much like the correlation function:

19. The output of a convolution is the sum of the products of the input signal (x) multiplied by the lagged filter function (h). This is written hx=y, where the ‘’ sign means convolution. Diagramatically:

31. Note that you need to be careful how you apply the indices if h(t) is not symmetric. If you put it in backwards, you get the wrong answer. Convolution and Fourier transforms go hand-in-hand. This is because the Fourier transform transforms the convolution operation into plain multiplication. This is the Convolution Theorem. In words, I can convolve two functions in the time domain or take their Fourier transforms and MULTIPLY them together in the frequency domain. So what? Doing the Fourier transform is actually faster than doing a long convolution. I believe the beak-even point is about a 32-point convolution.

32. Matlab uses the function “conv(h,x) to convolve two series. EXERCISE: Fire up Matlab and convolve the two functions used in the earlier example. h=[0 1 -2 1 ] and x=[0 -1 3 0 -3 -1 -2 -1 0] Instead of the h above, convolve the impulse I=[0 0 1 0 0] with x. What do you get? Convolve h with a boxcar: b=[001111111111100] What do you get?