DFiltInt (M, NCof, PSD, Delay, SFreq, g)

function [h, mse] = DFiltInt (M, NCof, PSD, Delay, SFreq, g)
% [h, mse] = DFiltInt (M, NCof, PSD, Delay, SFreq, g)
% h:  Filter coefficients (column vector)
% mse: Normalized mean-square error for each subfilter (M values)
%
% M:  Interpolation ratio
% NCof: Number of interpolating filter coefficients.
% PSD: Structure containing the power spectral values. Piecewise monotonic
%     cubic interpolation is used between the given values of the power
%     spectral density.
%     PSD.f: Frequency vector for the power spectral density (positive,
%       increasing values)
%     PSD.psd: Power spectral density values
%     PSD.fcos: Cosine frequencies
%     PSD.Pcos: Cosine powers
% Delay: Filter delay in samples at the sampling frequency after
%     interpolation. This value is not restricted to be integer. The
%     default value is (NCof-1)/2 - (Ng-1)/2, where Ng is the number of
%     coefficients in g. Changing the delay allows for non-symmetric
%     interpolating filters or filters which approximate a delay. Leave
%     this parameter empty to get the default value.
% SFreq: Sampling frequency for the sequence to be interpolated. The
%     sampling frequency for the interpolated sequence is this value
%     multiplied by the interpolating ratio. If the sampling frequency is
%     not specified, a normalized sampling frequency of one is used.
% g:  Reference filter, default single tap, value 1. The default is
%     appropriate for normal interpolating filters.
%
% This routine designs an interpolating filter which minimizes the mean-
% square interpolation error. The input to the interpolation process is a
% basic rate sampled signal. The output of the interpolation process is an
% increased rate sampled signal. The design of the interpolating filter is
% specified in terms of the power spectrum of the underlying continuous-
% time signal. This power spectrum is used to derive the auto-correlation
% functions for the sampled signals. The auto-correlation determines the
% mean-square error between the increased-rate signal and the samples of
% the underlying continuous-time signal.
%
% Conceptually, interpolation has two steps. In the first, the sampling
% rate  of the input sequence is increased by a factor of M by inserting
% M-1 zeros between each input sample. In the second step, the increased
% rate sequence is applied to the interpolating filter to form the
% interpolated sequence.
%
% The block diagram of the system is shown below. The lower branch is the 
% reference branch. It consists of a delay, a sampler and a reference
% filter g[n]. For ordinary interpolation, the reference filter is an
% identity filter. The upper branch interpolates a sampled signal. The
% interpolating filter h[n] is chosen to minimize the mean-square error.
%
%               -----      ------   ------       --------
%   x(t) --.--->| \ |----->| \M |-->| xM |------>| h[n] |---.--> y[n]
%          |    ----- x[n] ------   ------ xs[n] --------   |
%          |   sample    downsample upsample  interpolating |
%          |                                     filter     |
%          |  ---------   -----         --------          + v
%          -->| Delay |-->| \ |-------->| g[n] |----------> + --> e[n]
%             ---------   -----  xd[n]  --------   r[n]   - 
%                         sample        reference
%                                        filter
%
% This program generalizes the method outlined by Oetken, Parks, and
% Schussler.
% Reference:
%     G. Oetken, T. W. Parks and H. W. Schussler, "New Results in the
%     Design of Digital Interpolators", IEEE Trans. Acoustics, Speech,
%     Signal Processing, vol. 23, pp. 301-309, June 1975.
%
% The one-sided power spectrum of the input signal is specified in terms
% of a continuous power spectral density (PSD) and discrete (sinusoidal)
% components. The PSD is given as tabulated values at positive frequency
% values. For the continuous pat, piecewise monotonic cubic interpolation
% of the PSD is used between the given values. The sinusoidal components
% are specified by the frequencies and powers of the components.
%
% Values of the power spectrum beyond the last point specified are
% assumed to be zero. Note that if the last power spectrum point is
% explicitly set to zero, the slope at that point is also set to zero.
% To specify a power spectrum which is discontinuous, do not specify the
% frequency at which the spectrum goes to zero.
%
% Example: 23 tap linear phase interpolating filter:
%   The input spectrum is modelled by a power spectrum which is flat from
%   dc to 2400 Hz. The spectrum then decreases to become zero at 3200 Hz.
%   There is an additional dc component. The interpolating ratio is 4
%   (input 8000 Hz, output 32000 Hz).
%      PSD.f = [0 2400 3200];
%      PSD.psd = [1 1 0];    % Zero slope at the last point
%      PSD.fcos = 0;
%      PSD.Pcos = 1000;
%      h = DFiltInt (4, 23,PSD, [], 8000);

% $Id: DFiltInt.m,v 1.9 2009/07/06 17:45:24 pkabal Exp $

% Fill in missing parameters
if (~ exist ('g', 'var'))
  g = 1;
end
if (~ exist ('SFreq', 'var'))
  SFreq = 1;
end
if (~ exist ('Delay', 'var'))
  Delay = [];
end

% Reconcile the arguments
[PSD, Delay, g] = INPSD (M, NCof, PSD, Delay, SFreq, g);

% Design the filter
h = INdesInt (M, NCof, PSD, Delay, g);

% Calculate the mean-square error
mse = INintMSE (h, M, PSD, Delay, g);

% Summary printout
INprintSpec (M, NCof, PSD, Delay, g, mse);

return