%wienerIIR.m
%IIR wiener filter with the input signal modelled using ...
%second order Auto-Regressive (AR) model
order=4;
txsignal=sin(2*pi*16*(0:1:2047)*(1/1024))'+sin(2*pi*32*(0:1:2047)*(1/1024))'+...
sin(2*pi*64*(0:1:2047)*(1/1024))';
z=ar(txsignal(33:1:64),order);
%finite prediction error
U=fpe(z);
[p]=polydata(z)
generatedsignal=filter(1,p,randn(1,256)*sqrt(U));
figure
subplot(2,1,1)
plot(txsignal(257:1:512))
subplot(2,1,2)
plot(generatedsignal,'r')
%signal is added with the gaussian noise (variance 0.81).
rxsignal=txsignal'+randn(1,length(txsignal))*0.9;
%iir wiener filter is designed to estimate the tx signal from the rx
%signal.
SX_NUM_COEF=U*[1 zeros(1,order-1)];
SX_DEN_COEF=conv(p,p(length(p):-1:1));
L=length(conv(p,p(length(p):-1:1)));
temp=zeros(1,L);
temp(L-order+1)=1;
SY_NUM_COEF=U*temp+0.81*conv(p,p(length(p):-1:1));
SY_DEN_COEF=conv(p,p(length(p):-1:1));
SIGMARSQUARE=SY_NUM_COEF(1)/SY_DEN_COEF(1);
SY_NUM_COEF_MOD=SY_NUM_COEF/SY_NUM_COEF(1);
SY_DEN_COEF_MOD=SY_DEN_COEF/SY_DEN_COEF(1);
SYNUMR=roots(SY_NUM_COEF_MOD);
[P1,Q1]=find(abs(SYNUMR)<1);
GZ_NUM_COEF=poly(SYNUMR(P1));
SYDENR=roots(SY_DEN_COEF_MOD);
[P2,Q2]=find(abs(SYDENR)<1);
GZ_DEN_COEF=poly(SYDENR(P2));
G_ZINV_NUM_COEF=prod(SYNUMR(P1))/prod(SYDENR(P2))...
*poly(1./SYDENR(P1));
G_ZINV_DEN_COEF=poly(1./SYDENR(P2));
SXTOG_INV_NUM_COEF=conv(SX_NUM_COEF,G_ZINV_DEN_COEF);
SXTOG_INV_DEN_COEF=conv(SX_DEN_COEF,G_ZINV_NUM_COEF);
%Computing causal part of SXTOG_INV
[R,P,K]=residue(SXTOG_INV_NUM_COEF,SXTOG_INV_DEN_COEF);
[XPOS,YPOS]=find(abs(P)<1);
[B,A]=residue(R(XPOS),P(XPOS),[])
%%final transfer function
transferfunction_num=conv(GZ_DEN_COEF,[real(B) 0]);
transferfunction_den=SIGMARSQUARE*conv(GZ_NUM_COEF,real(A));
%filtering using the designed filter
b=transferfunction_num;
a=transferfunction_den;
estimatedsignal=filter(b,a,rxsignal);
figure
subplot(3,1,1)
plot(txsignal(1:1:256))
subplot(3,1,2)
plot(rxsignal(1:1:256))
subplot(3,1,3)
plot(real(estimatedsignal(1:1:256)))
figure
subplot(3,1,1)
temp1=fft(txsignal);
plot(abs(temp1(1:1:1024)))
subplot(3,1,2)
temp2=fft(rxsignal);
plot(abs(temp2(1:1:1024)))
subplot(3,1,3)
temp3=fft(real(estimatedsignal));
plot(abs(temp3(1:1:1024)))
