function [varargout] = dcblock(varargin);
%
% DCBLOCK   determines the filter coefficient a for the dc blocking/high-
%           pass filter y(t) = x(t) - x(t-1) + a*y(t-1).
%
%           [a] = dcblock(Fc) returns the filter coefficient a for
%           normalized cut-on frequency Fc.
%
%           [a] = dcblock(fc,fs)returns the filter coefficient a for cut-on
%           frequency fc, and sampling frequency fs, when both are given in
%           Hz.
%
%           [aQ] = dcblock(fc,fs,B) returns the quantized filter
%           coefficient aQ when the quantization level B is given.
%
%           dcblock(...) plots the frequency and phase response of the
%           filter.
%
%           Note that if two parameters are passed to the function without
%           an output argument, the plot routine always interprets the
%           first as fc and the second as fs.
%
% SYNTAX    [a] = dcblock(Fc);
%           [a] = dcblock(fc,fs);
%           [aQ] = dcblock(fc,fs,B);
%           [Fc,fc] = dcblock(a,fs);
%           dcblock(Fc)
%           dcblock(fc,fs)
%           dcblock(fc,fs,B)
%
% INPUT     Fc      Normalized frequency cut-on of the filter.
%                   Fc = fc[Hz]/(fs[Hz]/2); fs/2 is the Nyquist frequency.
%                   Note that 0 < Fc < 1/3;
%
%           fc      filter cut-on frequency in Hz (or other reciprocal
%                   units of frequency e.g. 1/metre).
%
%           fs      sampling frequency in Hz (or other reciprocal
%                   units of frequency e.g. 1/metre).
%
%           B       quantization level for fixed point processing e.g. B is
%                   16 for 16-bit operation.
%
% OUTPUT    a       filter coefficient; 0 < a < 1 for floating point
%                   operation.
%
%           aQ      quantized filter coefficient a; 0 < aQ < 2^(B-1) - 1.
%                   The quantization level is set by B which is assumed to
%                   be signed integer.
%
% EXAMPLES  (1) Calculate the filter coefficient a when the normalized
%               frequency is 1E-5:
%
%           >> format long                      % set display format
%           >> a = dcblock(1E-5)                % Fc = 0.00001
%           >> a = 0.99994558749947             % filter coefficient
%
%           (2) Calculate the normalized and true cut-on frequency when the
%               coefficient a is 0.9987 and the sampling frequency is 1MHz:
%
%           >> [Fc,fc] = dcblock(0.9987,1E6)    % a = 0.9987; fs = 1MHz
%           >> Fc = 2.3906E-004                 % this output is normalized
%           >> fc = 119.5323                    % this output is in Hz
%
%           Note the difference in inputs between [Fc,fc] = dcblock(a,fs),
%           the function call [a] = dcblock(fc,fs), and its plot facility
%           dcblock(fc,fs). You can plot the filter response using
%
%           >> dcblock(119.5323,1E6)            % fc = 119.5323; fs = 1E6
%
%           (3) Carry out filtering operation using the dc blocking filter
%               with cut-on frequency of 70Hz and sampling rate of 10kHz.
%
%           t = -0.1:1/10000:0.1-1/10000;       % set up time base
%           s = sin(2*pi*2.5*t+pi/6)+0.01*randn(size(t))+1.25;
%           x = gauspuls(t,250,0.6) + s;        % noisy Gaussian test pulse
%           p = dcblock(70,10000);              % get filter coefficient
%           b = [1 -1];                         % set up differentiator
%           a = [1 -p];                         % set up integrator
%           y = filter(b,a,x);                  % filter test data
%           plot(t,x,t,y)                       % plot results to screen
%           legend('Original','DC Filtered','location','northwest');
%
%           The forward pass of the filter creates a distortion in the
%           Gaussian pulse. This distortion can be removed by using the
%           filtfilt form of the filter. If you now try
%
%           y2 = filtfilt(b,a,x);               % filter test data
%           plot(t,x,t,y,t,y2)                  % plot results to screen
%           legend('Original','DC Filtered','Undistorted','location','northwest');
%
%           you will notice that the symmetry of the Gaussian pulse has
%           been restored.
%

% Author:   J M De Freitas
%           QinetiQ Ltd, Winfrith Technology Centre
%           Winfrith, Dorchester DT2 8XJ, UK.
%
% Email     jdefreitas@qinetiq.com
%
% Revision  1.0
%
% Date      11 October 2006.
%

% Perform checks on input and output
Lin = nargin;
Lout = nargout;
if Lin > 3
error('dcblock error: Too many input arguments.');
end
if Lin < 1
errror('dcblock error: No input argument used.');
end
if Lout > 2
error('dcblock error: Too many output arguments.');
end

