[SP2, mean_array2]=sub_mean(SP, Fs, pps)

function [SP2, mean_array2]=sub_mean(SP, Fs, pps)
% % This program calculates the running average for sound and vibrations
% % data.  The output has the same size as the input.  The sampling rate 
% % and number of data averages per seconds greatly affect the output. 
% % 
% Example='';
% SP=rand(2,5000);  % Pa sound pressure time record 
%                   % 
% Fs=50000;         % Hz smpling rate for sound data
% pps=25;           % stands for points per seconds
%                   % number of data averages in data points per second 
%                   % milliseconds per data point is 1000/pps
%                   % default is pps=25; 40 ms per data point 
% 
% [SP2, mean_array2]=sub_mean(SP, Fs, pps);
%
% % 
% % Output Variables
% % 
% % SP2 is the sound pressure with the running average subtracted  
% % SP2 has the same size as the input array
% % mean_array2 is the running mean.  
% % mean_array2 has the same size as the input array
% % 
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Program Written by Edward L. Zechmann
% %    
% %             Date    15 August 2007
% %         modified    18 December 2007    
% %                     added comments
% %                     added transpose SP to have
% %                     size(num_channel,num_data)
% %                     converted pps into an input variable       
% %                     
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Please Feel Free to Modify This Program
% %    

% Make sure SP has the correct order of the 
% channels and sound pressure data
[m1, n1]=size(SP);
flag1=0;

if m1 > n1
    flag1=1;
    SP=SP';
    [m1 n1]=size(SP);
end

SP2=zeros(m1, n1);

if nargin < 3
    pps=25;
end

% limit the number of bins to speed up processing the data
if nargin < 2
    bins=1000;
else
    bins=ceil(pps*n1/Fs); %25 Hz sampling rate, 40 ms per bin
    if bins > 100000
        bins=100000; % limit the number of bins to 10000.
    end
end

r1=floor(n1/bins);
mean_array2=zeros(m1, n1);

if r1 > 10
    
    for e1=1:m1;
        
        mean1=zeros(1, bins);
        array_bins=zeros(1, bins);
        for e2=1:(bins-1);
            
            SPbuffer=SP(e1, (e2-1)*r1+([1:r1]));
            mean1(e2)=mean(SPbuffer');
            array_bins(e2)=(e2-1)*r1+1;
            
        end
        
        e2=bins;
        SPbuffer=SP(e1, ((e2-1)*r1+1):end);
        mean1(e2)=mean(SPbuffer');
        array_bins(e2)=(e2-1)*r1+1;
        
        mean1=[mean1(1)*[1 1 1 1] mean1 mean1(end)*[1 1 1 1] ];
        array_bins=[[-3 -2 -1 0] array_bins (array_bins(end)+[1 2 3 4]) ];
        mean_array = interp1(array_bins,mean1,1:n1,'pchip');  
        SP2(e1,:)=SP(e1,:)-mean_array;
        mean_array2(e1,:)=mean_array;
        
    end
    
else
    SP2=SP-mean(SP')'*ones(1, n1);
end

% Make sure that the ouput has the same size as it was input
% channels and data have the smae indices 
if isequal(flag1, 1)
    SP2=SP2';
end