How can we make this code shorter?

1 view (last 30 days)
Alice Zurock
Alice Zurock on 13 Apr 2020
Commented: Alice Zurock on 13 Apr 2020
function APSM
figure
L=200;%Dimension of the unknown vector
N=3500; %Number of Data
theta=randn(L,1);%Unknown parameter
IterNo=100;
MSE1=zeros(N,IterNo);
MSE2=zeros(N,IterNo);
MSE3=zeros(N,IterNo);
MSE4=zeros(N,IterNo);
noisevar=0.3;%For the high noise scenario set this value equal to 0.3
epsilon=sqrt(2)*sqrt(noisevar); %Epsilon parameter used in the APSM
X = randn(L,N);
inputvec = @(n) X(:,n);
noise = randn(N,1)*sqrt(noisevar);
y = zeros(N,1);
y(1:N) = X(:,1:N)'*theta;
y = y + noise;
for It=1:IterNo
It
xcorrel = randn(N+L-1,1);
xcorrel = xcorrel/std(xcorrel);
X = convmtx(xcorrel,L)';%Generate noise process
X(:,1:L-1) = [];
inputvec = @(n) X(:,n);
noise = randn(N,1)*sqrt(noisevar);%Generate noise process
y = zeros(N,1);
y(1:N) = X(:,1:N)'*theta;
y = y + noise;
w=zeros(L,1);%Estimate vector
mu=0.2;
delta=0.001;
q=30;% Number of window used for the APA and the APSM.
for i=1:N
if i>q
qq=i:-1:i-q+1;
yvec=y(qq);
Xq=inputvec(qq);
Xq=Xq';
e=yvec-Xq*w;
eins=y(i)-w'*inputvec(i);
w=w+mu*Xq'*inv(delta*eye(q)+Xq*Xq')*e;%APA update
MSE1(i,It)=eins^2;
end
end
w=zeros(L,1);
for i=1:N
if i>q
qq=i:-1:i-q+1;
yvec=y(qq);
Xq=inputvec(qq);
Xq=Xq';
sum1=zeros(L,1);
sum2=0;
for jj=1:q
p=projection(yvec(jj),Xq(jj,:),w,epsilon);%Compute projection onto hyperslab
sum1=sum1+(1/q)*p;
sum2=sum2+(1/q)*(((p-w)'*(p-w)));
end
if sum1==w%Extrapolation Parameter Computation
Mn=1;
else
Mn=sum2/((sum1-w)'*(sum1-w));
end
eins=y(i)-w'*inputvec(i);
w=w+Mn*0.5*(sum1-w);%APSM update
MSE2(i,It)=eins^2;
end
end
w=zeros(L,1);%RLS recursion
delta=0.001;
Pi=delta*eye(L);
P=(1/delta)*eye(L);
for i=1:N
gamma=1/(1+inputvec(i)'*P*inputvec(i));
gi=P*inputvec(i)*gamma;
e=y(i)-w'*inputvec(i);
w=w+gi*(e);%RLS update
P=P-(gi*gi')/gamma;
MSE3(i,It)=e^2;
end
w=zeros(L,1);%NLMS Recursion
delta=0.001;
mu=1.2;
for i=1:N
e=y(i)-w'*inputvec(i);
mun=mu/(delta+inputvec(i)'*inputvec(i));
w=w+mun*e*inputvec(i);%NLMS update
MSE4(i,It)=e^2;
end
end
MSEav1=sum(MSE1')/IterNo;
MSEav2=sum(MSE2')/IterNo;
MSEav3=sum(MSE3')/IterNo;
MSEav4=sum(MSE4')/IterNo;
plot(10*log10(MSEav1),'r');hold on
plot(10*log10(MSEav2),'g');hold on
plot(10*log10(MSEav3),'b');hold on
plot(10*log10(MSEav4),'k')
function [fnew]=projection(y,x,f,epsilon)%Projection onto the hyperslab function
x=x';
inner=f'*x;
if((inner-y<-epsilon))
fnew=f+((y-inner-epsilon)/(x'*x))*x;
else if (inner-y>epsilon)
fnew=f+((y-inner+epsilon)/(x'*x))*x;
else
fnew=f;
end
end

Answers (1)

KALYAN ACHARJYA
KALYAN ACHARJYA on 13 Apr 2020
One Suggestion use array or cell array, where it can be, for example-
#Instead of this
MSE1=zeros(N,IterNo);
MSE2=zeros(N,IterNo);
MSE3=zeros(N,IterNo);
MSE4=zeros(N,IterNo);
Use
MSE=cells(1,4)
Instead of this
MSEav1=sum(MSE1')/IterNo;
MSEav2=sum(MSE2')/IterNo;
MSEav3=sum(MSE3')/IterNo;
MSEav4=sum(MSE4')/IterNo;
#DO Array indexing
MSEav(i)=sum(MSE{1}')/IterNo;
See more about vectorization.
  2 Comments
Alice Zurock
Alice Zurock on 13 Apr 2020
Could you replace in my code your advice, because ı confused how can I type correctly?

Sign in to comment.

Categories

Find more on Descriptive Statistics in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!