help debug my program.

Prev: nlinfit regression
Next: ode45 - fail to run

From: Gabriel Chan on 14 Oct 2009 10:09

---

Hey all,

I'm currently trying to implement a program from a published journal.
the steps can be found here. http://i34.tinypic.com/6sxkbl.jpg

where y(n) = signal + noise. H is the transpose of the matrix.

AIC and MDL are here http://i35.tinypic.com/24wvfw5.jpg

which brings me to my code.

clear all;
close all;

% Proposed Scheme A
% AIC & MDL

%Setting up the simulation

N = 100; %Number of samples.
Q = 40; %Number of channels

% Step 1 calculate sampling covariance matrix
t = 0;
for n = 1:N
s = [0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 1 1 0 1 1 1 1 0 1 1 1 0 0 0 0 0 0 0 1 0]'; %signal source. 1 means there is a signal in that channel.
e = randn(1,Q)';
Pnoise = e * sqrt(0.1); %10dB SNR
y = s + Pnoise;
G = y*y';
t = t + G;
end

Y = t/N; %Sampling convariance matrix Y = 1/N sum[y(n)*y(n)']

% Step 2 Perform eigenvalue decomposition

F = eig(Y); %F is the eigenvalue decomposition
c = dsort(F); %to get the diagonal values into (p1>p2>...>pq)
f = 1; %initial counter
p = 0;
z = 0;

%Step 3 use AIC or MDL to estimate number of occupied channels
H = 0;
for h = 2:Q
H = H + 1;
z = z + 1;
f = f * c(h) %Product part of AIC equation
p = p + c(h) %Summation part of AIC equation
AIC(z) = -2*N*log10(f/((1/(Q-H)*p)^(Q-H))) + 2*H*(2*Q - H);
MDL(z) = -N*log10(f/((1/(Q-H)*p)^(Q-H))) + 0.5*H*(2*Q-H)*log10(N);

end


[U,K] = min(AIC);
[u,k] = min(MDL);

% Step 4 choose indexes with largest diagonal elements
x = diag(Y);
X = diag(Y);
for b = K:-1:1 % AIC
[r,R] = max(x);
x(R) = 0;

end

for P = k:-1:1 % MDL
[l,L] = max(X);
x(L) = 0;
end

pdcounter = 0;
pd1counter = 0;
for B = 1:Q
if x(B) == 0 && s(B) == 1
pdcounter = pdcounter + 1;
end
end

for Bb = 1:Q
if X(Bb) == 0 && s(Bb) == 1
pd1counter = pd1counter + 1;
end
end


For those that are running the code, you'll find that both k values are 1, which is obviously wrong. can anyone debug it for me? thanks in advance!

From: Gabriel Chan on 15 Oct 2009 00:00

---

up!


"Gabriel Chan" <rafale12(a)hotmail.com> wrote in message <hb4m1v$pj1$1(a)fred.mathworks.com>...
> Hey all,
>
> I'm currently trying to implement a program from a published journal.
> the steps can be found here. http://i34.tinypic.com/6sxkbl.jpg
>
> where y(n) = signal + noise. H is the transpose of the matrix.
>
> AIC and MDL are here http://i35.tinypic.com/24wvfw5.jpg
>
> which brings me to my code.
>
> clear all;
> close all;
>
> % Proposed Scheme A
> % AIC & MDL
>
> %Setting up the simulation
>
> N = 100; %Number of samples.
> Q = 40; %Number of channels
>
> % Step 1 calculate sampling covariance matrix
> t = 0;
> for n = 1:N
> s = [0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 1 1 0 1 1 1 1 0 1 1 1 0 0 0 0 0 0 0 1 0]'; %signal source. 1 means there is a signal in that channel.
> e = randn(1,Q)';
> Pnoise = e * sqrt(0.1); %10dB SNR
> y = s + Pnoise;
> G = y*y';
> t = t + G;
> end
>
> Y = t/N; %Sampling convariance matrix Y = 1/N sum[y(n)*y(n)']
>
> % Step 2 Perform eigenvalue decomposition
>
> F = eig(Y); %F is the eigenvalue decomposition
> c = dsort(F); %to get the diagonal values into (p1>p2>...>pq)
> f = 1; %initial counter
> p = 0;
> z = 0;
>
> %Step 3 use AIC or MDL to estimate number of occupied channels
> H = 0;
> for h = 2:Q
> H = H + 1;
> z = z + 1;
> f = f * c(h) %Product part of AIC equation
> p = p + c(h) %Summation part of AIC equation
> AIC(z) = -2*N*log10(f/((1/(Q-H)*p)^(Q-H))) + 2*H*(2*Q - H);
> MDL(z) = -N*log10(f/((1/(Q-H)*p)^(Q-H))) + 0.5*H*(2*Q-H)*log10(N);
>
> end
>
>
> [U,K] = min(AIC);
> [u,k] = min(MDL);
>
> % Step 4 choose indexes with largest diagonal elements
> x = diag(Y);
> X = diag(Y);
> for b = K:-1:1 % AIC
> [r,R] = max(x);
> x(R) = 0;
>
> end
>
> for P = k:-1:1 % MDL
> [l,L] = max(X);
> x(L) = 0;
> end
>
> pdcounter = 0;
> pd1counter = 0;
> for B = 1:Q
> if x(B) == 0 && s(B) == 1
> pdcounter = pdcounter + 1;
> end
> end
>
> for Bb = 1:Q
> if X(Bb) == 0 && s(Bb) == 1
> pd1counter = pd1counter + 1;
> end
> end
>
>
> For those that are running the code, you'll find that both k values are 1, which is obviously wrong. can anyone debug it for me? thanks in advance!

From: Sebastiaan on 15 Oct 2009 03:48

---

Do you understand the algorithm you are implementing?

If so, make the problem smaller so you can write out the problem by hand. (e.g. N=5, Q=2).

Then use the Matlab debugger to see where the result in your code differs from your solution on paper.

Really, Mathematics is nowhere without pencil and paper.

Sebastiaan


"Gabriel Chan" <rafale12(a)hotmail.com> wrote in message <hb4m1v$pj1$1(a)fred.mathworks.com>...
> Hey all,
>
> I'm currently trying to implement a program from a published journal.
> the steps can be found here. http://i34.tinypic.com/6sxkbl.jpg
>
> where y(n) = signal + noise. H is the transpose of the matrix.
>
> AIC and MDL are here http://i35.tinypic.com/24wvfw5.jpg
>
> which brings me to my code.
>
> clear all;
> close all;
>
> % Proposed Scheme A
> % AIC & MDL
>
> %Setting up the simulation
>
> N = 100; %Number of samples.
> Q = 40; %Number of channels
>
> % Step 1 calculate sampling covariance matrix
> t = 0;
> for n = 1:N
> s = [0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 1 1 0 1 1 1 1 0 1 1 1 0 0 0 0 0 0 0 1 0]'; %signal source. 1 means there is a signal in that channel.
> e = randn(1,Q)';
> Pnoise = e * sqrt(0.1); %10dB SNR
> y = s + Pnoise;
> G = y*y';
> t = t + G;
> end
>
> Y = t/N; %Sampling convariance matrix Y = 1/N sum[y(n)*y(n)']
>
> % Step 2 Perform eigenvalue decomposition
>
> F = eig(Y); %F is the eigenvalue decomposition
> c = dsort(F); %to get the diagonal values into (p1>p2>...>pq)
> f = 1; %initial counter
> p = 0;
> z = 0;
>
> %Step 3 use AIC or MDL to estimate number of occupied channels
> H = 0;
> for h = 2:Q
> H = H + 1;
> z = z + 1;
> f = f * c(h) %Product part of AIC equation
> p = p + c(h) %Summation part of AIC equation
> AIC(z) = -2*N*log10(f/((1/(Q-H)*p)^(Q-H))) + 2*H*(2*Q - H);
> MDL(z) = -N*log10(f/((1/(Q-H)*p)^(Q-H))) + 0.5*H*(2*Q-H)*log10(N);
>
> end
>
>
> [U,K] = min(AIC);
> [u,k] = min(MDL);
>
> % Step 4 choose indexes with largest diagonal elements
> x = diag(Y);
> X = diag(Y);
> for b = K:-1:1 % AIC
> [r,R] = max(x);
> x(R) = 0;
>
> end
>
> for P = k:-1:1 % MDL
> [l,L] = max(X);
> x(L) = 0;
> end
>
> pdcounter = 0;
> pd1counter = 0;
> for B = 1:Q
> if x(B) == 0 && s(B) == 1
> pdcounter = pdcounter + 1;
> end
> end
>
> for Bb = 1:Q
> if X(Bb) == 0 && s(Bb) == 1
> pd1counter = pd1counter + 1;
> end
> end
>
>
> For those that are running the code, you'll find that both k values are 1, which is obviously wrong. can anyone debug it for me? thanks in advance!

 | 
Pages: 1
Prev: nlinfit regression
Next: ode45 - fail to run