%This implements denoising a noisy Wiener process using 
%(a) deterministic and (b) stochastic Wiener denoising filters.
%The deterministic denoised signal is just a scaled version of the noisy signal.
%The stochastic denoised signal is a smoothed version of the original signal.
%Note: signals differ from text, since all signals are randomly generated.
clear;dt=0.001;fmax=1/dt;df=1;sigma=0.01;%df=1/tmax=1;
X=cumsum(randn(1,1000))*dt;%Original Wiener process
N=sigma*randn(1,1000);Y=X+N;%Noisy Wiener process
t=[0:dt:0.999];f=[1:df:1000];FY=fft(Y);
figure,subplot(211),plot(t,X),title('Noiseless Wiener process')
subplot(212),plot(t(401:500),X(401:500)),title('Segment of noiseless Wiener process')
FX=abs(fft(X)).^2/1000;%Periodogram estimate of psd S_w(f) for Wiener process.
figure,plot(log10(f(2:400)),log10(FX(2:400)),...
log10(f(2:400)),log10(1/4/pi/pi./f(2:400).^2),'r'),...
title('Plot of log of periodogram estimate of S_w(f) vs. log(f). Slope is -2.')
figure,subplot(211),plot(t,X),title('Noiseless Wiener process')
subplot(212),plot(t,Y),title('Noisy Wiener process')
FZ=FX./(FX+sigma^2).*FY;Z=real(ifft(FZ));%Stochastic Wiener filter.
SX=sum(X.^2);Q=SX/(SX+sigma^2)*Y;%Deterministic Wiener filter (in time domain).
figure,subplot(211),plot(t,X),title('Noiseless Wiener process')
subplot(212),plot(t,Z),title('Denoised using stochastic Wiener filter')
figure,subplot(211),plot(t,X),title('Noiseless Wiener process')
subplot(212),plot(t,Q),title('Denoised using deterministic Wiener filter')
