MATLAB: How to get correct values of a(n) from this equation. Attached is the code but it gives wrong values. where n=1…5,

for looploops

clc; close all; clear d=0.5; N=10; for n=1:N/2 sum=0; fderiv = @(f,x) (f(x+eps)-f(x))/eps; % Derivative of a function xderiv = @(x) ((x+eps)-x)/eps; % Derivative of an expression for k = 1:(N/2) x(k)=sin(2*(k+n-1)*pi/N)/(k+n-1) + sin(2*(k-n)*pi/N)/(k-n); if isnan( x(k) ) % Test for ‘NaN’ x(k) = fderiv(@sin, 2*(k-n)*pi/N)/(xderiv(k-n)); % L'Hospital's Rule end end sum=sum+x; an = (2/(N*pi))*(pi-sum); end display(an);% this is not giving the correct values????? an_manually=[-0.5736 -0.4197 -0.1352 0.2383 0.6486] % these are correct values

Best Answer

Here is what I have understood of your problem:

When you calculate a(n), you have a problem in the sum because when k=n the expression sin(2*(k-n)*pi/N)/(k-n) is undefined.

So you try to replace it by its limit when k tends toward n. To evaluate this limit you use L'Hospital's rule.

What I suggest you:

Maybe you should not test for NaN values, but test if k==n, it would be clearer.
As Jan told you, change the name of your sum, something like mySum for example.In the loop on k, you add x to mySum, so define x and not x(k):

...
for n=1:Nele/2
    mySum = 0;
    for k = 1:(Nele/2)
        x = sin(2*(k+n-1)*pi/Nele)/(k+n-1) + sin(2*(k-n)*pi/Nele)/(k-n);
        ...
        mySum = mySum + x;
    end
    an = ...
end
...

Then, when k=n, you forgot the first part of your expression in your code and it is unclear what is your variable. We consider the two functions corresponding to the numerator and the denominator of the expression you want to find the limit of:

f_num = @(k) sin(2*(k-n)*pi/Nele);
f_den = @(k) k-n;

You will use the derivative when k=n, so I think you don't need to evaluate it numerically, because the expression of the derivative taken in n is really simple: it is 2*pi/Nele for the numerator and 1 for the denominator. So your x is:

if k==n
    x = sin(2*(k+n-1)*pi/Nele)/(k+n-1) + 2*pi/Nele;
end

I hope it helps

Related Solutions

MATLAB: Radix 2 disspation in frequency

Zhao - I think the problem may be due to trying to calculate the even and odd terms in your inner for loop (perhaps to improve computation time). If I remember correctly, the DFT (without a shift) can be simplified to just

    X = zeros(N,1);
    for k=0:N-1
        innerSum = 0;
        for n=0:N-1
            innerSum = innerSum + x(n+1)*(cos(2*pi*n*k/N) - 1i*sin(2*pi*n*k/N));
        end
        X(k+1) = innerSum;
    end

where X is your sum_all. But note how as we iterate from 0 to N - 1 in the inner for loop, that the factor x(n+1) is used but use n within the cos and sin function calls. That is not the case in your code where we have factors of 2n or 2n-1 but use just n within the cos and sin. I also don't understand why divide the N by two or the other factor for the sum_odd(k) terms of

 (cos(2*pi*(k-1)/N)-1j*sin(2*pi*(k-1)/N))

Please clarify why these extra steps are being made.

If you do use the above simplified DFT code, then we will observe the frequency at 300 Hz if we do

 f = (0:N-1)*(Fs/N);
 subplot(2,1,2);
 plot(f,abs(X)/(N/2));

which is still very similar to what you have written.

MATLAB: How can i use L’hospital rule i.e sin(x)/x for x tends to zero

By convention, 0/0 evaluates to NaN.

To be absolutely formal (if a bit kludgy) about using L’Hospital's rule in this situation:

fderiv = @(f,x) (f(x+eps)-f(x))/eps;                    % Derivative of a function
xderiv = @(x) ((x+eps)-x)/eps;                          % Derivative of an expression
for n = 0:7
    x(n+1)=sin((pi*n)/2)./(pi*n);
    if isnan( x(n+1) )                                  % Test for ‘NaN’
        x(n+1) = fderiv(@sin, pi*n/2)./(xderiv(pi*n));  % L'Hospital's Rule
    end
end

There are better and more numerically stable derivative expressions (such as the five-point algorithm), but this will do for an illustration.

Also, using L'Hospital's rule,

x=sin((pi*n)/2)./(pi*n)

will be 1 at n=0, however:

x=(sin(pi*n)/2)./(pi*n)

with the fderiv call then becoming:

x(n+1) = (fderiv(@sin, pi*n)/2)./(xderiv(pi*n));

will be 0.5 for n=0.

Best Answer

Related Solutions

MATLAB: Radix 2 disspation in frequency

MATLAB: How can i use L’hospital rule i.e sin(x)/x for x tends to zero

Related Question