%This compares deterministic and stochastic deconvolution of a Wiener process
%from observation of its summation (discrete time) or integration (continuous time)
%The deterministic Wiener filter fails. The stochastic Wiener filter works well.
%Note: signals differ from text, since all signals are randomly generated.
clear;dt=0.001;sigma=0.01;
X=cumsum(randn(1,1000))*dt;%Original Wiener process.
H=ones(1,101);%Impulse response of summation.
HX=conv(H,X);N=sigma*randn(1,1100);Y=HX+N;%Noiseless and noisy convolved Wiener process.
FY=fft(Y);FH=fft(H,1100);%Zero-pad DFTs to lengths of convolved Wiener process.
FX=abs(fft(X,1100)).^2/1100;SX=sum(X.^2);%S_w(f) for stochastic and deterministic.
FZ=conj(FH).*FX./(abs(FH).^2.*FX+sigma^2).*FY;Z=real(ifft(FZ));%Stochastic
FW=conj(FH).*SX./(abs(FH).^2*SX+sigma^2).*FY;W=real(ifft(FW));%Deterministic
figure,subplot(211),plot(HX(1:1000)),title('Convolved (summed) Wiener process')
subplot(212),plot(Y(1:1000)),title('Noisy Convolved (summed) Wiener process')
figure,subplot(211),plot(X),title('Noiseless original Wiener process')
subplot(212),plot(Z(1:1000)),title('Deconvolved using stochastic Wiener filter')
figure,subplot(211),plot(X),title('Noiseless original Wiener process')
subplot(212),plot(W(1:1000)),title('Deconvolved using deterministic Wiener filter')

