Help with Coding problem

I am trying to code the following in matlab.

<<

Following is my matlab code. But this doesn't give me the right result. Can somebody point out what am I doing wrong?

function [J] = funtc(x_t,c_mn,theta,M,N )
L = length(x_t );

l = 0:L-1;
m = 0:M-1; 
n = 0:N-1;

J = 0 ; 
for t = 1:L
    xs = 0 ;
    xp = 0 ;
    for j = 1: length(n)
        for k = 1: length(m)
            xs = xs + c_mn(k,j)*exp(-pi*(l(t)-n(j))^2/theta^2)*exp(1i*2*pi*m(k)*(l(t)- n(j))/M);
            xp = xp + c_mn(k,j)*exp(-pi*(l(t)-n(j))^2/theta^2)*exp(1i*2*pi*m(k)*(l(t)- n(j))/M)*(l(t)-n(j))^2/theta^3;
        end
    end
    J = J + ((x_t (t) -xs) * conj (xp));
end
J = real(-4*J);
end

[Geoff Hayes]

christina - in your code you consider

 c_mn(m,n)

but aren't m and n arrays? Is this intentional? Or do you want to consider

 c_mn(m(k),n(j))

instead?

[christina]

Yes, I want to consider c_mn (k, j).

[Geoff Hayes]

Does that lead to a better (correct) solution?

[christina]

No, I am afraid not. Do you think I have written the code correctly?

[Geoff Hayes]

I would try removing your arrays for n and m as I don't think they make the code easier to follow. Then your two inner loops would become simply

 for n=0:N-1
     for m=0:M-1
         % etc.
     end
 end

Your first sum would (I think) be more like

 for t = 0:L-1
    xs = 0 ;
    xp = 0 ;
    for n = 0:N-1
        for m = 0:M-1
            xs = xs + c_mn(m+1,n+1) * exp(-(t - n)/theta) * exp(1i * 2 * pi * m * (t - n) / M);
            xp = xp + c_mn(m+1,n+1) * exp(-(t - n)/theta) * exp(1i * 2 * pi * m * (t - n) / M) * (((t - n)^2) / theta^3);
        end
    end
    % J = J + ...;
 end

I'm assuming your theta input is the sigma in the equation. The above looks a little different than what you have so I'm not sure if it is better or worse.

I don't understand where the (1/sqrt(theta)) or a comes from (in your code). Or why a conjugate is used in the update to J...

[christina]

I use conj because c_mn can be complex valued. Also please ignore (1/sqrt(theta)) and a. I have removed it from my code above. Am I writing the update of J correctly?

[Geoff Hayes]

But why just use conj on the one sum and not both?

Answer 1

[Abraham Boayue]

Try this

function [J] = funtc(x_t,c,theta,M,N,L)
J = zeros(1,L); 
for t = 1:L
    xs = 0 ;
    xp = 0 ;
    for n = 1: N
        for m = 1: M
            xs = xs + c(n,m)*exp(-pi*((t-n)./theta)).^2*exp(1i*2*pi*m*(t- n)/M);
            xp = xp + c(n,m)*exp(-pi*((t-n)./theta).^2)*exp(1i*2*pi*m*(t- n)/M)*(t-n)^2./theta.^3;
        end
    end
    J = J + ((x_t (t) -xs).*conj (xp));
end
J = -4*real(J);
end

I tested this function using the mfile below.

L = 20;
N = 10;
M = 5;
c =round(10*rand(N,M)); 
theta0 = 0.96;
x0 = 4;
theta = 0:theta0/(L-1):theta0;
x_t = 0:x0/(L-1):x0; 

J =  funtc(x_t,c,theta,M,N ,L);
plot(theta,J);