1 view (last 30 days)

Show older comments

clear all; close all; clc

fs=200; %sampling freq.

dt =1/fs;

N0=fs/50; %number of samples/cycle

m=3; %no. of cycles

t = dt*(0:200); %data window

fi=50; %Frequency test

ww=wgn(201,1,-40);

size(transpose(ww))

X=sin(2*pi*fi*t + 0.3)

v = bsxfun( @plus, X , ww );

tmax=1;

% v : as function of time

% fs : sampling frequency (Hz)

% tmax : time of final estimation

% to test: [t,f]=ZC(@(t)(220*sin(2*pi*50.1*t+pi/2)),50*512,1)

% to test: [t,f]=ZC(@(t)(220*sin(2*pi*50.1*t+pi/2)+randn(1)*.1),50*512,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=1: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

plot(t_est,f_est,'red')

o=rms(fi);

c=rms(f_est)

RMSE = sqrt(mean(c - o).^2)

t_est;

f_est

plot(t_est, f_est,'red')

hold on

RMSE = sqrt(mean((f_est-fi).^2))

xlabel('time')

ylabel('frequency')

title('three LDFT white noise')

plot (t_est,fi*ones(size(t_est)))

hold off

Subscript indices must either be real positive integers or logicals.

Error in threeDFT (line 35)

x=v((j:j+N0-2)*dt);

VBBV
on 2 Dec 2020

for j=2:j_max+1

x=v((j-1:j+N0-2)*dt);

Use a positive non zero integer subscripts. In your case it is zero

Walter Roberson
on 2 Dec 2020

for j=1:j_max+1

x=v((j-1:j+N0-2)*dt);

j starts at 1. j-1 starts at 0. 0 multiplied by anything is 0. Therefore your first index into v would be 0, which is not a permitted index.

If you add one to the indices to avoid this problem, such as starting with j = 2, then j-1:j+N0-2 would be 1:2+4-2 ==> 1:4 which is plausible. Then you multiply that by dt which is 1/fs so dt = 1/200. This would give you an index of (1:4)/200 which would not be integer indices.

Therefore your problem is not just an "off-by-one" problem: you are confusing indexing (relative positions into an array) with the value that the location is intended to denote. You cannot index at 1/fs: you can only index at a value that is associated with 1/fs .

In most cases, your first index into an array should work out to be either 1 or 2 (unless you are reading from a static data array based upon user input or device readings or date)

Walter Roberson
on 4 Dec 2020

When I say reading from a static data array based upon user input or device readings or date, I mean things like:

temp = input('Enter temperature ');

property_idx = floor(temp*10) + 1;

property_value = material_density(property_idx);

or

basetime = datetime('now');

hidx = hms(basetime) + 1;

proj_elec = electrical_projections_by_hour(hidx);

-- cases where the data is not being stepped through in sequence.

Your code, though, is trying to start at the beginning of the data and step through all of the data eventually. And in such a case, your first array index should almost always work out as either 1 or 2.

You need it to be the case that the index of the "next" item in the array is one more (or a constant integer more) than the current index. The index of the "next" item in the array should never be 1/120 more.

Do not confuse the index of an item and the time associated with the index.

timestamp = (0:149) / 120;

y = sin(timestamp * 2 * pi * 5);

Now y(2) is the data value associated with the time given in timestamp(2), which is 1/120. You would not index the array y with y(1/120) to get the y value associated with the time 1/120: you would index y(2) to get the y value associated with the time given in timestamp(2) where timestamp(2) = 1/120 .

Do not confuse array indexing with formulas. You might have a formula

Y = @(t) sin(t * 2 * pi * 5)

such that Y(1/120) would compute the result associated with the value 1/120, but a formula is not an array, where the numbers in () tell you about relative location of the information.

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

Start Hunting!