The Discrete Fourier Transform and Noisy Signals

The discrete Fourier transform and noisy signals The objective of this lab is to explore how to uncover a signal buried in noise by manipulating it in the frequency domain via the discrete Fourier transform. Recall that the definition of the discrete Fourier transform for a discrete signal {x(n)} ∈ SN is x̂(k) = N −1 X x(n)e−2πikn/N for k = 0, 1, . . . , N − 1. (1) n=0 In Matlab, we compute the discrete Fourier transform efficiently using fft (fast Fourier transform) and ifft (inverse fast Fourier transform). Exercise 1 (DFT and shifts) As in the previous lab, use Matlab to construct a signal of length 213 = 8192 that contains only three frequencies: x=0.5*sin(2*pi*200*t)+0.2*sin(2*pi*455*t)-0.3*sin(2*pi*672*t) Compute and plot the absolute value of the DFT, also as in the previous lab. Next shift the signal by replacing t with (t-t0): xshifted=0.5*sin(2*pi*200*(t-t0))+0.2*sin(2*pi*455*(t-t0))-0.3*sin(2*pi*672*(t-t0)) Using different values of t0 like 0.1 or 0.4, compute and plot the absolute value of the DFT for the shifted signal. What changes? Submit a figure with 3 subplots with the top plot showing both the original and shifted signals, and the two lower plots showing the DFT of each signal. Clearly label axes, etc. Exercise 2 (Detecting a signal above the noise) Add noise to your signal from Exercise 1: xnoisy=x+randn(size(x)); Plot the noisy signal and then find and plot the absolute value of its DFT. We can’t make out the signal very well in the time domain–it looks buried in noise–but it’s clear in the frequency domain. The ”white noise” we added has a uniform distribution in frequency, so we can try to remove it by keeping only the DFT coefficients whose absolute value rises above the noise threshold. Examine your DFT plot and determine a threshold value (that most of the noisy DFT coefficients fall below). Then clean up the signal as follows: Xcleaned=X.*(abs(X)>=threshold); % set to zero if below the threshold value signal=real(ifft(Xcleaned)); % take real part in case some residual imaginary part Submit a figure with 3 subplots showing the noisy signal, the absolute value of its DFT, and the recovered signal plotted with the original (non-noisy) signal (set xlim to give a clear view in each subplot, label the axes, and add a legend in the 3rd subplot). Exercise 3 (Tracking a moving object embedded in noise) The last task is to uncover a signal hidden in noise with an extra twist: suppose we have observations of the temperature along a 60cm-circumference ring taken at different times. We know there is a hot spot on the rod that is moving back and forth, and we want to know where the hot spot is at each time. Download the file TrackingData.csv and load it into Matlab: U=dlmread(’TrackingData.csv’,’,’). Initialize some variables we need as follows: N=512; L=60; t=1:20; freq=[0:(N/2-1) -N/2:-1]/L; x=L/N*[-N/2:N/2-1]’; [T,S]=meshgrid(x,t); % look up meshgrid in the Matlab help to see what it does [F,S]=meshgrid(freq,t); The matrix M in this case is a matrix whose rows contain the temperature data for equally spaced positions along the ring for each time point. Plot the (noisy) temperature profiles at each time: figure; waterfall(T,S,U) xlabel(’Position along ring (cm)’); ylabel(’Time’); title(’Noisy signals’) Click on the symbol at the top of the figure window that looks like an arrow bending counterclockwise; this allows you to move the waterfall plots around so you can get a full 3D feeling. Your job is to find the signal buried in the noise. Experiment by plotting the summed temperature profiles to see if the noise “averages out” well enough to see what’s happening: figure;plot(x,mean(U)) Summing in the frequency domain usually works better in this setting, because we’re expecting a similar pattern in each temperature profile, but shifted (as the hotspot changes its location on the ring). From Exercise 1, we know this doesn’t change the position of the dominant frequencies in the DFT, so we have a better chance of reinforcing the desired frequencies while averaging out the noise-induced frequencies by summing in the frequency domain: for k=1:20 Uhat(k,:)=fft(U(k,:)); end figure;plot(fftshift(freq),fftshift(abs(mean(Uhat)))) Use this figure to determine which frequencies correspond to the signal that we want the filter to keep nonzero. The brute force approach of a step function is not ideal here; using a Gaussian will work better to more gently taper off the undesired frequencies: filter=exp(-damp*freq.^2). Choose a value for damp so that when you plot this Gaussian filter on the plot with the mean Uhat values, it covers the frequency band of interest (you can plot the filter multiplied by the max value of mean(Uhat) to get a better view). Use a for-loop as above to generate a “cleaned” version of Uhat, multiplying each Uhat(k,:) by the Gaussian filter, then take the inverse transform to obtain the cleaned up signals. For your third figure to submit, make a waterfall plot of the cleaned up signals, rotated for a good view of the locations of the hotspots (which should look like red peaks).