科研心酸历程

first

clear;
clc;
fs = 4000;
sigduan = 10;
N = 100;%% set the length of original signal 
whaledata = wavread('11.wav');
y = whaledata(1:N);
%N = length(y);
n = N;
y = y-mean(y);
y = y/max(abs(y));
dictionary = dctmtx(N+20);
%%% analysis frequency
sigfft = abs(fft(y,n));
f = 0:n-1;
f = f*fs/N;
sigfft(1) = sigfft(1)/2;
figure(1);
plot(f,sigfft*2/N);
title('distribution of original frequency','fontsize',20);
xlabel('frequency/Hz','fontsize',18);
ylabel('amplitude/v','fontsize',18);
%%% deciding which scale should the signal use
wavfc = centfrq('sym3');
wscal1 = fs*wavfc/max(f(1:N/2));
wscal2 = ceil(wscal1);
if(wscal2<ceil(log2(sigduan)))
    wscal3 = ceil(log2(sigduan));
else wscal3 = wscal2;
end
%%% using wavelet decomposition to find the position of weak energy
sigtree = wpdec(y,wscal3,'sym3');
plot(sigtree);
for i = 1:2^wscal3
    sigE(i) = norm(wpcoef(sigtree,[wscal3,i-1]),2);%energy in node of [wscal3,i]
end
sigsumE = sum(sigE);%totel energy
for i = 1:2^wscal3
    sigpfir(i) = sigE(i)/sigsumE;%the probability of every node 
end
%%% deviding original signal
sigord = sort(sigpfir);
for i = 1:sigduan
    [row,column] = find(sigpfir==sigord(i));
    sigrow(i) = row;
    sigcolumn(i) = column;
end
wavef = scal2frq(1:2^wscal3,'sym3',1/fs);
sigidx = wavef(sigrow,sigcolumn);
sigidx = sigidx(1,:);
sigduanidx = sort(sigidx);
%%%sparse coding by subsection
res = [];
sols = [];
iters = [];
activ = [];
tempfd = 1;
j = 1;
done = 1;
tic,
%for i = 1:sigduan
for i = 1:sigduan
   while done % j = 1:N
       if(sigduanidx(i)<f(j))
            sigfd = y(tempfd:j);
            tempfd = j;
            j = j+1;
            done = 0;
       else 
           j = j+1;
       end

   end
    % erro:which place does this end in? that is nothing impotance
    len = length(sigfd);
    fmax = max(f(1:N/2));
    dicstandcol = 1/(2-4*fmax/fs);
    tempdic = num2str(dicstandcol);
    dicstandcol = str2num(tempdic);
    if (len<dicstandcol)
        A1 = dictionary(1:len,1:dicstandcol);  %(1:len,j:sigduanidx(i));%erro :extract atoms
    else 
        A1 = dictionary(1:len,:);
    end
    [temprow,tempcol] = size(A1);
    [res1,sols1,iters1,activ1] = solveOMP(A1,sigfd,tempcol);%%%the 163rd line
    res = [res;res1];
    sols = [sols;sols1];
    iters = [iters;iters1];
    activ = [activ;activ1];
    %%j = sigduanidx(i);   %after extract dictionary
    i = i+1;   
end
fprintf('calculate time of subsection = %f\n',toc);
%%% sparse coding by whole original signal 
tic,
A2 = dictionary(1:N,:);
[resw,solsw,itersw,activw] = solveOMP(A2,y,N);
fprintf('calculate time of whole signal = %f\n',toc)

%%%comparison between subsection and whole signal
MSE_sub = sum(abs(res))/N;
MSE_whole = sum(abs(resw))/N;
xl = 0:N*1000/fs;
figure(1);
set(1,'Name','res_comparison');
plot(xl,MSE_sub,'k--',xl,MSE_whole,'r-');
title([' MSE of subsection = %f,MSE of whole signal = %f',MSE_sub,MSE_whole]);
xlabel('time/ms');
ylabel('amplitude');
set(gca,'xtick',1000/fs);
legend('MSEsub','MSEwhole');

sparksub = 0;
sparkwhole = 0;
for i = 1:N
    if (sols(i)<0.0001)
        sparksub = sparksub+1;
    if(solsw(i)<0.0001)
        sparkwhole = sparkwhole+1;
    end
    end
end

figure(2);
set(2,'Name','comparison of sols ');
subplot(211);
title(['spark of sub = %d',sparksub]);
stem(xl,sols);
xlabel('time/ms');
ylable('the value of coefficients');
subplot(212);
title(['spark of whole = %d',sparkwhole]);
stem(xl,solsw);
xlabel('time/ms');
ylable('the value of coefficients');
set(gca,'xtick',1000/fs);

so many bugs make me angry

second

