function [clean,filter]=opt_filt(spec,amp,width,slope)

% OPTIMUM FILTER function
%
% Syntax:
% CLEAN=OPT_FILT(SPEC)                 Interactive
% CLEAN=OPT_FILT(SPEC,AMP,WIDTH,SLOPE)
% [CLEAN,FILTER]=OPT_FILT(...)         Returns the FFT filter used
%
% This function cleans a regularly-spaced data series of high frequency noise
% by low-pass filtering the Fourier transform of the data.  The low pass filter
% is an optimal filter in the sense of Brault and White (1971, Astronomy and
% Astrophysics, Vol 13, p. 169).
%
% CLEAN is the cleaned version of SPEC
% AMP, WIDTH, and SLOPE are input parameters if the shape of the optimal
% filter to use is known, ahead of time.  Otherwise, the routine is
% interactive.
%
% The relationship between CLEAN and SPEC is given by the following lines of 
% code:
% CLEAN=real(ifft(fft(SPEC).*filter));
% CLEAN=CLEAN+mean(SPEC-CELAN);
%
% Eric Tittley (95 07 10) etittley@phobos.astro.uwo.ca

if (nargin~=1 & nargin~=4) error('Incorrect number of arguments passed to OPT_FILT'), end
if nargin==4 interact_flag=0;
else interact_flag=1;
end

spec=spec(:)';  % Make sure it's a row vector.

n=length(spec);
m=ceil(n/2);
o=ceil(n/4);

wave=[1:m];
wave2=[1:n];

power=abs(fft(spec));

noise=exp(mean(log(power(o:m+o))))*(1&wave2);
if interact_flag
 disp('First we will set the noise level.')
 ok=0;
 while ok~=1
  plot(wave,log([power(1:m);noise(1:m)]))
  axis([0 m -5 5])
  ok=input('Good enough? (1 for yes, else a new value): ');
  if ok~=1 noise=exp(ok)*(1&[1:n]); end
 end

 disp('Good. Now lets find the region of signal.')
 disp('Three parameters are required.')
 disp('1 The Amplitude')
 disp('2 The width')
 disp('3 The slope')
 disp('The last two affect each other.')
 amp=power(2);
 width=ceil(n/22);
 slope=1.2;
 signal=amp*(1&wave2)./exp((wave2/width).^slope);
 ok=[0,0];
 maxy=ceil(max(log(power)));
 miny=floor(min(log(power)));
 while length(ok)>1
  plot(wave,log([power(1:m);noise(1:m);signal(1:m)]));
  axis([0 m -5 5])
  disp(sprintf('[amp=%3.1f,width=%4.0f,slope=%3.1f]',log(amp),width,slope));
  ok=input('New params [amp,width,slope] or 1 if ok: ');
  if length(ok)==3
   amp=exp(ok(1));
   width=ok(2);
   slope=ok(3);
   signal=amp*(1&wave2)./exp((wave2/width).^slope);
  end
 end
else
 signal=exp(amp)*(1&wave2)./exp((wave2/width).^slope);
end


signal(2:n)=signal(2:n)+signal(n:-1:2);
filter=signal./(signal+noise);
clean=real(ifft(fft(spec).*filter));
clean=clean+mean(spec-clean);