% Calculate a when Fc is given and plot response if required
if Lin == 1
Fc = varargin{1};
if (Fc < 0)|(Fc > 1/3)
error('dcblock error: Normalized cut-on frequency must be between 0 and 1/3.');
end
a = coef(Fc);
if Lout == 1
varargout{1} = a;
elseif Lout == 0
PlotResp(1,a);
end
end

% Calculate a when fc and fs are given and plot response if required
% or Fc and fc when a and fs are given
if Lin == 2
if (Lout == 1)|(Lout == 0)
fc = varargin{1};
fs = varargin{2};
Fc = 2*fc/fs;
if (Fc < 0)|(Fc > 1/3)
error('dcblock error: Cut-on frequency must be between 0 and fs/6.');
end
a = coef(Fc);
if Lout == 1
varargout{1} = a;
else
PlotResp(fs/2,a);
end
elseif Lout == 2
a = varargin{1};
fs = varargin{2};
if (a < 0)|(a > 1)
error('dcblock error: Filter coefficient a must be between 0 and 1.');
end
Fc = atan(sqrt(3)*(1-a^2)/(1+4*a+a^2))/pi;
varargout{1} = Fc;
varargout{2} = Fc*fs/2;
end
end

% Calculate a when B, fc and fs are given and plot if required
if Lin == 3
fc = varargin{1};
fs = varargin{2};
B = floor(varargin{3});
Fc = 2*fc/fs;
if (Fc < 0)|(Fc > 1/3)
error('dcblock error: Cut-on frequency fc must be between 0 and fs/6.');
return
end
if (B <= 2)|(B > 300)
error('dcblock error: Quantization level B must be between 2 and 300.');
return
end
a = coef(Fc);
if Lout == 1
varargout{1} = floor(a*(2^(B-1)-1));
elseif Lout == 0
aQ = floor(a*(2^(B-1)-1))/(2^(B-1)-1);
PlotResp(fs/2,aQ);
end
end

function a = coef(Fc);
% get up coefficient
a = (sqrt(3) - 2*sin(pi*Fc))/(sin(pi*Fc) + sqrt(3)*cos(pi*Fc));

function PlotResp(Fn,a);
% Plot frequency and phase response
f = 0.000001:0.000001:1-0.000001;
h = (1+a)*sin(pi*f/2)./sqrt(1+a^2-2*a*cos(pi*f));
phi = atan((1-a)*(sin(pi*f)./(1-cos(pi*f)))/(1+a));
indigo = [.48 .06 .89];
subplot(2,1,1);
loglog(Fn*f,h,'LineWidth',1.5,'Color',indigo);
title('\bf Frequency Response: IIR DC Blocking Filter','Color','k');
if Fn == 1
xlabel('Normalized Frequency');
else
xlabel('Frequency (Hz)');
end
ylabel('Magnitude ');
minY = min(get(gca,'YTick'));
minX = min(get(gca,'XTick'));
maxX = max(get(gca,'XTick'));
if minY >= 0.1
TickY = [0.1 1];
else
for j = 1:abs(log10(minY))+1 TickY(j) = 10^(log10(minY)+j-1); end
end
for j = 1:(log10(maxX)-log10(minX))+1 TickX(j) = 10^(log10(minX)+j-1); end
ylim([TickY(1) 1.1]);
set(gca,'XGrid','on','YGrid','on','XMinorGrid','on',...
'YMinorGrid','on','GridLineStyle','-',...
'MinorGridLineStyle','-','Xtick',TickX,'Ytick',TickY,...
'Color',[0.835 0.918 0.918]);
subplot(2,1,2);
loglog(Fn*f,phi*180/pi,'LineWidth',1.5,'Color',indigo);
title('\bf Phase Response: IIR DC Blocking Filter','Color','k');
if Fn == 1
xlabel('Normalized Frequency');
else
xlabel('Frequency (Hz)');
end
ylabel('Phase (degree) ');
ylim([0.1 100]);
minX = min(get(gca,'XTick'));
maxX = max(get(gca,'XTick'));
for j = 1:(log10(maxX)-log10(minX))+1
TickX2(j) = minX*10^(j-1) ;
end
set(gca,'XGrid','on','YGrid','on','XMinorGrid','on',...
'YMinorGrid','on','GridLineStyle','-',...
'MinorGridLineStyle','-','Xtick',TickX2,'YTick',[0.1 1 10 100],...
'Color',[0.835 0.918 0.918]);