Graphing Fourier magnitude and phase spectra using fft, error: Array indices must be positive integers or logical values.

I am writing a script to get the frequency spectra of a sawtooth function with a period of 2pi

The problem I am facing is I cannot find a way to use the frequency domain in my code. I tried fft() and fourier() but have been met with errors either way. any tips appreciated!

Goal:

Plotting magitude and phase of Cn and Dn with respect to frequency (see attached images for what I mean)

fft=[];ann=[];cnplot=[];cnt=[];jnt=[];
T = 2*pi;
A=1;

funct1 = @(t) t/T;
funct = @(t) ones(size(funct1));

w0=2*pi*(1/T);

a0 = (1/(T))*integral(funct1,0,T);
for n=1:1:20 % Change the range from 3, 10, 100 
  for t=1:1:360*5
  anfunct= @(t) funct1(t).*cos(n*t*w0);
  bnfunct= @(t) funct1(t).*sin(n*t*w0);
  an = (1/(T/2))*integral(anfunct,0,T);
  bn = (1/(T/2))*integral(bnfunct,0,T);

  cn= (an*an+bn*bn)^(1/2);
  theta=atan((-bn)/an);

  dn= (.5*cn)*theta;

  
  ft(t)=an*cos(n*w0*deg2rad(t))+bn*sin(n*w0*deg2rad(t));
  x(t)=-cn*cos(n*w0*deg2rad(t)+theta);
  d(t)=dn*exp(1i*n*w0*t);

  k(w)=fft(x(t)); %error here: Array indices must be positive integers or logical values.
  end
  fft=[fft;ft];
  cnt=[cnt;x];  
  jnt=[jnt;d];
end

F_A = sum(fft)+a0;
C_A = sum(cnt)+a0;
D_A=sum(jnt)+a0;

tt=deg2rad(1:1:t);

subplot(3,1,1)
plot(tt,F_A)
title("Signal X(t) derived from trigonomtric form")
xlabel("t")
ylabel("f(t)")

subplot(3,1,2)
plot(tt,C_A) 
title("Signal X(t) derived from compact trigonomtric form")
xlabel("t")
ylabel("f(t)")

ww=(1:1:w);
nn=(1:1:n);
%any code below here does not work
subplot(3,1,3);
plot(w,k(w))
title("amplitude spectra")
xlabel("w")
ylabel("k(w)")

Answer 1

[Paul]

Hi Luke,

At the line in question, we see that x(t) is non-integer numeric value, which can't be used to index into the array fft. Also, at that point in the code, w is not defined so that will be a problem as well.

It looks like fft is intended to become a 2D array; are you sure you mean to index into it with a single number?

Finally, consider renaming the variable fft so as not to shadow the built-in function fft.

fft=[];ann=[];cnplot=[];cnt=[];jnt=[];
T = 2*pi;
A=1;
funct1 = @(t) t/T;
funct = @(t) ones(size(funct1));
w0=2*pi*(1/T);
a0 = (1/(T))*integral(funct1,0,T);
for n=1:1:20 % Change the range from 3, 10, 100 
  for t=1:1:360*5
      anfunct= @(t) funct1(t).*cos(n*t*w0);
      bnfunct= @(t) funct1(t).*sin(n*t*w0);
      an = (1/(T/2))*integral(anfunct,0,T);
      bn = (1/(T/2))*integral(bnfunct,0,T);
      cn= (an*an+bn*bn)^(1/2);
      theta=atan((-bn)/an);
      dn= (.5*cn)*theta;
      
      ft(t)=an*cos(n*w0*deg2rad(t))+bn*sin(n*w0*deg2rad(t));
      x(t)=-cn*cos(n*w0*deg2rad(t)+theta);
      d(t)=dn*exp(1i*n*w0*t);
      x(t)
      k(w)=fft(x(t)); %error here: Array indices must be positive integers or logical values.
  end
  fft=[fft;ft];
  cnt=[cnt;x];  
  jnt=[jnt;d];
end

ans = -0.0056

Array indices must be positive integers or logical values.

[Luke McDevitt]

Would a solution be to do something like finding the nearest integer of abs(1000*x(t)) then plotting that in order to find the general shape of the amplitude graph?

[Paul]

I did not look closely, if at all, at what the code is actually trying to compute, so I'm afraid I can't comment on that.