DSPRelated.com
   Code Snippets

Interleave Zeros

Jeff T ● July 10, 2011●1 comment ● Coded in Matlab 
function out = interleave_zeros(a)

% Interleave vector with zeros.
% Only 1D matrices are supported.
%
% Usage: b = INTERLEAVE_ZEROS(a);
%
%        A is the vector that will have zeros inserted
%        B will contain the vector set A and the interleaved
%          zeros (twice the length)
%
% Author: sparafucile17 2/13/04

my_zeros = zeros(1, length(a));
out = matintrlv([a  my_zeros],2,length(a));

Zero Crossing Counter

Jeff T ● July 9, 2011●9 comments ● Coded in Matlab 
function count = zero_crossings(x)
% Count the number of zero-crossings from the supplied time-domain
% input vector.  A simple method is applied here that can be easily
% ported to a real-time system that would minimize the number of
% if-else conditionals.
%
% Usage:     COUNT = zero_crossings(X);
%
%            X     is the input time-domain signal (one dimensional)
%            COUNT is the amount of zero-crossings in the input signal
%
% Author:    sparafucile17 06/27/04
%

% initial value
count = 0;

% error checks
if(length(x) == 1)
    error('ERROR: input signal must have more than one element');
end

if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
    error('ERROR: Input must be one-dimensional');
end
    
% force signal to be a vector oriented in the same direction
x = x(:);

num_samples = length(x);
for i=2:num_samples

    % Any time you multiply to adjacent values that have a sign difference
    % the result will always be negative.  When the signs are identical,
    % the product will always be positive.
    if((x(i) * x(i-1)) < 0)
        count = count + 1;
    end
    
end

ML Estimation of Several Frequencies Using the Nonuniform FFT

Jesus Selva ● July 8, 2011 ● Coded in Matlab 
function st = Demo

% Author: J. Selva.
% Date: July 2011.
% Version: 0.

st = [];

n = input('Options: 1 perform single estimation, 2 estimate RMS error: ');

switch n
  case 1
    st = TestSingleTrial;
  case 2
    TestRMSError;
end

function st = TestSingleTrial

%Scenario with 5 frequencies. 

M = 2048;
SNRInSpecdB = 27; %Signal-to-noise ratio after correlation. 

Nf = 5;
df = 1/M; %Minimum separation among frequencies. 

MaxAmp = 1; %Maximum amplitude.
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;

[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);

y = DataRealization(vf,va,M,nVar);

p = Preprocessing(y,nVar,PFA);

disp('True frequencies: ')
disp(vf)

disp('True amplitudes: ')
disp(va)

disp(['Post-correlation Signal-to-noise ratio: ',...
      num2str(SNRInSpecdB),' dB'])

s = input('To compute ML estimate, use (S)erial  or (M)anual method? ','s');

if strcmp(lower(s),'s')
  st = MethodSerial1(p);
else
  st = MethodManual1(p);
end

disp('ML frequency estimates: ')
disp(st.Freqs)

disp('ML amplitude estimates: ')
disp(st.Amps)

function TestRMSError

%Estimates the RMS error and computes the Cramer-Rao bounds.

format compact
format short g

M = 2048;
SNRInSpecdB = 27;

Nf = 5;
df = 1/M;

MaxAmp = 1;
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;

[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);

crb = CRBBound(vf,va,M,nVar);

