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 |