clear;
clc;
fs = 4000;
sigduan = 10;
N = 100;%% set the length of original signal 
whaledata = wavread('11.wav');
y = whaledata(1:N);
%N = length(y);
n = N;
y = y-mean(y);
y = y/max(abs(y));
dictionary = dctmtx(N+20);
%%% analysis frequency
sigfft = abs(fft(y,n));
f = 0:n-1;
f = f*fs/N;
sigfft(1) = sigfft(1)/2;
figure(1);
plot(f,sigfft*2/N);
title('distribution of original frequency','fontsize',20);
xlabel('frequency/Hz','fontsize',18);
ylabel('amplitude/v','fontsize',18);
%%% deciding which scale should the signal use
wavfc = centfrq('sym3');
wscal1 = fs*wavfc/max(f(1:N/2));
wscal2 = ceil(wscal1);
if(wscal2<ceil(log2(sigduan)))
    wscal3 = ceil(log2(sigduan));
else wscal3 = wscal2;
end
%%% using wavelet decomposition to find the position of weak energy
sigtree = wpdec(y,wscal3,'sym3');
plot(sigtree);
for i = 1:2^wscal3
    sigE(i) = norm(wpcoef(sigtree,[wscal3,i-1]),2);%energy in node of [wscal3,i]
end
sigsumE = sum(sigE);%totel energy
for i = 1:2^wscal3
    sigpfir(i) = sigE(i)/sigsumE;%the probability of every node 
end
%%% deviding original signal
sigord = sort(sigpfir);
for i = 1:sigduan
    [row,column] = find(sigpfir==sigord(i));
    sigrow(i) = row;
    sigcolumn(i) = column;
end
wavef = scal2frq(1:2^wscal3,'sym3',1/fs);
sigidx = wavef(sigrow,sigcolumn);
sigidx = sigidx(1,:);
sigduanidx = sort(sigidx);
%%%sparse coding by subsection
res = [];
sols = [];
iters = [];
activ = [];
tempfd = 1;
j = 1;
done = 1;
tic,
%for i = 1:sigduan
for i = 1:sigduan
   while done % j = 1:N
       if(sigduanidx(i)<=f(j)/2&&i<sigduan)
            sigfd = y(tempfd:j);
            tempfd = j;
            j = j+1;
            done = 0;
       else if(i==sigduan&&j<N)  
               sigfd = y(j:N);
           else
           j = j+1;
           end
       end


   end
    % erro:which place does this end in? that is nothing impotance
    len = length(sigfd);
    fmax = max(f(1:N/2));
    dicstandcol = 1/(2-4*fmax/fs);
    tempdic = num2str(dicstandcol);
    dicstandcol = str2num(tempdic);
    if (len<dicstandcol)
        A1 = dictionary(1:len,1:dicstandcol);  %(1:len,j:sigduanidx(i));%erro :extract atoms
    else 
        A1 = dictionary(1:len,:);
    end
    [temprow,tempcol] = size(A1);
    [res1,sols1,iters1,activ1] = solveOMP(A1,sigfd,tempcol);%%%the 163rd line
    res = [res;res1];
    sols = [sols;sols1];
    iters = [iters;iters1];
    %activ = [activ;activ1];
    %%j = sigduanidx(i);   %after extract dictionary
    i = i+1; 
    done = 1;
end
fprintf('calculate time of subsection = %f\n',toc);
%%% sparse coding by whole original signal 
tic,
A2 = dictionary(1:N,:);
[resw,solsw,itersw,activw] = solveOMP(A2,y,N+20);
fprintf('calculate time of whole signal = %f\n',toc)

%%%comparison between subsection and whole signal
MSE_sub = sum(abs(res))/length(res);
MSE_whole = sum(abs(resw))/N;
xl = (1*1000/fs):(1*1000/fs):(N*1000/fs);
xls = (1*1000/fs):(1*1000/fs):(length(res)*1000/fs);
figure(3);
set(1,'Name','res_comparison');
plot(xls,res,'k--',xl,resw,'r-');
title([' MSE of subsection = ',num2str(MSE_sub),'MSE of whole signal = %f',num2str(MSE_whole)]);
xlabel('time/ms');
ylabel('amplitude');
set(gca,'xtick',1000/fs);
legend('MSEsub','MSEwhole');

sparksub = 0;
sparkwhole = 0;
for i = 1:N
    if (sols(i)<0.0001)
        sparksub = sparksub+1;
    if(solsw(i)<0.0001)
        sparkwhole = sparkwhole+1;
    end
    end
end

figure(4);
set(2,'Name','comparison of sols ');
subplot(211);
title(['spark of sub = ',num2str(sparksub)]);
stem(sols);
xlabel('time/ms');
ylabel('the value of coefficients');
subplot(212);
title(['spark of whole = ',num2str(sparkwhole)]);
stem(solsw);
xlabel('time/ms');
ylabel('the value of coefficients');
set(gca,'xtick',1000/fs);

