How to find the RMSE?

I have the following code:

clear;clc
ula = phased.ULA('NumElements',10,'ElementSpacing',0.5);
angs = [40 -20 20; 0 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2); 
c = physconst('LightSpeed');
fc = 300e6;  % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator('SensorArray',ula,...
        'OperatingFrequency',fc,'ScanAngles',-90:90,...
        'DOAOutputPort',true,'NumSignalsSource','Property','NumSignals',NumSignals);
[~,ang] = musicspatialspect(signal)

ang = 1×3

40 -20 20

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ymusic = musicspatialspect(signal);
helperPlotDOASpectra(musicspatialspect.ScanAngles,ymusic)% Changed by Me

function helperPlotDOASpectra(x2,y2)% Changed by Me
  % Plot spectra in dB normalized to 0 dB at the peak 
  y2_dB = 20*log10(y2) - max(20*log10(y2));
  plot(x2,y2_dB)
  xlabel('Broadside Angle (degrees)');
  ylabel('Power (dB)');
  title('DOA Spatial Spectra')
 end

This estimates the angles and displays them in the command window. It also gives me the plot. Now I want to find the RMSE plot i.e., the root mean square error vs number of runs plot , where number of runs=100. The RMSE is between the actual angles i.e., 40 -20 20 assigend to "angs" and the estimated angles in variable "ang". But I don't know how to do it?

Answer 1

[aditi bagora]

Hi @Sadiq Akbar,

To calculate the RMSE between the values of 'angs' and 'ang' over multiple runs, you can store the RMSE values in an array for each run and then you can plot these errors against the run number using the 'plot' function

Please find the sample code attached below:

num_runs = 100;
rmse_values = zeros(1, num_runs);
for idx = 1:num_runs    
    error = angs - ang;
    rmse = sqrt(mean(error.^2));    
    rmse_values(idx) = rmse;
end
figure;
plot(1:num_runs, rmse_values, '-o');
xlabel('Run Number');
ylabel('RMSE');
title('RMSE vs Number of Runs');
grid on;

Hope this helps!

[Sadiq Akbar]

@aditi bagora

Thanks a lot for your kind response. Yes, I know that. But I don't know how to run the MUSIC algorithm 100 times as each time it gives the same values of estimated angles as are the desired one.

[Torsten]

Then the RMSE is

zeros(1,100)

ans = 1×100

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Sadiq Akbar]

@Torsten

Thanks a lot for your kind response. No this is not correct because I have seen RMS plots in papers but here all values are zeros which is not correct.

[Torsten]

But if you make no error in the angle estimation, RMSE is 0, isn't it ?

[Sadiq Akbar]

@Torsten

Sorry for being too late in response. Yes you are right but if we add noise in the signal say for example -15dB, then the RMSE should not all be zeros but should be some very very small values though non-zero like 10^-5, 10^-4, 10^-3 and so on. For this consider the following code:

clear;clc
ula = phased.ULA('NumElements',4,'ElementSpacing',0.5);
angs = [40 -20; 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2); 
c = physconst('LightSpeed');
fc = 300e6;  % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 3;%1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
signal1 = awgn(signal,-15); % Add -15dB Noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator('SensorArray',ula,...
        'OperatingFrequency',fc,'ScanAngles',-90:90,...
        'DOAOutputPort',true,'NumSignalsSource','Property','NumSignals',NumSignals);
[~,ang] = musicspatialspect(signal)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RMSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num_runs = 100;
rmse_values = zeros(1, num_runs);
for idx = 1:num_runs    
    error = angs - ang;
    rmse = sqrt(mean(error.^2));    
    rmse_values(idx) = rmse;
end
figure;
plot(1:num_runs, rmse_values, '-o');
xlabel('Run Number');
ylabel('RMSE');
title('RMSE vs Number of Runs');
grid on;

But this gives error.

[Torsten]

I don't understand what you are doing.

Technically, "angs" is 2x2, "ang" is 1x2 and thus "rmse" is 1x2. You can't save a vector in a scalar component, thus

rmse_values(idx) = rmse;

throws an error.

But the complete loop

for idx = 1:num_runs    
    error = angs - ang;
    rmse = sqrt(mean(error.^2));    
    rmse_values(idx) = rmse;
end

doesn't make sense because you compute the same thing "num_runs" times (you only have data for num_runs = 1 times).

And don't call a variable "rmse" - you overwrite the MATLAB function with the same name:

https://uk.mathworks.com/help/matlab/ref/rmse.html

Answer 2

[Steven Lord]

If you're using release R2022b or later, use the rmse function on your data.

[Sadiq Akbar]

@Torsten

Thanks a lot for your response. I modified the code as below:

clear;clc
num_runs = 1000;
rmse_values = zeros(1, num_runs);
for idx = 1:num_runs

ula = phased.ULA('NumElements',4,'ElementSpacing',0.5);
angs = [40 -20; 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2); 
c = physconst('LightSpeed');
fc = 300e6;  % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 3;%1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
signal1 = awgn(signal,-15);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator('SensorArray',ula,...
        'OperatingFrequency',fc,'ScanAngles',-90:90,...
        'DOAOutputPort',true,'NumSignalsSource','Property','NumSignals',NumSignals);
[~,ang] = musicspatialspect(signal)


%%%%%%%%%%%%%%%%%%%%
% Swapping vector b % Don't activate swapping for GA as it makes result wrong
%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(angs(1,:)); % u is my desired vector
[~, ix1(ix)] = sort(ang);  % temp stores the randomly generated vector "best" by algorithm
ang = ang(ix1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RMSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    error(idx,:) = angs(1,:) - ang;
    rmse(idx) = sqrt(meanabs((error.^2)));    
    rmse_values(idx) = rmse(idx);
end
figure;
plot(1:num_runs, rmse_values, '-o');
xlabel('Run Number');
ylabel('RMSE');
title('RMSE vs Number of Runs');
grid on;

But now it gives all values of MSE as zeros which shouldn't be the case because now we have added the noise.

[Torsten]

Shouldn't it be

[~,ang] = musicspatialspect(signal1)

instead of

[~,ang] = musicspatialspect(signal)

?

And again: Don't name a variable "rmse" as you do - you overwrite the MATLAB function "rmse".

[Sadiq Akbar]

@Torsten

Thanks a lot for your guiadance. Yes you are right. It should be:

[~,ang] = musicspatialspect(signal1)

instead of

[~,ang] = musicspatialspect(signal)

but still the values of RMSE are not like 10^-30, 10^-31, 10^-32,..... but instead all the values are "7.3824" which is wrong.

[Torsten]

I cannot help you further because background knowledge about your application would be necessary to give advice.