function [yhat H] = wienerFilter(ideal,observation,R,graphicsFlagOn,Fs); % % filtdata = wienerFilter(ideal,observation); % % FFT based Wiener filter in one dimension % % Given a ideal of our perfect underlying signal that % we wish to recover, we estimate the noise from % noise = observation-ideal; % The filtering is then performed in the frequency % domain by constructing the optimal (Wiener) filter % for this noise/ideal estimate ... H=Sf2./(Sf2/Nf2) % See Ref: Numerical Recipes in C, chapter 13.3, Press % Ref: http://www.phys.uni.torun.pl/nrbook/c13-3.pdf % % G.D. Clifford 2004 gari AT mit DOT edu % Further info: % R=real(fft(r)) where r is the smearing function % in the time domain. Default==1 ... no smearing. % To prevent a graphical display of the results, % pass an extra argument==0 % For plotting, a sampling freq Fs=256 is assumed. % You may pass a fifth arguement to change this % (although it must be integer % % Licensing: GNU GPL applies if nargin < 5 Fs=256; end Fs=round(Fs); if nargin < 4 graphicsFlagOn=1; end if nargin < 3 R=1; % response function % ... no smearing end % estimate noise from ideal noise = observation-ideal; % work out how long to make FFT N=length(observation); % Wiener filter Sf2=real(fft(ideal,N*2-1)).^2; % Smeared ideal Nf2=real(fft(noise,N*2-1)).^2; % noise Cf=real(fft(observation,N*2-1)); % ~= sqrt(Sf2+Nf2); % Corrupted ideal H=Sf2./(Sf2+Nf2); % Optimal filter Yhat=H.*Cf/R; % 'uncorrupted' ideal estimate ... yhat=real(ifft(Yhat)); % ..... in time domain % ...compensate for FFT being two sided in matlab yhat=yhat(1:length(observation)); if graphicsFlagOn==1 plotWiener(ideal, noise, observation, yhat, Fs) end