%Waveform Filtering
%just the computations -- no GUI

clear all
figure(1)
clf
figure(2)
clf

%choose a waveform 'mepp|pulse|chirp"
%choose noise 'on|off' and the amplitude
%choose a filter type 'butter|cheby|Parks'
wavetype='chirp'
noise='off'
filtertype='parks'

%1000 samples in 1/10 second =
%sample rate of 10000 /sec
time = linspace(0,.1,1000);

%make the chosen waveform
switch wavetype
   
case 'mepp'
   %simulate a noiseless mepp
   gain = 550;    %sets the mepp amplitude
   tconst = .005; %sets the decay time
   begintime = .02   
   waveform = gain * (time-begintime) ... 
      .* exp(-(time-begintime)/tconst) ...
      .* (time > begintime);
case 'pulse'
   %simulate a noiseless mepp
   gain = 1;    %sets the pulse amplitude
   begintime = .02 ;
   endtime = .07 ;
   waveform = gain * (time>begintime & time<endtime);
case 'chirp'
   gain = 1;    %sets the pulse amplitude
   beginfreq = 10 ;
   endfreq = 5000;
   waveform = gain * chirp(time,beginfreq,.1,endfreq); 
end

switch noise
case 'on'
   %now add some gaussian noise
   noisegain = .10; %sets the noise amplitude
   waveform = waveform + ...
      noisegain*randn(size(time));
end

figure(1)
%plot the original waveform in blue
hold on
plot(time,waveform,'b','linewidth',2);

%design a filter for the mepp
%the cutoff is (bandwidth freq)/(1/2 samplerate)
%--increasing the cutoff lets more noise through the
%  filter
%--decreasing the cutoff too far chops out the mepp
cutoff = 500/(10000/2); %=0.1

switch filtertype
case 'butter'
   numSamples = 4;
   [b,a]=butter(numSamples,cutoff);
case 'cheby'
   numSamples = 4;
   %bandpass ripple for Chebychev
   ripple = .1;
   [b,a]=cheby1(numSamples,ripple,cutoff);
case 'parks'
   numSamples = 121;
   stopband = .2;
   freq = [0 cutoff stopband 1]; amp = [1 1 0 0];
   b = remez(numSamples,freq,amp);
   a = [ 1 zeros(length(b)-1,1)'];   
end

%now do the actual filtering
filterwave = filter(b,a,waveform);

%and display the filtered result
plot(time,filterwave,'r','linewidth',2)

figure(2)
%the filter frequency response
[fresponse, ffreq] = freqz(b,a,1000);
plot(ffreq/pi,abs(fresponse), 'b', 'linewidth',2);
xlabel('frequency'); ylabel('filter amplitude');