white noise case in 3L dft

1 view (last 30 days)
common fernando on 6 Apr 2021
Commented: Walter Roberson on 10 Apr 2021
hello guys, maybe someone of you is able to help me...I tried to fix my problem with the solution above but it wasn`t possible.
Why do I get the following error message:
Subscript indices must either be real positive integers or logicals.
Error in threeDFT (line 33)
x=v((j-1:j+N0-2)*dt);
when I edit code
for j=0:j_max+1
x=v(((50*512):j+N0-2)*dt);
c(j)=x.*Hc';
s(j)=x.*Hs';
I get new error :
Error using *
Inner matrix dimensions must agree.
Error in threeDFT (line 33)
c(j)=x*Hc';
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
and I changed * to .*
Error using .*
Matrix dimensions must agree.
Error in threeDFT (line 33)
c(j)=x.*Hc';
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
what can I do ? I need the solution code for it
when start at J =2 >>>>>>>>>>>>>>>> I get :
Subscript indices must either be real positive integers or logicals.
Error in threeDFT (line 37)
x=v((((j-1):(j+N0-2))*dt));
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
I run this code in other cases not " white noise" and no errors in other cases under the same code parameters
the code according to attached documentation
• The code :
function [t_est,f_est]=threeDFT(v,fs,tmax,N0)
% v : volt as function of time
% fs : sampling frequency (Hz)
% tmax : time of final estimation
% N0 : number of samples in the window
% to test: [t,f]=threeDFT(@(t)(220*sin(2*pi*50.1*t+pi/2)),50*512,1,512)
% to test: [t,f]=threeDFT(@(t)(220*sin(2*pi*50.1*t+pi/2)+randn(size(t))*.1),50*512,1,512)
fs=50*512 ; %sampling freq.
dt =1/fs;
N0=fs/50; %number of samples/cycle
m=3; %no. of cycles
t = dt*(0:m*N0); %data window
fi=50; %Frequency test
ww=wgn(1+N0*m,1,-40);
size(transpose(ww))
y=sin(2*pi*fi*t + 0.3);
v=sin(2*pi*fi*t + 0.3)+transpose(ww);
tmax=1;
n=N0-1:-1:0;
f0=50;
f=50.88;
Hc=2/N0*cos(2*pi*n/N0+pi/N0);
Hs=-2/N0*sin(2*pi*n/N0+pi/N0);
t_est=[];
f_est=[];
j_max=tmax*fs;
for j=2:j_max+1
x=v((((j-1):(j+N0-2))*dt);
c(j)=x*Hc';
s(j)=x*Hs';
if(j>N0)
Ac(j-N0)=sqrt(sum(c(end-N0+1:end).^2)/N0);
As(j-N0)=sqrt(sum(s(end-N0+1:end).^2)/N0);
cc(j-N0)=c(end-N0+1:end)*Hc';
ss(j-N0)=c(end-N0+1:end)*Hs';
if(j>2*N0)
Acc(j-2*N0)=sqrt(sum(cc(end-N0+1:end).^2)/N0);
Ass(j-2*N0)=sqrt(sum(ss(end-N0+1:end).^2)/N0);
ccc(j-2*N0)=cc(end-N0+1:end)*Hc';
ccs(j-2*N0)=cc(end-N0+1:end)*Hs';
ssc(j-2*N0)=ss(end-N0+1:end)*Hc';
sss(j-2*N0)=ss(end-N0+1:end)*Hs';
ff=f0*N0/pi*atan(tan(pi/N0)*((ccc(j-2*N0).^2+ccs(j-2*N0).^2)./(ssc(j-2*N0).^2+sss(j-2*N0).^2)).^.25);
t_est=[t_est;(j-1)*dt];
f_est=[f_est;ff];
end
end
end
t_est;
f_est
Result = bsxfun(@minus, f_est, fi);
R = nansum(Result)
R1 = sum(isfinite(Result))
RS = R/R1
RMSE = sqrt(mean((RS).^2))
plot(t_est, f_est,'red')
hold on
xlabel('time')
ylabel('frequency')
title('three level DFT white noise ')
plot (t,fi)
hold off

Walter Roberson on 6 Apr 2021
for j=1:j_max+1
So j starts at 1
x=v((j-1:j+N0-2)*dt);
j starts at 1, so the first time that is executed it is
x = v((1-1:1+512-2)*1/25600)
which is v((0:511)/25600) . But not one of those subscripts is a positive integer, with the first of them, 0/25600 not even positive.
for j=0:j_max+1
x=v(((50*512):j+N0-2)*dt);
The first time that would be v((50*512):0+512-2)*1/25600) which would be v((25600:510)/25600) . But 25600:510 is an empty row vector, size 1 by 0, so that is v([]) .
When you try to use x (1 by 0) .* anything that is not empty, you have a size mismatch.
When you try to use x (1 by 0) * something, then the inner dimensions must agree -- so the something would have to have 0 rows and some number of columns; if that does somehow turn out to be the case, then the result would be 1 x that many columns, all set to 0.
c(j)=x*Hc'
j starts at at 0, remember, so if the x*Hc' or x.*Hc' worked, you would be trying to assign into c(0) which would fail. If you were to assign to c(j+1) then that would designate a scalar location. With x obviously expected to be a vector, .* cannot produce a scalar result (except when x is accidentally scalar.) With x a vector, you have more opportunity for x*Hc' to result in a scalar -- if Hc was a row vector with the same number of columns as the x had, then x*Hc' would be (1 x something) * (1 x something).' which would be (1 x something) * (something x 1) which would produce a scalar. All it needs is for Hc to increase in size as your v(((50*512):j+N0-2)*dt) increases in size.
Back to v(((50*512):j+N0-2)*dt) : as you increase j large enough, eventually (50*512):j+512-2 would not be empty. But then you have the problem that you multiply (50*512):j+N0-2 by dt, which is 1/25600 and since your j is increasing 1 at a time, you can be sure that the result of the multiplication will not be exact integers at some point, and then you would be in trouble again.
Why are you multiplying indices by dt to try to get new indices?
Walter Roberson on 10 Apr 2021
Our role as volunteers here is more as Tutor than as unpaid consultant. We point out mistakes, and we provide information on how to deal with that category of mistakes... but we not at all likely to go through your paper and write code to match the paper.
Again, I would recommend that you read https://www.mathworks.com/matlabcentral/answers/395816-why-do-i-get-array-indices-must-be-positive-integers-or-logical-values-error-when-using-angle-fu#comment_882605 about the difference between formulas and indexing, and how to convert formulas into indexing. Multiplying indices by dt is a formula for converting indices to time, but indexing by the result is a mismatch, an attempt to index by a formula.