bdist.m

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function [ddb] = bdist(x,y,fs)
% bdist - Evaluates the distortion introduced by the
%         Blind Signal Separation algorithm
%
% usage: ddb = bdist(x,y,fs);
%
% x(:,i,j) - The i^th microphone signal with only the j^th source active
% y(:,i,j) - The i^th separated signal with only the j^th source active
%       fs - Sample frequency - used for plotting only
%      ddb - The distortion caused by BSS in dB
%
% The sizes of `x' and `y'
% Note: it is assumed that the j^th source signal is recovered at the
%       j^th output.
%
% See also: http://www.ele.tue.nl/ica99/
%      and: "Evaluation of Blind Signal Separation", in proc.
%            int. workshop on Independent Component Analysis and Blind
%            Signal Separation"

% Author: D.W.E. Schobben ([e-mail])
% rev. 1.03 @ 990108

if nargin ~=3,
  error('Please supply 3 arguments')
end

if length(size(x)) ~= 3 | length(size(y)) ~= 3,
  error('Please supply two 3D arrays')
end

J = size(x,2);

if size(x,3) ~= J | size(y,2) ~= J | size(y,3) ~= J | size(x,1) ~= size(y,1),
  error('Incorrect array size')
end

reduc=4;			% reduce the number of points in the plot with `red'

for j=1:J,
  for k=1:J,
    xplot=abs(fft(x(:,j,k)));
    yplot=abs(fft(y(:,j,k)));
    plotdata(:,1,j,k)=resample(20*log10(xplot(1:length(xplot)/2)),1,reduc);
    plotdata(:,2,j,k)=resample(20*log10(yplot(1:length(yplot)/2)),1,reduc);
  end
end

p_x_axis=linspace(0,fs/2,size(plotdata,1));
p_axis=[10 fs/2 min(min(min(min(plotdata(2:length(plotdata)))))) max(max(max(max(plotdata))))];
p_axis=[10 5000 85 125];

figure; shg; hold;

% Plot the data

% Plot contribution of source j to microphone
for j=1:J,
  subplot(3,J,j);
  semilogx(p_x_axis,plotdata(:,1,j,j));
  axis(p_axis);
%  xlabel('Frequency (Hz)');
  ylabel('Magnitude (dB)');
  title(['mic ' char('0'+j) ' due to source ' char('0'+j)]);
  grid
end

% Plot contribution of source j to BSS output j
for j=1:J,
  subplot(3,J,j+J);
  semilogx(p_x_axis,plotdata(:,2,j,j));
  axis(p_axis);
%  xlabel('Frequency (Hz)');
  ylabel('Magnitude (dB)');
  title(['BSS output ' char('0'+j) ' due to source ' char('0'+j)]);
  grid
end

% Plot unwanted contributions from other sources to BSS output j
for j=1:J,
  subplot(3,J,j+2*J);
  plotdatatemp=sum(plotdata,4)-plotdata(:,:,:,j);
  semilogx(p_x_axis,plotdatatemp(:,2,j));
  axis(p_axis);
  xlabel('Frequency (Hz)');
  ylabel('Magnitude (dB)');
  title(['BSS output ' char('0'+j) ' due to all other sources']);
  grid
end