For Loop - Getting error - Index in position 2 exceeds array bounds

I am trying to run this loop for snr and then plot snr wrt pmusic.

However, I am getting this array indexing error.

I have tried to remove ii from some variables inside the loop, re-inserted, placed differently, but the errors keep on coming.

I dont know how to solve this.

Getting this Error:

Index in position 2 exceeds array bounds (must not exceed 1).

Error in MUSIC_RMSE_wrt_SNR_08may22 (line 39)

NN(ii) = N(:,1:M-P()); %Estimate noise subspace

D=zeros(P,M); %To create a matrix with P row and M column
for k=1:P
D(k,:)=exp(-1j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %Assignment matrix
end
D=D';

xx=2*exp(j*(w*[1:N])); %Simulate signal
x=D*xx;

theta=30;
snr = 0:0.5:40;
for ii=1:numel(snr)
    x1(ii)=x(ii)+awgn(x(ii),snr(ii));
    R(ii) = x1(ii)*x1(ii)' %Data covarivance matrix
    [N,V] = eig(R(ii)); %Find the eigenvalues and eigenvectors of R
    NN(ii) = N(:,1:M-P); %Estimate noise subspace
    SS=zeros(1,length(M));
    for jj=0:M-1
        SS(1+jj)=exp(-j*2*jj*pi*d*sin(theta/180*pi)/lambda);
    end
    PP=SS*NN*NN'*SS';
    Pmusic(ii)=abs(1/ PP);
end

Pmusic=10*log10(Pmusic/max(Pmusic))%Spatial spectrum function

plot(snr,Pmusic,'-k')

[Walter Roberson]

    R(ii) = x1(ii)*x1(ii)' %Data covarivance matrix
    [N,V] = eig(R(ii));

what is the eigenvectors of a scalar?

[Umer Khalid]

Thanks Walter.

Which scalar are you referring to please?

[Umer Khalid]

I have successfully implemented this code wrt theta.

Actually now I am trying to modify this code to run and plot it wrt SNR.

The code wrt theta is given below:

xx=2*exp(j*(w*[1:N])); %Simulate signal
x=D*xx;

x1=x+awgn(x,snr);
R=x1*x1'; %Data covarivance matrix
[N,V]=eig(R); %Find the eigenvalues and eigenvectors of R
NN=N(:,1:M-P); %Estimate noise subspace
theta=-90:0.5:90; %Peak search
for ii=1:length(theta)
    SS=zeros(1,length(M));
    for jj=0:M-1
        SS(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);
    end
    PP=SS*NN*NN'*SS';
    Pmusic(ii)=abs(1/ PP);
end
Pmusic=10*log10(Pmusic/max(Pmusic))%Spatial spectrum function

plot(theta,Pmusic,'-k')

[Walter Roberson]

What would you want the outputs to be in that case?

[Umer Khalid]

I want to find RMS Error after that wrt SNR.

I have to run Monte Carlo simulations and then calculate RMSE wrt SNR.

[KSSV]

Many variables are not defined in the given code. Show us the full code, so that we can help you.

[Umer Khalid]

Well this is the full code.

I have implemented MUSIC Algo for finding RMSE in direction finding.

Now the code works like this, but I am not sure whether I am doing it right, especially the Monte Carlo Iterations part.

%% MUSIC RMSE for Wideband Signals versus SNR 
%% Fixed DOA - Theta angle
%% Uniform Circular array - 4 elements antenna

clc 
clear
close all
%%  design the UCA
M = 4; % 4 antenna elements
N = 250; % Number of Samples - Keeping at 250 now - The higher the better - But depends on receiving system
Rad = 250e-3; % Radius of circular array = 250mm
f = 225e6; %Frequency of operation - 225 MHz
c = 3e8; % Speed of light - constant
lambda = c/f; % Wavelength
d = lambda/4.28; % Inter Element spacing equals 0.312mm in diagonals


w=[pi/6]';%Frequency
doadeg = [30]
doa = doadeg/180*pi;
P=length(w); %The number of signal

snr1=-20;
snr2=0;
snr3=20;
D=zeros(P,M);
for k=1:P
    D(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]);
end

D=D';
xx=2*exp(j*(w*[1:N]));
x=D*xx;

x1=x+awgn(x,snr1);
x2=x+awgn(x,snr2);
x3=x+awgn(x,snr3);

R1=x1*x1';
R2=x2*x2';
R3=x3*x3';

[N1,V1]=eig(R1);
[N2,V2]=eig(R2);
[N3,V3]=eig(R3);

NN1=N1(:,1:M-P);
NN2=N2(:,1:M-P);
NN3=N3(:,1:M-P);

theta=-90:0.5:90;

Pmusic1=zeros(1,length(theta));
Pmusic2=zeros(1,length(theta));
Pmusic3=zeros(1,length(theta));

for ii=1:length(theta)
    SS=zeros(1,length(M));
    for jj=0:M-1
    SS(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);
    end
    
    PP1=SS*NN1*NN1'*SS';
    PP2=SS*NN2*NN2'*SS';
    PP3=SS*NN3*NN3'*SS';
    
    Pmusic1(ii)=abs(1/ PP1);
    Pmusic2(ii)=abs(1/ PP2);
    Pmusic3(ii)=abs(1/ PP3);
end

Pmusic1=10*log10(Pmusic1/max(Pmusic1));
Pmusic2=10*log10(Pmusic2/max(Pmusic2));
Pmusic3=10*log10(Pmusic3/max(Pmusic3));


plot(theta,Pmusic1,'--k','linewidth',2.0)
hold on
plot(theta,Pmusic2,'-k','linewidth',1.0)
hold on
plot(theta,Pmusic3,'-.k','linewidth',0.1)
hold off

xlabel('angle \theta/degree')
ylabel('spectrum function P(\theta) /dB')
title('DOA estimation based on MUSIC algorithm')
grid on

[pks,locs] = findpeaks(Pmusic1, theta); % TO find peak points at their locations
RMSE1 = sqrt(mean((locs - doadeg).^2))



%% Test Monte carlo for 1000 iterations
N = 1000;
for n = 1:N
   y = locs*rand(1,N)
   yest = doadeg*rand(1,N)
   RMSE5(n) = sqrt(mean((y - yest).^2))
end

cumRMSE5 = cumsum(RMSE5) ./ [1:N];
meanRMSE5=mean(cumRMSE5)
figure(3)
plot(cumRMSE5)

%% Finding Error (RMSE)

[pks,locs] = findpeaks(Pmusic1, theta); % TO find peak points at their locations
RMSE1 = sqrt(mean((locs - doadeg).^2))


[pks,locs] = findpeaks(Pmusic2, theta); % TO find peak points at their locations
RMSE2 = sqrt(mean((locs - doadeg).^2))

[pks,locs] = findpeaks(Pmusic3, theta); % TO find peak points at their locations
RMSE3 = sqrt(mean((locs - doadeg).^2))




x=[snr1 snr2 snr3]
y=[RMSE1 RMSE2 RMSE3];
figure(2)
plot(x,y)
grid on