bug is still exist,so……

third

clear;
clc;
%%%%%use original signal to divide energy rather than signal after wavedec
fs = 4000;
sigduan = 10;
N = 100;%% set the length of original signal 
whaledata = wavread('11.wav');
y = whaledata(1:N);
%N = length(y);
n = N;
y = y-mean(y);
y = y/max(abs(y));
dictionary = dctmtx(N+20);
%%% analysis frequency
sigfft = abs(fft(y,n));
f = 0:n-1;
f = f*fs/N;
sigfft(1) = sigfft(1)/2;
figure(1);
plot(f,sigfft*2/N);
title('distribution of original frequency','fontsize',20);
xlabel('frequency/Hz','fontsize',18);
ylabel('amplitude/v','fontsize',18);
%%% deciding which scale should the signal use
E = abs(y).^2;
Esort = sort(E);
divis = sigduan-1;
for i = 1:divis
    [row,column] = find(E==Esort(i));
    if (size(row)==1)
    sigrow(i) = row;
  %  sigcolumn(i) = column;
    else
        sigrow(i) = row(1);
        i = i+2;
     %   sigcolumn(i) = column(1);
        divis = divis+2;
    end

end
sigrow = sort(sigrow);
%sigcolumn = sort(sigcolumn);
%%% using wavelet decomposition to find the position of weak energy
%%% deviding original signal
%%%sparse coding by subsection
res = [];
sols = [];
iters = [];
activ = [];
tempr = 1;
%tempc = 1;
j = 1;
tic,
%for i = 1:sigduan
for i = 1:sigduan-1
    if(i<sigduan)
   sigfd = y(tempr:sigrow(i),1);
   if(sigrow(i)>1)
       tempr = sigrow(i)+1;
   end
    else
        sigfd = y(sigrow(i-1):N);
    end
%    if(sigcolumn(i)>1)
%        tempc = sigcolumn(i)+1;
%    end

    % erro:which place does this end in? that is nothing impotance
    len = length(sigfd);
    fmax = max(f(1:N/2));
    dicstandcol = 1/(2-4*fmax/fs);
    tempdic = num2str(dicstandcol);
    dicstandcol = str2num(tempdic);
    if (len<dicstandcol)
        A1 = dictionary(1:len,1:dicstandcol);  %(1:len,j:sigduanidx(i));%erro :extract atoms
    else 
        A1 = dictionary(1:len,:);
    end
    [temprow,tempcol] = size(A1);
    %[res1,sols1,iters1,activ1] = solveOMP(A1,sigfd,tempcol);%%%the 163rd line
    [res1, sols1] = SolveBP(A1,sigfd,tempcol);%, maxIters, lambda, OptTol)
    res = [res;res1];
    sols = [sols;sols1];
    %iters = [iters;iters1];
    %activ = [activ;activ1];
    %%j = sigduanidx(i);   %after extract dictionary
    i = i+1; 
end
fprintf('calculate time of subsection = %f\n',toc);
%%% sparse coding by whole original signal 
tic,
A2 = dictionary(1:N,:);
[resw,solsw] = solveBP(A2,y,N+20);
fprintf('calculate time of whole signal = %f\n',toc)

%%%comparison between subsection and whole signal
MSE_sub = sum(abs(res))/length(res);
MSE_whole = sum(abs(resw))/N;
xl = (1*1000/fs):(1*1000/fs):(N*1000/fs);
xls = (1*1000/fs):(1*1000/fs):(length(res)*1000/fs);
figure(3);
set(3,'Name','res_comparison');
plot(xls,res,'k--',xl,resw,'r-');
title([' MSE of subsection = ',num2str(MSE_sub),'     MSE of whole signal = ',num2str(MSE_whole)]);
xlabel('time/ms');
ylabel('amplitude');
set(gca,'xtick',1000/fs);
legend('MSEsub','MSEwhole');

sparksub = 0;
sparkwhole = 0;
for i = 1:N
    if (sols(i)<0.0001)
        sparksub = sparksub+1;
    if(solsw(i)<0.0001)
        sparkwhole = sparkwhole+1;
    end
    end
end

figure(4);
set(4,'Name','comparison of sols ');
subplot(211);
stem(sols);
title(['spark of sub = ',num2str(sparksub)]);
xlabel('time/ms');
ylabel('the value of coefficients');
subplot(212);
stem(solsw);
title(['spark of whole = ',num2str(sparkwhole)]);
xlabel('time/ms');
ylabel('the value of coefficients');

results were getten ,but it is not what I want to get,in other words ,it is not ideal