% Script for DFT low-pass filter (based on Timothy Sauer's book "Numerical
% Analysis")
clear all
clc
close all
% Read the data from a standard .wav file
fprintf('Reading data...\n')
[y,Fs,nbits] = wavread('curve');
fprintf('The sample rate is %d Hz, ', Fs)
fprintf('and is sampled using %d bits/sample.\n', nbits)
N = length(y);
c = y;

soundsample = c(1:2^16);   % good to take power of 2 since m needs to be even
n = length(soundsample);   % n is the length of signal

modeskept = 1;
while modeskept
    modeskept = input(sprintf('Enter fraction of Fourier modes to keep (power of 2 <= %d, 0 to stop): ',n));
    if modeskept == 0 break; end
    m = n/modeskept;          % an even integer <= n, the number of Fourier modes kept

    FTsample = fft(soundsample);                  % compute interpolation coefficients
    truncatedFT = zeros(n,1);             % will hold coefficients for ifft
    truncatedFT(1:m/2) = FTsample(1:m/2);        % keep only first m frequencies 
    truncatedFT(m/2+1) = real(FTsample(m/2+1));  % since m is even, keep cos term only
    if(m < n)                  % unless at the maximum frequency,
      truncatedFT(n-m/2+1) = truncatedFT(m/2+1);   %  add complex conjugate to
    end                        %  corresponding place in upper tier
    truncatedFT(n-m/2+2:n) = FTsample(n-m/2+2:n);% more conjugates for upper tier
    filteredsample = real(ifft(truncatedFT));   % invert fft to recover data
    disp('Original sound')
    sound(soundsample,Fs)
    disp('Filtered sound')
    sound(filteredsample,Fs)
    t = (0:n-1)/Fs;
    figure
    plot(t,soundsample,t,filteredsample)
end