% Function to do instrumental profile deconvolution with
% either Gaussian or Voigt noise filtering.  Uses cbell.m
% for cosine bell weighting.  To use:
%  d = decnv(r,p);
% r = stellar line profile data with continuum = 0.
% p = one half of symmetric instrumental profile.
function d = decnv(r,p);
% Set up filter using data power spectrum...
nfft = input('Enter desired FFT length: ');
Ps = abs(fft(r,nfft));
ampl = Ps(1); n2 = nfft/2;
disp('Amplitude = '); disp(ampl);
disp('Set the noise level by clicking on left mouse button...');
semilogy(Ps), xlabel('Frequency'), ylabel('Power')
[ix,Pn] = ginput(1);
disp('Noise power level = '); disp(Pn);
% Filter choice doesn't matter much...
ifilt = input('Gaussian or Voigt filter [1,2]? ');
disp('Set signal power...');
ok = 0; sig = (0:1:n2-1)'; x = (0:1:nfft-1)';
phi = zeros(n2,1); filter = zeros(nfft,1);
while ok ~= 1
 if ifilt == 1,
  alpha = input('Enter filter parameter alpha: ');
  ba = (Pn/ampl)^2; phi = 1 ./ (1 + ba .* 10 .^ (2 .*(alpha.*sig).^2));
 end
 if ifilt == 2,
  c1 = input('Enter filter parameter c1: ');
  c2 = input('Enter filter parameter c2: ');
  Psig = ampl .* exp(-(c1*c2).*sig - (c2.*sig).^2);
  phi = 1 ./ (1 + (Pn./Psig).^2);
 end
 filter(1:n2) = phi; filter(nfft-n2+2:nfft) = flipud(phi(2:n2));
 plot(x,fftshift(filter),x,fftshift(Ps)), xlabel('Frequency'), ylabel('Power')
 ok = input(' Is this OK [1,0]? ');
end
% Set up instrumental profile transform...
z = zeros(nfft,1); np = length(p);
p = p ./ sum(p);
z(1:np) = p; z(nfft - np + 2:nfft) = flipud(p(2:np)); Z = fft(z);
% Set up stellar data transform...
nr = length(r); w = cbell(nr); ave = mean(r); R = fft(w.*(r - ave),nfft);
d = real(ifft(filter .* R ./ Z));
x1 = (0:1:nr-1)';
d = ave+d(1:nr); disp('Deconvolved profile...');
plot(x1,r,x1,d), xlabel('Wavelength'), ylabel('Residual Flux')