Ntr = 0;
S2 = zeros(Nf,1);
for k = 1:1e5

  y = DataRealization(vf,va,M,nVar);

  p = Preprocessing(y,nVar,PFA);

  st = MethodSerial1(p);
  
  if length(st.Freqs) == 5
    Ntr = Ntr + 1;
    S2 = S2 + (vf-st.Freqs).^2;
  end
  
  if mod(k,100) == 0
    disp(['Number or trials: ', num2str(k)])
    disp(['             Frequencies:      ' num2str(vf.')])
    disp(['Deviation CRB bounds (dB):   ' num2str(20*log10(crb.'))])
    disp(['          RMS errors (dB):   ' num2str(10*log10(S2(:).'/Ntr))])
    disp(' ');
  end
end

function [vf,va] = RandomFreqAndAmpVectors(Nf,MinFreqSep,MinAmp,MaxAmp)

%Constructs random vectors of frequencies and amplitudes. The module of the largest
%amplitude is MaxAmp. There are no amplitudes with module below MinAmp, and the minimum
%separation among frequencies is MinFreqSep.

vf = sort(rand(Nf,1));
vf = vf + MinFreqSep*(cumsum(ones(Nf,1))-1);
vf = vf/(1+(Nf-1)*MinFreqSep) - 0.5;

pr = rand(Nf,1);
ind = floor(rand(1)/Nf)+1;
pr(ind) = 1;

va = (MinAmp+(MaxAmp-MinAmp)*pr) .* exp(1i*2*pi*rand(Nf,1));

function y = DataRealization(vf,a,M,nVar)

m1 = FirstIndex(M);

y = exp(1i*2*pi*(m1:m1+M-1).'*vf(:).')*a(:)+randn(M,2)*[1 ; 1i]*sqrt(nVar/2);

function p = Preprocessing(y,nVar,PFA)

%Performs the preprocessing. Fundamentally this consists in computing the FFT inside 
%"Correlation.m".
%The input parameters are
%
% y: Data vector.
% nVar: Noise variance estimate. 
% PFA: Probability of seeing a frequency when there isn't any. Usually a very low value
%      like 1e-5.

p = [];
p.M = length(y);
p.nVar = nVar;
p.PFA = PFA;
p.nThres = NoiseThres(p.M,p.nVar,p.PFA); %Noise threshold.

icf = Correlation(y,p.nThres); %Struct with data for interpolating the correlation.

p.cfd = Cost(real(y(:)'*y(:)),icf); %Struct with data for computing the ML cost function ...
                                    %and its derivatives. 

p.freqThres = 1/(10*p.M); %If two frequencies differ in less than ...
                          %freqThres  in  the  iterative  search, then   one  of  them  is
                          %discarded. This prevents ill-conditioning.

p.freqStopThres = 1/(p.M*1e4); %Threshold for stopping the ML cost function minimization.

function st = MethodSerial1(p)

%Adds frequencies and computes ML estimates sequentially.

vfEst = [];
ModifiedEstimation = logical(1);

while ModifiedEstimation
  
  [vf1,va1] = Correlation(p.cfd.icf,'gridMainLocalMaxima',vfEst);

  ModifiedEstimation = ~isempty(vf1);

  if ModifiedEstimation
    vfEst(end+1) = vf1(1);
    [listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,...
	10.^-8,p.freqThres,p.freqStopThres);
    vfEst = sort(listv(:,end));
  end 
end

st = Report(p,vfEst);

function st = MethodManual1(p,varargin)

%Allows the user to set the initial frequency estimates on a plot.

if nargin == 1
  vfEst = [];
else
  vfEst = varargin{1};
end

vfEst = vfEst(:);

figure
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
yl = get(gca,'YLim');
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

while 1
  pr1 = input('A/a -> Add freq., D/d -> Delete freq, X/x -> exit: ','s');
  
  switch lower(pr1)
    case 'a'
      [f1,pr] = ginput(1);
      vfEst = sort([vfEst; f1]);
      [listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,10.^-8,... 
	  p.freqThres,p.freqStopThres);
      vfEst = sort(listv(:,end));

      subplot(2,1,1)
      Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);

      subplot(2,1,2)
      Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

    case 'd'
      [f1,pr] = ginput(1);
      [pr,ind] = min(abs(f1-vfEst));
      vfEst(ind) = [];

      subplot(2,1,1)
      Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);

      subplot(2,1,2)
      Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

    case 'x'
      
      st = Report(p,vfEst);
      return
  end
end

function thres = NoiseThres(M,nVar,PFA)

%Computes the noise threshold. The method will discard any correlation local maximum below
%this threshold.

thres = sqrt((M*nVar/2)*icdf('chi2',(1-PFA)^(1/M),2));

function m1 = FirstIndex(M)

m1 = -ceil(M/2);

function varargout = Correlation(varargin)

%Computes the correlation value and its derivatives by calling "Barycentric". It does
%other things like plotting. 

if ~isstruct(varargin{1})
  st = [];
  temp = {[],[],10.^-15,1.5};
  [temp{1:nargin}] = deal(varargin{:});
  [a,st.nThres,epsilon,OvFMin] = deal(temp{:});
  a = a(:).';
  
  st.M = length(a);
  st.K = 2^nextpow2(st.M*OvFMin);
  st.bary = Barycentric(0,SampleFourierCorrelation1DVer1(a,st.K),...
           -2*FirstIndex(st.M),1/st.K,whatP(1/st.K,-2*FirstIndex(st.M),epsilon));
       

  st.IndActive = find(abs(st.bary.a) > st.nThres);
  st.IndActive = unique([st.IndActive-1,st.IndActive,st.IndActive+1]);
  if st.IndActive(1) == 0
    st.IndActive(1) = [];
  end

  if st.IndActive(end) == st.K+1
    st.IndActive(end) = [];
  end
  
  st.aTActive = st.bary.a(st.IndActive);
  st.vfActive = centreFreq((st.IndActive-1)/st.K);
  varargout = {st};
  return
end 

[st,msg] = deal(varargin{1:2});
  
switch msg
    
  case 'value'
    varargout = {Barycentric(st.bary,varargin{3:end})};
    
  case 'gridValueInRange'
      
    [ind1,ind2] = deal(varargin{3:4});
    varargout = {st.bary.a(mod(ind1:ind2,st.K)+1)};
     
  case 'ampEstimates'
      
    vf = varargin{3};
    if length(vf) == 0
      varargout = {[]};
    end

    varargout = {RegGeoSum(st.M,vf,vf,0,0) \ Barycentric(st.bary,vf)};
      
  case 'gridValue'
    
    lInd = varargin{3};
    if nargin == 4
      vf = varargin{4};
    else
      vf = [];
    end
      
    if isempty(vf)
      varargout = {st.bary.a(mod(lInd,st.K)+1)};
    else
      vf1 = centreFreq(lInd/st.K);
      
      VRaa = RegGeoSum(st.M,vf,vf,0,0);
      VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
      
      pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
      pr = reshape(pr,size(lInd));

      varargout = {st.bary.a(mod(lInd,st.K)+1)-pr};
    end

  case 'gridFit'
    
    lInd = varargin{3};
    if nargin == 4
      vf = varargin{4};
    else
      vf = [];
    end
      
    if isempty(vf)
      varargout = {zeros(size(lInd))};
    else
      vf1 = centreFreq(lInd/st.K);
      
      VRaa = RegGeoSum(st.M,vf,vf,0,0);
      VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
      
      pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
      pr = reshape(pr,size(lInd));

      varargout = {pr};
    end
    
  case 'plotErrorIndB'

    if nargin == 2
      vf = [];
    else 
      vf = varargin{3};
    end
      
    s = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vf)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s = s(ord);
      
    plot(gf,s,[-0.5,0.5],20*log10(abs(st.nThres)*[1,1]))
      
    Ms = max(s);
    set(gca,'YLim',[Ms-50,Ms+10])
    xlabel('Frequency (f)')
    ylabel('Error spectrum (dB)')
    grid on
    
    varargout = {};
      
  case 'PlotPeriodogramFit'
    
    vfEst = varargin{3};
    
    if isempty(vfEst)
      
       s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
       gf = centreFreq((0:st.K-1)/st.K);
       [gf,ord] = sort(gf);
       s0 = s0(ord);
      
       plot(gf,s0)
       hold on
       nt = 20*log10(abs(st.nThres));
       plot([-0.5,0.5],nt([1,1]).',...
	   'LineWidth',2,'LineStyle','--','Color',[1 0 0])
       hold off
       text(-0.49,nt+2,'Noise threshold')
       
       xlabel('Frequency (f)')
       ylabel('Amplitude (dB)')
       title('Periodogram data, fit, and est. frequencies (stems)')
       return
    end
    
    vaEst = Correlation(st,'ampEstimates',vfEst);
    
    stem(vfEst,20*log10(st.M*abs(vaEst)),'g')
    
    hold on 
    
    s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s0 = s0(ord);
      
    plot(gf,s0,'b')
    hold on
    nt = 20*log10(abs(st.nThres));
    plot([-0.5,0.5],nt([1,1]).',...
	'LineWidth',2,'LineStyle','--','Color',[1 0 0])
      
    s1 = 20*log10(abs(Correlation(st,'gridFit',0:st.K-1,vfEst)));
    s1 = s1(ord);
      
    plot(gf,s1,'r')

    Ms = max(s0);
    set(gca,'YLim',[Ms-50,Ms+10])
    
    xlabel('Frequency (f)')
    ylabel('Amplitude (dB)')
    title('Periodogram fit')
    
    hold off
    
    varargout = {};
    
  case 'PlotPeriodogramFitError'
    
    if nargin == 2
      vfEst = [];
    else
      vfEst = varargin{3};
    end
    
    vaEst = Correlation(st,'ampEstimates',vfEst);
    
    s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vfEst)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s0 = s0(ord);
     
    nt = 20*log10(abs(st.nThres));
    
    plot(gf,s0)
    hold on 
    plot([-0.5,0.5],nt([1,1]).',...
	'LineWidth',2,'LineStyle','--','Color',[1 0 0])
    hold off
    
    Ms = max(s0);
    set(gca,'YLim',[Ms-50,Ms+10])
    
    xlabel('Frequency (f)')
    ylabel('Amplitude (dB)')
    title('Periodogram error')
    
    varargout = {};
    
  case 'gridValueGivenAmps'
      
    [lInd,vf,aEst] = deal(varargin{3:5});
    vf1 = centreFreq(lInd/st.K);
    VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
    
    pr1 = st.bary.a(mod(lInd,st.K)+1);
    pr1 = reshape(pr1,size(lInd));
    pr2 = VRa1a*aEst;
    pr2 = reshape(pr2,size(lInd));
      
    varargout = {pr1-pr2};
      
  case 'gridLocalMaxima'
      
    if nargin == 2

      pr = [-2,-1,0].';
      
      pr = abs(st.bary.a(mod(st.IndActive([1;1;1],:) + ...
             pr(:,ones(1,length(st.IndActive))),st.K)+1));

      pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
      
      vf = st.vfActive(pr);
      va = st.aTActive(pr);
   
      varargout = {vf,va};
      return
    end
      
    vf = varargin{3};
      
    if length(vf) == 0
      [varargout{1:2}] = Correlation(st,'gridLocalMaxima');
      return
    end
   
    amp = Correlation(st,'ampEstimates',vf);

    tmp = abs(Correlation(st,...
	    'gridValueGivenAmps',st.IndActive-1,vf,amp)) >= st.nThres;
     
    ind = st.IndActive(tmp);
    
    pr = [-2,-1,0].';
    
    ind1 = ind([1,1,1].',:) + pr(:,ones(1,length(ind)));
    
    pr = abs(Correlation(st,'gridValueGivenAmps',ind1,vf,amp));

    pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
    
    ind = ind(pr)-1;
   
    if isempty(ind)
      varargout = {[],[]};
    else
      varargout = {centreFreq(ind/st.K),Correlation(st,...
	  'gridValueGivenAmps',ind,vf,amp)};
    end

  case 'gridMainLocalMaxima'
      
    [vf,va] = Correlation(st,'gridLocalMaxima',varargin{3:end});
    if length(vf) == 0
      varargout = {[],[]};
      return
    end
	  
    [varargout{1:2}] = selectLocalMaximaVer3(vf,va,st.nThres,st.M);
	
      
  case 'gridMaximum'
      
    pr = Correlation(st,'gridValue',0:st.K-1,varargin{3:end});
    [ma,ind] = max(abs(pr));
      
    varargout = {(ind-1)/st.K,ma};
      
end

function tab = SampleFourierCorrelation1DVer1(a,K)  

M = length(a);
m1 = FirstIndex(M);

a = a(:).';
tab = fft([a(-m1+1:end),zeros(1,K-M),a(1:-m1)]);

function out = centreFreq(f)

out = f + ceil(-0.5-f);

function P = whatP(T,B,epsilon)

P = ceil(acsch(epsilon)/(pi*(1-B*T)));

function d = dist(x,y)

if x<y
  d = min([y-x,x+1-y]);
else
  d = min([x-y,y+1-x]);
end

function v = sincMask(x,y,M)

v = zeros(size(x));
for k = 1:length(x)
  x1 = max([dist(x(k),y)-1/(2*M),0]);
  v(k) = 1/max([M*sin(pi*x1),1]); 
end

function v = slBoundVer2(x,vam,vfe,nThres,M1)

vam = vam(:);
vfe = vfe(:);

sll = sincMask(vfe,x,M1).'*vam;

v = max([sll+nThres,1.1*sll]);

function varargout = selectLocalMaximaVer3(vf,va,nThres,M)

va1 = abs(va(:));
[pr,ord] = sort(va1);
ord = flipud(ord);

vf1 = vf(ord);
va1 = va1(ord);

n1 = length(va1);

if va1(1) < nThres 
  varargout = {[],[]};
  return
end

va2 = va1(1);
va2Abs = va1(1);
vf2 = vf1(1);

for n=2:n1
  if va1(n)<nThres
    continue
  end
  if va1(n) < slBoundVer2(vf1(n),va2Abs,vf2,nThres,M)
    continue
  end
  vf2(end+1) = vf1(n);
  va2(end+1) = va1(n);
  va2Abs(end+1) = abs(va(ord(n)));
end

varargout = {vf2,va2};

function out = Barycentric(varargin)

%Performs  barycentric interpolation using the method in [Bary].

if ~isstruct(varargin{1})
  st = [];
  [st.n0,st.a,st.B,st.T,st.P] = deal(varargin{:});
  
  st.N = length(st.a);
  st.alpha = baryWeights(st.T,st.B,st.P);
  st.alpha = st.alpha([1:st.P,st.P+2:end,st.P+1]);
  
  st.vt = (-st.P:st.P)*st.T; 
  st.vt = st.vt([1:st.P,st.P+2:end,st.P+1]);

  out = st;
  
else
  
  [st,v] = deal(varargin{1:2});
  
  if nargin == 3
    nd = varargin{3};
  else
    nd = 0;
  end
 
  v = v(:);
  Nv = numel(v);

  n = floor(v/st.T+0.5);
  u = v - n*st.T;
  n1 = n - st.n0;
    
  pr = [-st.P:-1,1:st.P,0];
  da = st.a(mod(n1(:,ones(1,2*st.P+1)) + pr(ones(Nv,1),:),st.N)+1);
  
  out = DerBarycentricInterp3Vec(st.alpha,da,st.vt,u,nd);
     
end

function w = baryWeights(T,B,P)

vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;

N = length(vt);
LD = ones(1,N);

for k = 1:N
  LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end

w = gam ./ LD;
w = w / max(abs(w));

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)

NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
             %usage.

NBlock = 1+floor(NtauL/LBtau);

t = t(:).';
Nt = length(t);

out = zeros(NtauL,nd+1);

for r = 0:NBlock-1
  ind1 = 1+r*LBtau;
  ind2 = min([(r+1)*LBtau,NtauL]);
 
  Ntau = ind2-ind1+1;
  z = zL(ind1:ind2,:);
  tau = tauL(ind1:ind2);
 
  vD = zeros(Ntau,1);
 
  LF = [1,zeros(1,Nt-1)];
  LF = LF(ones(Ntau,1),:);

  for k = 1:Nt-1
    vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
    LF(:,k+1) = LF(:,k) .* (tau-t(k));
  end

  vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);

  for kd = 0:nd

    pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);

    for k = 1:Nt-1
      cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
    end
    cN = cN ./ vD;
 
    ztau = z(:,end)+(tau-t(end)) .* cN;
    out(ind1:ind2,kd+1) = ztau;
 
    if kd < nd
      pr1 = ztau;
      pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
     
      pr2 = t(1:end-1);
      pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
     
      z = [pr1 ./ pr2, cN] * (kd+1);
     
    end
  end
 
end

function varargout = Cost(varargin)

%Computes the ML cost function differentials by calling "Correlation" and "RegGeoSum". 

if ~isstruct(varargin{1})
  st = [];
  [st.EnergyRyy,st.icf] = deal(varargin{:});
  varargout = {st};
  return
end

st = varargin{1};

if nargin > 2
  vf = varargin{3};
else
  vf = [];
end

switch varargin{2}
  
  case 'value'
    
    vf = varargin{3};
    Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
    Rsd = RegGeoSum(st.icf.M,vf,vf,0,1);
    Rdd = RegGeoSum(st.icf.M,vf,vf,1,1);
    
    pr = Correlation(st.icf,'value',vf,1);
    Rsy = pr(:,1);
    Rdy = pr(:,2);
    
    pr = Rss\[Rsy,Rsd];
    MRsy = pr(:,1);
    MRsd = pr(:,2:end);
    
    varargout = {st.EnergyRyy-RealTrace(Rsy',MRsy),-2*RealDiag(MRsy,Rdy'-MRsy'*Rsd),...
	2*RealHadamard((Rdd'-Rsd'*MRsd).',MRsy*MRsy')};

  case 'value0'
  
    if length(vf) == 0
      varargout = {st.EnergyRyy};
      return
    end

    vf = varargin{3};
    Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
    Rsy = Correlation(st.icf,'value',vf);
    
    [varargout{1}] = st.EnergyRyy-RealTrace(Rsy',Rss\Rsy);
    
end

function v = RealHadamard(A,B)

v = real(A).*real(B)-imag(A).*imag(B);

function C = RealDiag(A,B)

% Same as real(diag(A*B)) but using less flops.

C = RealHadamard(A.',B);

if size(C,1) == 1
  C = C.';
else
  C = sum(C).';
end

function C = RealTrace(A,B)

% Same as real(trace(A*B)) but using less flops.

C=RealHadamard(A.',B);
C=sum(C(:));

function [listv,listLv,listCL,ErrRatio] = ...
    LookForMinimum(L,v0,nThres,freqThres,freqStopThres,varargin)

%Executes a descent algorithm using a method similar to that in [ML].

if ~isempty(varargin) 
  nItMax = varargin{1};
else
  nItMax = 20;
end

v0 = v0(:);
Nf = numel(v0);

listv = [v0,zeros(Nf,nItMax)];
listLv = [Cost(L,'value0',v0),zeros(1,nItMax)];
listCL = [NaN,zeros(1,nItMax)];

succ = 1;
nIt = 1;

while succ && nIt <= nItMax
  [succ,v,Lv,CL] = NewtonIt(L,listv(:,nIt),freqThres);
  
  if succ 
    nIt = nIt + 1;
    listv(:,nIt) = v;
    listLv(nIt) = Lv;
    listCL(nIt) = CL;
  else
    break
  end
    
  if max(abs(listv(:,nIt)-listv(:,nIt-1))) < freqStopThres
    break 
  end
end

ErrRatio = (Cost(L,'value0')-listLv(end))/listLv(end);
listv(:,nIt+1:end) = [];
listLv(nIt+1:end) = [];
listCL(nIt+1:end) = [];

function [v,vf] = validVector(vf,thres)

vf = sort(vf(:));
if length(vf)==1 
  v = 1;
  return
end

m = vf(2:end)-vf(1:end-1) < thres;

if any(m)
  v = 0;
  vf([logical(0); m]) = [];
else
  v = 1;
end

function [succ,v1,Lv1,nCutsLeft] = NewtonIt(L,v0,thres)
  
Lv1 = [];
nCutsLeft = [];

[vv,v1] = validVector(v0,thres);

if ~vv
  succ = 0;
  return
end

[Lv0,g0,h0] = Cost(L,'value',v0);

d = -h0\g0;

nCutsLeft = 20;

while nCutsLeft > 0
  [vv,v1] = validVector(v0+d,thres);
  if ~vv
    succ = 0;
    v1 = v0;
  end
  
  Lv1 = Cost(L,'value0',v0+d);
  
  if Lv1 < Lv0
    succ = 1;
    v1 = v0 + d;
    break;
  end
  
  nCutsLeft = nCutsLeft - 1;
  d = d/2;
end

if nCutsLeft == 0
  succ = 0;
end

function out = RegGeoSum(varargin)

%Computes the value of a geometric sum or its differentials using the exact formula. It
%takes care of the preventable singularity at zero.

switch nargin
  case 2
    [M,vf] = deal(varargin{:});
    nd = 0;
    msg = 'value';
    
  case 3
    [M,vf,nd] = deal(varargin{:});
    msg = 'value';
  case 5
    [M,vf1,vf2,nd1,nd2] = deal(varargin{:});
    msg = 'corr';
end

switch msg
  
  case 'value'
    out = RegGeoSum1(M,vf,nd);
    
  case 'corr'
    vf1 = vf1(:);
    Nf1 = length(vf1);
    vf2 = vf2(:);
    Nf2 = length(vf2);
    
    out = RegGeoSum1(M,-vf1(:,ones(1,Nf2))+vf2(:,ones(1,Nf1)).',nd1+nd2)*(-1)^nd1;
end

function v = RegGeoSum1(M,vf,nd)

m1 = -ceil(M/2);

ind = abs(vf)>eps^(1/3);

vf1 = vf(ind);

v = zeros(size(vf));

switch nd
  
  case 0

    v(ind) = -(exp(1i*2*pi*m1*vf1)-exp(1i*2*pi*(m1+M)*vf1))./(-1+exp(1i*2*pi*vf1));
    
    v(~ind) = M;
  
  case 1
    
    v(ind) = (-2*1i*pi*(  exp(1i*2*pi*vf1*(1 + m1))*(-1 + m1) - exp(1i*2*pi*vf1*m1)*m1  ...
	-exp(1i*pi*2*vf1*(1 + M + m1))*(-1 + M + m1) + ...
       exp(1i*2*pi*vf1*(M + m1))*(M + m1)))./(-1 + exp(1i*2*pi*vf1)).^2;
   
    v(~ind) = pi*1i*M*(-1+M+2*m1);
    
  case 2
    
    v(ind) = (4*(exp(1i*2*pi*vf1*(2 + m1))*(-1 + m1)^2 + exp(1i*2*pi*vf1*m1)*m1^2 - ...
       exp(1i*2*pi*vf1*(2 + M + m1))*(-1 + M + m1)^2 - ...
       exp(1i*2*pi*vf1*(M + m1))*(M + m1)^2 + ...
       exp(1i*2*pi*vf1*(1 + m1))*(1 - 2*(-1 + m1)*m1) + ... 
       exp(1i*2*pi*vf1*(1 + M + m1))* ...
        (-1 + 2*M^2 + 2*(-1 + m1)*m1 + M*(-2 + 4*m1)))*pi^2)./ ...
       (-1 + exp(1i*2*pi*vf1)).^3;
   
    v(~ind) = -(2*pi^2/3)*M*(1+6*m1^2+6*m1*(-1+M)-3*M+2*M^2);
    
end

function st = Report(p,vfEst)

st = [];
st.Freqs = vfEst;
st.Amps = Correlation(p.cfd.icf,'ampEstimates',vfEst);

function FreqsDevBounds = CRBBound(vf,va,M,nVar)

%Computes the Cramer-Rao bound of the frequencies.

va = va(:);
vf = vf(:);

Raa = RegGeoSum(M,vf,vf,0,0);
Rad = RegGeoSum(M,vf,vf,0,1);
Rdd = RegGeoSum(M,vf,vf,1,1);

FreqsDevBounds = sqrt((nVar/2)*diag(inv(real((conj(va)*va.') .* (Rdd-Rad'*(Raa\Rad))))));

function st = MLEstimate(y,nVar,PFA)

p = Preprocessing(y,nVar,PFA);

st = MethodSerial1(p);

st.CRBBoundsOnFreqDevs = CRBBound(st.Freqs,st.Amps,length(y),nVar);

Phaser audio effect

Gabriel Rivas ● July 4, 2011●2 comments ● Coded in C 
/*
Phaser audio effect:

Xin                                                                       Yout
-------------------------------[dir_mix]--------------------->[+]--------->
      |                                                        ^  
      |                                                        |
      |-->[VNS1]-->[VNS2]-->[VNS3]...-->[VNSN]-->[pha_mix]------      

            ^        ^        ^           ^ 
            |        |        |           |
            |--------|--------|---...------
            |
          [LFO]         

VNS = Variable notch stage

*/

#include "br_iir.h"
#include "Phaser.h"

/*This defines the phaser stages
that is the number of variable notch blocks
*/
#define PH_STAGES 20

static short center_freq;    /*Center frequency counter*/
static short samp_freq;      /*Sampling frequency*/
static short counter;        /*Smaple counter*/
static short counter_limit;  /*Smaple counter limit*/
static short control;        /*LFO Control*/
static short max_freq;       /*Maximum notch center frequency*/
static short min_freq;       /*Minimum notch center frequency*/ 
static double pha_mix;       /*Filtered signal mix*/
static short f_step;         /*Sweep frequency step*/ 
static double dir_mix;       /*Direct signal mix*/
static struct br_filter H[PH_STAGES]; /*Array of notch filters stages*/

/*
This funtion initializes the phaser control variables
and the variable notch filter coefficients array
*/
void Phaser_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,double pha_mixume,short freq_step, double dmix) {
    /*Initialize notch filter coefficients set array*/
    br_iir_init(sampling,gainfactor,Q,freq_step, minf);
   
    /*Initializes the phaser control variables*/
    center_freq = 0;
    samp_freq = sampling;
    counter = effect_rate;
    control = 0;
    counter_limit = effect_rate;

    /*Convert frequencies to integer indexes*/
    min_freq = 0;
    max_freq = (maxf - minf)/freq_step;

    pha_mix = pha_mixume;
    f_step = freq_step;
    dir_mix = dmix;
}

/*
This function does the actual phasing processing
1. It takes the input sample and pass it trough the
cascaded notch filter stages
2. It takes tha output of the cascaded notch filters
and scales it, scales the input sample and generate
the output effect sample.
*/
double Phaser_process(double xin) {
    double yout;
    int i;

    yout = br_iir_filter(xin,&H[0]);

    for(i = 1; i < PH_STAGES; i++) {
        yout = br_iir_filter(yout,&H[i]);
    }
     
    yout = dir_mix*xin + pha_mix*yout;
      
    return yout;
}

/*
This function makes vary the center notch frequency
in all the cascaded notch filter stages by a simulated
triangle wave LFO that goes up and down
*/
void Phaser_sweep(void) {
    int i;

    if (!--counter) {
        if (!control) {
            center_freq+=f_step;

            if (center_freq > max_freq) {
                control = 1;
            }
        }
        else if (control) {
            center_freq-=f_step;

            if (center_freq == min_freq) {               
                control = 0;
            }
        }
        for(i = 0; i < PH_STAGES; i++) {
            br_iir_setup(&H[i],center_freq);
        }       
        counter = counter_limit;
    }
}

/************

Phaser.h

***********/

#ifndef __PHASER_H__
#define __PHASER_H__

extern void Phaser_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,double mix_volume,short freq_step, double dmix);
extern double Phaser_process(double xin);
extern void Phaser_sweep(void);

#endif

Variable 2nd order band pass IIR filter

Gabriel Rivas ● July 4, 2011●2 comments ● Coded in C 
#include "bp_iir.h"
#include <math.h>
#define BP_MAX_COEFS 120
#define PI 3.1415926

/*This is an array of the filter parameters object
defined in the br_iir.h file*/
static struct bp_coeffs bp_coeff_arr[BP_MAX_COEFS];

/*This initialization function will create the
band pass filter coefficients array, you have to specify:
fsfilt = Sampling Frequency
gb     = Gain at cut frequencies
Q      = Q factor, Higher Q gives narrower band
fstep  = Frequency step to increase center frequencies
in the array
fmin   = Minimum frequency for the range of center   frequencies
*/
void bp_iir_init(double fsfilt,double gb,double Q,short fstep, short fmin) {
      int i;
      double damp;
      double wo;
      
       damp = gb/sqrt(1 - pow(gb,2));

	for (i=0;i<BP_MAX_COEFS;i++) {
            wo = 2*PI*(fstep*i + fmin)/fsfilt;
		bp_coeff_arr[i].e = 1/(1 + damp*tan(wo/(Q*2)));
		bp_coeff_arr[i].p = cos(wo);
		bp_coeff_arr[i].d[0] = (1-bp_coeff_arr[i].e);
		bp_coeff_arr[i].d[1] = 2*bp_coeff_arr[i].e*bp_coeff_arr[i].p;
		bp_coeff_arr[i].d[2] = (2*bp_coeff_arr[i].e-1);
	}
}

/*This function loads a given set of band pass filter coefficients acording to a center frequency index
into a band pass filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
void bp_iir_setup(struct bp_filter * H,int ind) {
	H->e = bp_coeff_arr[ind].e;
	H->p = bp_coeff_arr[ind].p;
	H->d[0] = bp_coeff_arr[ind].d[0];
	H->d[1] = bp_coeff_arr[ind].d[1];
	H->d[2] = bp_coeff_arr[ind].d[2];
}

/*This function loads a given set of band pass filter coefficients acording to a center frequency index
into a band pass filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
double bp_iir_filter(double yin,struct bp_filter * H) {
	double yout;

	H->x[0] =  H->x[1]; 
	H->x[1] =  H->x[2]; 
	H->x[2] = yin;
	
	H->y[0] =  H->y[1]; 
	H->y[1] =  H->y[2]; 

	H->y[2] = H->d[0]* H->x[2] - H->d[0]* H->x[0] + (H->d[1]* H->y[1]) - H->d[2]* H->y[0];
	
	yout =  H->y[2];

	return yout;
}

/******************
bp_iir.h
**********************/

#ifndef __BP_IIR_H__
#define __BP_IIR_H__

struct bp_coeffs{
	double e;
	double p;
	double d[3];
};

struct bp_filter{
	double e;
	double p;
	double d[3];
	double x[3];
	double y[3];
};

extern void bp_iir_init(double fsfilt,double gb,double Q,short fstep, short fmin);
extern void bp_iir_setup(struct bp_filter * H,int index);
extern double bp_iir_filter(double yin,struct bp_filter * H);

#endif

Digital Comb Filter

Jeff T ● June 21, 2011 ● Coded in Matlab 
% Usage:     [B,A] = COMB(order, scalar);
%
%             ORDER is the number of samples delayed prior to add
%             SCALAR is the coefficient that will be applied to
%                the delayed signal path at the final summation block.
%
% Note, there are two types of comb filters.  A DC-blocker and a DC-passer.
% To get a DC-Blocker (tooth at DC), pass in a -1 for the scalar.  
% To get a DC-Passer (+6dB at DC), pass in a +1 for the scalar.
%
% By default, if the scalar is not passed, a DC-Passer is assumed.
%
% Author: sparafucile17 03/16/04

% Validate that the proper argument count was supplied
error(nargchk(1, 2, nargin));

% Use scalar is passed as an argument, otherwise assume scalar=1;
if (nargin == 1)
    scalar = 1; 
else
    scalar = varargin{1};
end

% Input has zeros to simulate a single delay of N samples
a = [ 1 ];
b = [ 1 zeros(1, N-1) scalar*1];

Equal Loudness Curves (ISO226)

Jeff T ● June 21, 2011●1 comment ● Coded in Matlab 
function [spl, freq] = iso226(phon);
%
% Generates an Equal Loudness Contour as described in ISO 226
%
% Usage:  [SPL FREQ] = ISO226(PHON);
% 
%         PHON is the phon value in dB SPL that you want the equal
%           loudness curve to represent. (1phon = 1dB @ 1kHz)
%         SPL is the Sound Pressure Level amplitude returned for
%           each of the 29 frequencies evaluated by ISO226.
%         FREQ is the returned vector of frequencies that ISO226
%           evaluates to generate the contour.
%
% Desc:   This function will return the equal loudness contour for
%         your desired phon level.  The frequencies evaulated in this
%         function only span from 20Hz - 12.5kHz, and only 29 selective
%         frequencies are covered.  This is the limitation of the ISO
%         standard.
%
%         In addition the valid phon range should be 0 - 90 dB SPL.
%         Values outside this range do not have experimental values
%         and their contours should be treated as inaccurate.
%
%         If more samples are required you should be able to easily
%         interpolate these values using spline().
%
% Author: sparafucile17 03/01/05

%                /---------------------------------------\
%%%%%%%%%%%%%%%%%          TABLES FROM ISO226             %%%%%%%%%%%%%%%%%
%                \---------------------------------------/
f = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 ...
     1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500];

af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
      0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
      0.243 0.242 0.242 0.245 0.254 0.271 0.301];

Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
       -2.0  -1.1  -0.4   0.0   0.3   0.5   0.0 -2.7 -4.1 -1.0  1.7 ...
        2.5   1.2  -2.1  -7.1 -11.2 -10.7  -3.1];

Tf = [ 78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5  22.1  17.9  14.4 ...
       11.4   8.6   6.2   4.4   3.0   2.2   2.4   3.5   1.7  -1.3  -4.2 ...
       -6.0  -5.4  -1.5   6.0  12.6  13.9  12.3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

%Error Trapping
if((phon < 0) | (phon > 90))
    disp('Phon value out of bounds!')
    spl = 0;
    freq = 0;
else
    %Setup user-defined values for equation
    Ln = phon;

    %Deriving sound pressure level from loudness level (iso226 sect 4.1)
    Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
    Lp=((10./af).*log10(Af)) - Lu + 94;

    %Return user data
    spl = Lp;  
    freq = f;
end

Variable 2nd order Notch IIR filter

Gabriel Rivas ● June 21, 2011●2 comments ● Coded in C 
#include "br_iir.h"
#include <math.h>

#define BR_MAX_COEFS 120
#define PI 3.1415926

/*This is an array of the filter parameters object
defined in the br_iir.h file*/
static struct br_coeffs br_coeff_arr[BR_MAX_COEFS];

/*This initialization function will create the
notch filter coefficients array, you have to specify:
fsfilt = Sampling Frequency
gb     = Gain at cut frequencies
Q      = Q factor, Higher Q gives narrower band
fstep  = Frequency step to increase center frequencies
in the array
fmin   = Minimum frequency for the range of center   frequencies
*/
void br_iir_init(double fsfilt,double gb,double Q,double fstep, short fmin) {
	int i;
      double damp;
      double wo;
      
      damp = sqrt(1 - pow(gb,2))/gb;

	for (i=0;i<BR_MAX_COEFS;i++) {
                wo = 2*PI*(fstep*i + fmin)/fsfilt;
		br_coeff_arr[i].e = 1/(1 + damp*tan(wo/(Q*2)));
		br_coeff_arr[i].p = cos(wo);
		br_coeff_arr[i].d[0] = br_coeff_arr[i].e;
		br_coeff_arr[i].d[1] = 2*br_coeff_arr[i].e*br_coeff_arr[i].p;
		br_coeff_arr[i].d[2] = (2*br_coeff_arr[i].e-1);
	}
}

/*This function loads a given set of notch filter coefficients acording to a center frequency index
into a notch filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
void br_iir_setup(struct br_filter * H,int ind) {
	H->e = br_coeff_arr[ind].e;
	H->p = br_coeff_arr[ind].p;
	H->d[0] = br_coeff_arr[ind].d[0];
	H->d[1] = br_coeff_arr[ind].d[1];
	H->d[2] = br_coeff_arr[ind].d[2];
}

/*
This function does the actual notch filter processing
yin = input sample
H   = filter object 
*/
double br_iir_filter(double yin,struct br_filter * H) {
	double yout;

	H->x[0] = H->x[1]; 
	H->x[1] = H->x[2]; 
	H->x[2] = yin;
	
	H->y[0] = H->y[1]; 
	H->y[1] = H->y[2]; 

	H->y[2] = H->d[0]*H->x[2] - H->d[1]*H->x[1] + H->d[0]*H->x[0] + H->d[1]*H->y[1] - H->d[2]*H->y[0];
	
	yout = H->y[2];
	return yout;
}

/******************/
#ifndef __BR_IIR_H__
#define __BR_IIR_H__

/*
Notch filter coefficients object
*/
struct br_coeffs {
	double e;
	double p;
	double d[3];
};

/*
Notch filter object
*/
struct br_filter {
	double e;
	double p;
	double d[3];
	double x[3];
	double y[3];
};

extern void br_iir_init(double fsfilt,double gb,double Q,double fstep, short fmin);
extern void br_iir_setup(struct br_filter * H,int index);
extern double br_iir_filter(double yin,struct br_filter * H);

#endif

Farrow 2-D Image Interpolation

Jesus Selva ● June 15, 2011 ● Coded in Matlab 
function Demo

%Author: J. Selva.
%Date: June 2011.

close all
format short g
format compact

T = 1; %Sampling period in both dimensions.
B = 1/(1.223*T); %This is the two-sided bandwidth in both dimensions. It is a typical ...
                 %value in SAR image co-registration.
P = 18; %Kernel semi-length.
Q = 10; %Polynomial order in Farrow structure.

%The next call to "ChebApproxPolys" computes all Farrow coefficients.
t0 = -P*T-T/2; %The parameters t0, Nt, 2*P+1 specify the 2*P+1 intervals in which Knab's ...
               %pulse is to be approximated by polynomials. 
Nt = 2*P+1;
Dt = T;
IOut = [-T/2,T/2];
  
%Perform Chebyshev interpolation. The parameters T, B and  P are passed to KnabAPPulse.
K = ChebApproxPolys('KnabAPPulse',t0,Nt,Dt,IOut,40,Q,T,B,P).';

%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = 100;

x0 = (0:Lx-1)';
x0 = x0(:,ones(1,Lx));

y0 = x0.';
x0 = x0(:);
y0 = y0(:);

disp('Computing image in a regular grid...')

tic;
Image = reshape(TestImage5(x0,y0,B),[Lx,Ly]);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

clear x0 y0

%Specify NP points on the image at random.
NP = 10000;

x = P*T+rand(1,NP)*(Lx-2*P)*T;
y = P*T+rand(1,NP)*(Ly-2*P)*T;

%Evaluate image directly.
disp('Computing image at selected positions exactly...')
tic;
vRef = TestImage5(x,y,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

%Interpolate the image using the 2D-Farrow structure. 
disp('Computing image at selected positions using Farrow 2-D...')
tic;
v = Farrow2DInterp(Image,K,x,y);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

%Show the error. 
disp(['Maximum interpolation error: ', ...
      num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])

disp(' ')
disp('Display image and interpolation error:')
disp(' ')

%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = Lx;

x1 = (P*T:0.5*T:(Lx-P)*T)';
x1 = x1(:,ones(1,length(x1)));

y1 = x1.';
x1 = x1(:);
y1 = y1(:);

tic;
disp('Computing image in over-sampled grid ...')
vRef =  TestImage5(x1,y1,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

tic;
disp('Computing signal in over-sampled grid using 2-D Farrow ...')
v = Farrow2DInterp(Image,K,x1,y1);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

x1 = (P*T:0.5*T:(Lx-P)*T)';
y1 = x1.';

disp(['Maximum interpolation error: ', ...
      num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])

vRef = reshape(vRef,[length(x1),length(y1)]);
v = reshape(v,[length(x1),length(y1)]);

mesh(x1,y1,abs(vRef))

xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Image''s abs. value')
title('Test Image (abs. value)')

figure

mesh(x1,y1,20*log10(abs(v-vRef)/max(abs(vRef(:)))))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Interp. error (dB)')
title('Interpolation error (dB)')

function P = ChebApproxPolys(Func,t0,Nt,Dt,IOut,NPoints,Q,varargin)

%This function constructs Chebyshev interpolation polynomials in a set of consecutive ...
%intervals.
%
% Func is the function to be interpolated.
% t0 is the first instant of the first interval.
% Nt is the number of intervals.
% IOut is the interval in which the output polynomials will be evaluated.
% Dt is the interval width. The intervals are [t0,t0+Dt], [t0+Dt,t0+2*Dt],....
%    [t0+(Nt-1)*Dt,t0+Nt*Dt].
% NPoints is the number of points in which Func is evaluated, for each interval.
% Q is the output polynomial order.
% varargin contains extra parameters which are passed to varargin.

%Construct matrix that maps function values to polynomial coefficients. 
%The algorithm is that in section 9.8 of
% 
% W. H. Press et al., Numerical Recipes in C. Cambridge University
% Press, 1997.

M = MapPolyMatrix([-1,1],IOut,Q) * ChebCoefMatrix(Q) * ...
    BlockChebInvMatrix(NPoints,Q);

% Get Chebyshev roots.
r = ChebRoots(NPoints);
r = r(:);

P = zeros(Q,Nt);

%Compute polynomial for each interval.
for k = 0:Nt-1,
  v1 = feval(Func,MapIaToIb(r,[-1,1],t0+[k,k+1]*Dt),varargin{:});
  P(:,k+1) = double(M * v1(:));
end

function M = MapPolyMatrix(IIn,IOut,Q)

%For a polynomial with Q coefficients which is evaluated in interval IIn, this function ...
%gives the matrix that transforms its coefficients, so that it can be evaluated in the ...
%associated points in interval IOut. 

ayx = (IIn(1)-IIn(2))/(IOut(1)-IOut(2));
byx = IIn(1)-ayx*IOut(1);

if byx == 0
  M = diag(ayx.^(0:Q-1));
else
  M = zeros(Q);
  for q = 0:Q-1
    for r = 0:q
      M(r+1,q+1) = nchoosek(q,r) * ayx^r * byx^(q-r);
    end
  end
end

function M = ChebCoefMatrix(N)

M = zeros(N);

M(1) = 1;

if N > 1
  M(2,2) = 1;
end

for n = 3:N,
  M(2:n,n) = 2*M(1:n-1,n-1);
  M(1:n-2,n) = M(1:n-2,n) - M(1:n-2,n-2);
end

function M = BlockChebInvMatrix(N,Q)

r = ChebRoots(N);
r= r(:);

M = zeros(N,Q);

if Q >=1
  M(:,1) = 1;
end

if Q >= 2
  M(:,2) = r;
end

for q = 3:Q,
  M(:,q) = 2 * r .* M(:,q-1) - M(:,q-2);
end

M = M.'/N;
M(2:end,:) = M(2:end,:) * 2;

function R = ChebRoots(N)

R = cos(pi*((1:N)-1/2)/N);

function y = MapIaToIb(x,Ia,Ib)

%Linear transformation that maps interval Ia onto interval Ib.

y = (Ib(1)+Ib(2))/2 + (x-(Ia(1)+Ia(2))/2) * ((Ib(2)-Ib(1))/(Ia(2)-Ia(1))); 

function [g,BoundAP] = KnabAPPulse(t,T,B,P)

%Knab interpolation pulse in 
%
%  J. J. Knab, 'Interpolation of band-limited functions using the Approximate
%  Prolate series', IEEE Transactions on Information Theory, vol. IT-25, 
%  no. 6, pp. 717-720, Nov 1979. 
%
%Grid points -P*T:T:P*T.
%Interpolation interval: [-T/2,T/2].

g = sinc(t/T) .* sinc((1/T-B)*sqrt(t.^2-(P*T)^2)) / sinc(j*(1-B*T)*P);

if nargout == 2
  BoundAP = 1/sinh(P*pi*(1-B*T));
end

function v = TestImage5(x,y,B)

%Test image formed by a sum of random sincs.

Seed = rand('state');

rand('state',362436069);

L = 100; %Image size in both dimensions.

Np = 5000;
p = rand(Np,2)*L;
A = exp(j*2*pi*rand(Np,1));

rand('state',Seed);

v = zeros(size(x));

for k = 1:length(x(:))
  v(k) = A(:).' * (sinc(B*(x(k)-p(:,1))) .* sinc(B*(y(k)-p(:,2))));
end

function vx = Farrow2DInterp(Image,Kernel,x,y)

%Farrow interpolation in two dimensions using interlaced FFTs.

vxSize = size(x);

[LK,NK] = size(Kernel);
[Lx,Ly] = size(Image);

nx = round(x);
ny = round(y);
ux = x-nx;
uy = y-ny;

np = 1+nx+ny*Lx;
np = np(:);

clear x y nx ny

P = (LK-1)/2;

NFFTx = 2^ceil(log2(Lx+LK-1));
NFFTy = 2^ceil(log2(Ly+LK-1));

KernelFx = fft(Kernel,NFFTx);
KernelFy = fft(Kernel,NFFTy);

ImageFx = fft(Image,NFFTx);

vx = zeros(size(np))+j*zeros(size(np));
for kx = NK:-1:1,
 
  ImageCx = ifft(ImageFx .* KernelFx(:,kx(ones(1,Ly))),NFFTx); 
  ImageCx = ImageCx(1+P:Lx+LK-1-P,:);
  
  ImageCxFy = fft(ImageCx,NFFTy,2);
  
  vy = zeros(size(np))+j*zeros(size(np));
  for ky = NK:-1:1,
    ImageCxCy = ifft(ImageCxFy .* KernelFy(:,ky(ones(1,Lx))).',NFFTy,2);
    ImageCxCy = ImageCxCy(:,1+P:Ly+LK-1-P);
  
    vy(:) = vy(:) .* uy(:) + ImageCxCy(np);
  end
  
  vx(:) = vx(:) .* ux(:) + vy(:);
end

vx = reshape(vx,vxSize);

Shelving Filter Design

Jeff T ● May 29, 2011●1 comment ● Coded in Matlab 
function [b, a]  = shelving(G, fc, fs, Q, type)

%
% Derive coefficients for a shelving filter with a given amplitude and
% cutoff frequency.  All coefficients are calculated as described in 
% Zolzer's DAFX book (p. 50 -55).  
%
% Usage:     [B,A] = shelving(G, Fc, Fs, Q, type);
%
%            G is the logrithmic gain (in dB)
%            FC is the center frequency
%            Fs is the sampling rate
%            Q adjusts the slope be replacing the sqrt(2) term
%            type is a character string defining filter type
%                 Choices are: 'Base_Shelf' or 'Treble_Shelf'
%
% Author:    sparafucile17 08/22/05
%

%Error Check
if((strcmp(type,'Base_Shelf') ~= 1) && (strcmp(type,'Treble_Shelf') ~= 1))
    error(['Unsupported Filter Type: ' type]);
end

K = tan((pi * fc)/fs);
V0 = 10^(G/20);
root2 = 1/Q; %sqrt(2)

%Invert gain if a cut
if(V0 < 1)
    V0 = 1/V0;
end

%%%%%%%%%%%%%%%%%%%%
%    BASE BOOST
%%%%%%%%%%%%%%%%%%%%
if(( G > 0 ) & (strcmp(type,'Base_Shelf')))
   
    b0 = (1 + sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
    b1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*K + K^2);
    b2 = (1 - sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
    a1 =                (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
    a2 =             (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%    BASE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Base_Shelf')))
    
    b0 =             (1 + root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
    b1 =                (2 * (K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
    b2 =             (1 - root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
    a1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
    a2 = (1 - root2*sqrt(V0)*K + V0*K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE BOOST
%%%%%%%%%%%%%%%%%%%%
elseif (( G > 0 ) & (strcmp(type,'Treble_Shelf')))

    b0 = (V0 + root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
    b1 =             (2 * (K^2 - V0) ) / (1 + root2*K + K^2);
    b2 = (V0 - root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
    a1 =              (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
    a2 =           (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE CUT
%%%%%%%%%%%%%%%%%%%%

elseif (( G < 0 ) & (strcmp(type,'Treble_Shelf')))

    b0 =               (1 + root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
    b1 =                  (2 * (K^2 - 1) ) / (V0 + root2*sqrt(V0)*K + K^2);
    b2 =               (1 - root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
    a1 =             (2 * ((K^2)/V0 - 1) ) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
    a2 = (1 - root2/sqrt(V0)*K + (K^2)/V0) / (1 + root2/sqrt(V0)*K + (K^2)/V0);

%%%%%%%%%%%%%%%%%%%%
%   All-Pass
%%%%%%%%%%%%%%%%%%%%
else
    b0 = V0;
    b1 = 0;
    b2 = 0;
    a1 = 0;
    a2 = 0;
end

%return values
a = [  1, a1, a2];
b = [ b0, b1, b2];