[yC, B, A]=Cweight_time_filter(y, Fs, B, A)

function [yC, B, A]=Cweight_time_filter(y, Fs, B, A)
% % This program applies an C-weighting filter to 
% % a sound pressure time record.  Then it outputs the C-weighted tme
% % record and the filter coefficients.
% % 
% % 
% Example='';
% y=rand(1, 50000); % (Pa) waveform
%                   % y should have the size (num_channels, num_datasample)
%                   % 
% Fs=50000;         % (Hz) Sampling rate 
% B=[];             % C-weighting filter coefficients
% A=[];             % C-weighting filter coefficients
% % 
% % Output Variables
% % yC is the C-weighted time record Pa
% % B is an array of C-weighting filter coefficients
% % A is an array of C-weighting filter coefficients
% % 
% % 
% % This program uses the signal processing toolbox 
% % 
% % This program uses functions from 
% % -- Christophe Couvreur ([e-mail]), August 1997. 
% % this is the octave set of filter functions from Matlab Central
% %
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Written by Edward L. Zechmann           
% %     date  June 2007
% %     modified 15 November 2007   Updated comments
% %     modified 19 December 2007   added a resampling routine 
% %                                 to improve accruacy at low and
% %                                 high freuencies 
% %                                 Updated comments
% %                                 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % 
% % Please feel free to modify this code.
% % 

% surpress an annoying warning due to a nearly singular matrix which almost always occurs with filtering.
warning('off','MATLAB:nearlySingularMatrix')

p=1;
q=1;
Fsn=Fs;
flag0=0;

if Fs > 50000;
    flag0=1;
    q=ceil(Fs/40000);
    Fsn=Fs/q;
end

[m1 n1]=size(y);
flag1=0;

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


if isequal(flag0, 1)
    buf=resample(y(1, :), p, q);
    n2=length(buf);
    y2=zeros(m1,n2);
    y2(1, :)=buf;
    clear('buf');
    if m1 > 1
        for e1=2:m1;
            y2(e1, :)=resample(y(e1, :), p, q);
        end
    end
else
    y2=y;
end

clear('y');

if (nargin < 4) || length(B) < 4 ||  length(A) < 4

    [B,A] = cdsgn(Fsn); 
end

yC=zeros(size(y2));
   
for e1=1:m1;
    yC(e1, :) = real(filter(B, A, y2(e1, :)));
end

if isequal(flag0, 1)
    buf=resample(yC(1, :), q, p);
    n2=length(buf);
    yCC=zeros(m1,n2);
    yCC(1, :)=buf;
    clear('buf');
    if m1 > 1
        for e1=2:m1;
            yCC(e1, :)=resample(yC(e1, :), q, p);
        end
    end
else
    yCC=yC;
end

clear('yC');

if isequal(flag0, 1)
    if n2 > n1
        n2=n1;
    end

    yC=yCC(1:m1, 1:n2);
else
    yC=yCC(1:m1, 1:n1);
end

%output yC with the same size as y was input
if isequal(flag1, 1)
    yC=yC';
end