M371 - Alexiades
                      Lab 8
        Signal Processing with FFT: Compression and Denoising
            (Read it over to see what is to be done, then do it...)

---

Let's try Matlab's FFT on (samples of) the Runge function   (or can use Python: numpy.fft ).
1. Create a Matlab script Lab8a.m to do the following:
  - set number of data points: M=8 to start with,
  - set a=−5; b=5; ,
  - define a vector x of M equispaced points in [a,b]   ( x = linspace(a,b, M) )
  - and a vector y of values of the Runge function f(x)=1/(1+x²) : ( y = 1 ./ (1 + x.*x) )
  Now y is an M-point equispaced sample of (values of) f(x).
  - Compute the FFT of y: Y = fft(y)   (components are complex numbers)
  To see what the signal looks like in the "frequency domain":
  - Compute the Power Spectrum Density: PSD = YY.*conj(YY); = contribution (power) at each frequency
  - set array of frequencies: freq = 1/(b-a) * (0:M-1);
  - and plot to see signal distribution in frequency domain:
    plot( freq, PSD );
    title('Power Spectrum Density (PSD) distribution vs frequency');
    pause
  - Now take inverse FFT of Y: Yinv = ifft(Y) to 'recover' the signal in original x-domain.
  - Plot both (x,y) and (x,Yinv) on the same plot:
      plot( x,y, x,Yinv,'o' )
      title('data y and their iFFT: Yinv (o)')
  Are they the same ?
  Suppress printing of x,y,Y,Yinv ( ; at the end), and run your script for M=64 and for M=128.

2. Data Compression:
  Copy Lab8a.m to Lab8b.m and modify it to do the following (and insert appropriate comments):
  - Pick a cutoff level, say cutoff=.05 (for M=8)
  - Compute the FFT of the data vector y: YY = fft(y);
    - and zero out components of amplitude less than the cutoff,
    - count how many are set to zero (in an integer variable zeroed),
    - compute the % compression ratio cratio = 100* zeroed/M
  - Invert the filtered transform and plot to see how it's doing:
      YYinv = ifft(YY);
      plot( x,y,'k', x,YYinv,'r' )
      title('data y (black) and compressed YYinv (red)')
  - Find the L² error errYYinv between y and YYinv (what norm corresponds to L²-norm for vectors?)
  - Print them to a file out-compress nicely labeled:
      id = fopen('out-compress','a');
      fprintf( id,'# M=%i samples of Runge function \n', M )
      fprintf( id,'# cutoff   errYYinv   zeroed   cratio \n' )
      fprintf( id,' %.2f   %.2e   %4i   %.0f \n', cutoff, errYYinv, zeroed, cratio)
      fclose(id);
  a. Debug on M=16 with cutoff = 0.2, 0.5 .
  b. Turn off all unnecessary screen printing (except cutoff, zeroed, errYYinv),
    and run your script for M=128 with cutoff = 0, 1, 5, 10. Make sure printing to the file works.
  c. Looking at the plots and at the output in the out file, comment about what you observe regarding data compression via FFT.
    Which cutoff seems "best" for our (M=128) data ? in what sense? what compression ratio did it achieve ?

3. Denoising:   Now let's try some cheap filtering of "noise" in data.
  Copy Lab8b.m to Lab8c.m and modify it to do the following:
  - Perturb the data y with a small amplitude random noise: yy = y + 0.1*( rand(1,M)−0.5 );
    What does this do to the values ? Why the −0.5 ?
    What is the amplitude of the noise here ?
      plot( x,y, x,yy, 'o')
      title('original y and perturbed yy (o)')
      pause
  - Compute the L² error eryy between y and yy
  Now try to filter out the noise (smoothen the noisy signal) by zeroing out coefficients at low frequencies:
  - Compute the FFT: YY=fft(yy); ,
  - Pick a cutoff level cutoff=0.2 (for M=8),
    and zero out components of amplitude less than the cutoff, counting how many are set to zero (in a variable zeroed).
    Now YY is "filtered" or "denoised" (in frequency domain).
  - Compute the compression ratio cratio (as percent).
  - To recover the signal in x-domain, invert the filtered transform and plot:
      YYinv = ifft(YY);
      plot( x,y,'k', x,yy,'b', x,YYinv,'r')
      title('original(b), perturbed yy(g), filtered(r)')
  - Find the error errYYinv between y and YYinv, and the relative reduction of noise   relErr = 1 - errYYinv / eryy.
  - Print them to a file out-denoise nicely labeled:
      id = fopen('out-denoise','w');
      fprintf( id,'# output from Lab8c.m on denoising \n' )
      fprintf( id,'# M=%i samples, perturbed, eryy=%g \n', M, eryy )
      fprintf( id,'# cutoff   errYYinv   relErr   zeroed   cratio \n')
      fprintf( id,' %.2f   %.2e   %.2f   %4i   %.0f \n', cutoff, errYYinv, relErr, zeroed, cratio )
      fclose(id);
  a. Debug on M=16 with cutoff = 0.5, 1, 2 .
  b. Turn off all unnecessary screen printing (except cutoff, errYYinv,relErr, zeroed, cratio).
    With M=64, try cutoff = 0.5, 1, 2 to see which one seems to do best.
    Remember, the noise is random, so make several runs with each.
  c. With M=128, try cutoff = 0.5, 1, 2, 3 .
  d. Looking at the plots and at the output in the out file, comment about what you observe regarding noise reduction via FFT.
    Which cutoff seems "best" ? in what sense? what compression ratio did it achieve ?

Create a file Lab8.txt containing:
  Name:     Lab8.txt     Date:
  ==================================
  answers to 2.c
  ==================================
  answers to 3.d   (for M=128)
  ==================================
  your Lab8c.m code CLEANED up
Submit it.

---