sparsemang1.zip_欠定_欠定盲_欠定盲分离_盲源分离_稀疏 盲

function [y,A]= sparseBSS1(X,L,langda,G,h,delta) %---------------------------------------------------------------- % 2009-04-15 YangZhicong % X: observed signal,each row correspond to a sensor observations % L: the length of FFT (or the length of Hanning window) % langda: adjust the desired angular width % G: % discretize the potential field by taking a sample of G points % h: threshold % y: estimated sourses % A: estimated mix matrix % here we just consider the special case,i.e, m = 2 [m, T] = size(X); if nargin < 2,L = 2048; end d = round(0.15*L)*2; % the hop distance overlap = L - d; % number of samples of overlap between adjacent windows w = hann(L)'; % Hann (Hanning) window function; Frame_X = bss_make_frames(X,w,overlap); % decompose X into frames. For instance, Frame_X(:,:,1) % correspond to 1st sensor,each row is a frame. % now, each column of Frame_X(:,:,i) is a frame, for convenience of fft operation for i=1:m Frame_X_temp(:,:,i) = Frame_X(:,:,i)'; end Frame_X = Frame_X_temp; clear Frame_X_temp; Frame_X_Fs = fft(Frame_X); K = L/2 +1; % preserve the positive half spectrum for i = 1:m temp_Matrix = Frame_X_Fs(1:K,:,i); realPart = real(temp_Matrix); imagPart = imag(temp_Matrix); temp_Matrix = [realPart;imagPart]; [row,column] = size(temp_Matrix); Xu(i,:) = reshape(temp_Matrix ,1,row*column); end if nargin < 4, G =60; end if nargin < 3,langda = 5; end % we consider the special case for m = 2 theta = rem(atan2(Xu(2,:),Xu(1,:))+2*pi,pi); radius = sqrt((Xu(1,:).*Xu(1,:)+Xu(2,:).*Xu(2,:))); radius = radius /max(radius); % normalize to one theta_k = pi/2/G + (1:G)*pi/G; if nargin < 5,h = 0.2;end % discarding the less reliable data points index = find(radius>=h); radius = radius(index); theta = theta(index); figure('Name','Scatter Plot '); polar(theta,radius,'.'); for k = 1:G potential(k) = sum(radius.*Tao(langda*(theta-theta_k(k)))); end figure('Name','Potential Funtction'); plot(theta_k/pi*180,potential); % A point is considered a maximum peak if it has the maximal % value, and was preceded (to the left) by a value lower by % DELTA. if nargin < 6, delta = max(abs(potential))*0.1; end maxtab = peakdet(potential, delta, theta_k); % find the maxima of potential function theta_esti = maxtab(:,1); % get the first m maxima of potential function theta_esti = sort(theta_esti,'descend'); theta_esti = theta_esti (1:m); A = [cos(theta_esti) sin(theta_esti)]'; W = inv(A); % separation matrix y = W*X; end %% % ------------------------------PEAKDET.m-------------------------------- function [maxtab, mintab]=peakdet(v, delta, x) %PEAKDET Detect peaks in a vector % [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local % maxima and minima ("peaks") in the vector V. % MAXTAB and MINTAB consists of two columns. Column 1 % contains indices in V, and column 2 the found values. % % With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices % in MAXTAB and MINTAB are replaced with the corresponding % X-values. % % A point is considered a maximum peak if it has the maximal % value, and was preceded (to the left) by a value lower by % DELTA. %%%%%%%% http://www.billauer.co.il/peakdet.html % Eli Billauer, 3.4.05 (Explicitly not copyrighted). % This function is released to the public domain; Any use is allowed. maxtab = []; mintab = []; v = v(:); % Just in case this wasn't a proper vector if nargin < 3 x = (1:length(v))'; else x = x(:); if length(v)~= length(x) error('Input vectors v and x must have same length'); end end if (length(delta(:)))>1 error('Input argument DELTA must be a scalar'); end if delta <= 0 error('Input argument DELTA must be positive'); end mn = Inf; mx = -Inf; mnpos = NaN; mxpos = NaN; lookformax = 1; for i=1:length(v) this = v(i); if this > mx, mx = this; mxpos = x(i); end if this < mn, mn = this; mnpos = x(i); end if lookformax if this < mx-delta maxtab = [maxtab ; mxpos mx]; mn = this; mnpos = x(i); lookformax = 0; end else if this > mn+delta mintab = [mintab ; mnpos mn]; mx = this; mxpos = x(i); lookformax = 1; end end end end %% % --------------------------------Tao.m-------------------------------- function [tao] = Tao(a) tao = 1 -abs(a) /(pi/4); tao(tao>1 | tao<0)=0; end