What is the Best Approach for Solving a Fourier Coefficient Problem in MATLAB?

[1Keenan]

Hi all,

I have a problem in evaluating Fourier coefficient.
I have an array of numbers (which is a magnetic field) taken on a circle of radius 8mm and I want to know the harmonics of this field.

I have written a code in matlab, I have used some in built function in MATLAB and mathcad, but I can get the correct results.

Any advise?

my code is:

x=Bx(:,4); 

passo=length(x);
period = 2*pi ;
t=linspace(0,period,passo);
t=t';   
Tc=t(2)-t(1); 
T = max(t)-min(t); 
C = 2/T; 
omega = (2*pi/T);



n_armoniche =16;

a0 = C*sum(x)*Tc;



for k=1:n_armoniche
  
    
    
    a(k)=C*sum(x.*cos(omega*k.*t))*Tc; 
    b(k)=C*sum(x.*sin(omega*k.*t))*Tc; 
    c(k)= sqrt(a(k)^2+b(k)^2);
    phi(k)=atand(b(k)/a(k));
      
end

    

xvalues = a0/2*ones(length(t),1);

for kval=1:n_armoniche

        xvalues=xvalues+a(kval)*cos(omega*kval.*t);
        xvalues=xvalues+b(kval)*sin(omega*kval.*t);
    
end

figure()
subplot(3,1,1)
plot(t,x, 'Linewidth', 2)
hold on
plot(t,xvalues,'r','Linewidth', 2)
xlabel('t')
ylabel('function')
axis([min(t) max(t) -inf inf])



subplot(3,1,2)
plot(t, a0/2*ones(1,length(t)), 'Linewidth', 2)
hold all
for kval = 1:n_armoniche
    plot(t,a(kval)*cos(omega*kval.*t), 'Linewidth', 2)
    plot(t,b(kval)*cos(omega*kval.*t), 'Linewidth', 2)
end
xlabel('t')
ylabel('function')
title('Armoniche')
grid on
legend('a0','a1','b1','a2','b2','a3','b3','a4','b4','...')
axis([min(t) max(t) -inf inf])


subplot(3,1,3)
plot(t, a0/2*ones(1,length(t)), 'Linewidth', 2)
hold all
for kval = 1:n_armoniche
    plot(t,c(kval)*cos(omega*kval.*t+phi(kval)), 'Linewidth', 2)
   
end
xlabel('t')
ylabel('function')
title('Armoniche in forma polare')
grid on
legend('c0','c1','c2','c3','c4','c5','c6','c7','c8','...')
axis([min(t) max(t) -inf inf])

 

[Claude Bile]

You need to explain your code better - it is impossible for a third party to understand what you are trying to do...

In my experience, errors in evaluating Fourier coefficients come down to calculating the units in the Fourier domain. Basic rule of thumb is that the domain size (in real space) determines the minimum frequency interval, and that the minimum spatial interval determines the maximum frequency (divided by 2 - because you need to include negative frequencies in your allowed bandwidth).

Claude.

 

[1Keenan]

YOu are right, I'm posting my code, the main part at least, with some comment, hopfully you may help me. I think in my case the error comes down from a bad definition of the period and, I'm not sure, problably it is also related to the fact that I'm analizing a set of data, not a continuous function.

x=Bx(:,4); %Vector of data to be analysed

passo=length(x); %number of element in x
period = 2*pi ; %period of the vector, it is actually evaluated on a circle with fixed radius
t=linspace(0,period,passo);
t=t';
Tc=t(2)-t(1); %sampling intervat
T = max(t)-min(t);
C = 2/T; %constant of integrals
omega = (2*pi/T);
n_armoniche =16; %number of coefficients to be evaluated

a0 = C*sum(x)*Tc; % a0 coefficient
for k=1:n_armoniche
a(k)=C*sum(x.*cos(omega*k.*t))*Tc; %an coefficients
b(k)=C*sum(x.*sin(omega*k.*t))*Tc; %bn coefficients
c(k)= sqrt(a(k)^2+b(k)^2); %modulus of comeplex coefficient
phi(k)=atand(b(k)/a(k)); %phase

end