欠定盲源分离经典程序

function  [y  A]= sparseBSS1（XLlangdaGhdelta）
%----------------------------------------------------------------
%                  2009-04-15 YangZhicong
% X: observed signaleach 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 casei.e m = 2
[m T] = size（X）;
figure（‘Name‘‘guance‘）;
plot（X（1）X（2）‘.‘）;
if nargin < 2L = 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（Xwoverlap）; % decompose X into frames. For instance frame_X（::1）
% correspond to 1st sensoreach 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];
    [rowcolumn] = size（temp_Matrix）;
    Xu（i:） = reshape（temp_Matrix 1row*column）;
end

if nargin < 4 G =60; end
if nargin < 3langda = 5; end 
% we consider the special case for m = 2
theta = rem（atan2（Xu（2:）Xu（1:））+2*pipi）;
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 < 5h = 0.2;end % discarding the less reliable data points
index = find（radius>=h）;
radius = radius（index）;
theta = theta（index）;

figure（‘Name‘‘Scatter Plot ‘）;
polar（thetaradius‘*‘）;

for k = 1:G
    potential（k） = sum（radius.*Tao（langda*（theta-theta_k（k））））;
end
figure（‘Name‘‘Potential Funtction‘）;
plot（theta_k/pi*180potential）;

%        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;
figure（‘Name‘‘source singal‘）;
plot（y（1）y（2）‘.‘）;
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.
%        MAX