Sign in

username:

password:



Not a member?

Search Online Books



Search tips

Free Online Books

Ads

Chapters

Spectral Audio Signal Processing
    FIR Digital Filter Design
       FIR Fractional Delay Filter Design
          Example 1: Fractionally Iterated Delay

Chapter Contents:

Search Spectral Audio Signal Processing

  

Book Index | Global Index

Would you like to be notified by email when Julius Orion Smith III publishes a new entry into his blog?

  

Example 1: Fractionally Iterated Delay

Figure B.25a shows the phase response of $ z^{-4.5}$ when computed using the matlab code shown in Fig.B.26. (The plot utility is listed separately in Fig.B.27.) The fractional delay of 4.5 samples is decomposed as $ 3\cdot 1.5$ so that the transfer function can be computed as $ (z^{-3})^{1.5}$, i.e., as a fractionally iterated convolution of $ \delta(n-3)$ with itself.

Figure: Phase response corresponding to a $ 4.5$ sample delay, computed by a variety of methods in Matlab. Only the last case (d) is correct.
\includegraphics[width=\textwidth ]{eps/phasewrapz}

Figure B.26: Simple matlab program for computing the frequency response of a fractional delay FIR filter.
  
power = 1.5;      % z^{-3} will be raised to this power

N = 64;           % FFT size (power of 2)
fs = 1;           % Sampling frequency
T = 1/fs;         % Sampling period
f = (0:N-1)/N;    % Frequency axis
w = 2*pi*f;       % rad/sec

S = exp(-j*3*w*T); % Starting spectrum (3-sample delay)
Stp = S .^ power;  % Go for it
Stpp = angle(Stp);
plottitle = sprintf(['angle(S^{%0.1f}), ',...
      'where S = exp(-j*3*w*T)',...
      '= spectrum of 3-sample delay'],power);
phasesubplot(Stpp,plottitle,1,4);

Figure: Plot utility for matlab program Fig.B.26.
  
function [ph] = phasesubplot(phase,ptitle,iplot,nplots)
% PHASESUBPLOT - Plot phase data in subplot iplot/nplots
%
% Phase data assumed to be in "FFT buffer format", \ie,
% ordered from bin 0 to N-1 where N = FFT size.

subplot(nplots,1,iplot);
N = length(phase);
f = (0:N-1)/N;  % Frequency axis
ph = stem(f-0.5,fftshift(phase),'k'); 
pic = chop(pi,3); % pi to three digits
% place Y grid on multiples of pi/2:
set(gca,'YTick',[-pic -pic/2 0 pic/2 pic]);
set(gca,'XTick',-0.5:0.1:0.5);
grid on;
title(ptitle);
axis([-0.5 0.5, -4, 4]);
if (iplot == nplots)
  xlabel('Normalized Frequency (cycles per sample)');
end
ylabel('Phase (rad)');
if (nplots>1)
  yl=get(gca,'YLim'); 
  th=text(-0.6,yl(2),['(','a'+iplot-1,')']);
  set(th,'FontWeight','bold');
end

We see that everything looks good between $ \omega = -\pi/3$ and $ \omega = \pi/3$ ( $ \omega\in[-1.0472,1.0472]$), but outside of these limits the phase response jumps by $ \pi$ radians. (For reference, the correct phase response is shown in Fig.B.25d.) A phase jump of $ \pi$ radians corresponds to a sign inversion. Thus, the entire high-frequency band above $ \omega = \pi/3$ has been negated! We are trying to design a fractional delay filter here, not a ``high-frequency inverting fractional delay filter''! Something went wrong.

To analyze this failure, let's look analytically at the expected phase response. The phase response of a unit delay $ z^{-1}$ is

$\displaystyle \Theta(\omega) \isdef \angle e^{-j\omega T} = -\omega T, \quad -\pi \leq \omega T < \pi $

which ranges from $ -\pi$ to $ \pi$. Similarly, the phase response of $ z^{-3}$ is

$\displaystyle \Theta_3(\omega) \isdef \angle e^{-j3\omega T} = -3\omega T, \quad -\pi \leq \omega T < \pi $

which ranges from $ -3\pi$ to $ 3\pi$. However, this an unwrapped phase response. While this is naturally obtained by hand analysis, it is not what is typically computed in software, such as by the following matlab expression:
angle(exp(-j*3*omega*T))
It is more typical that the computed phase response will be wrapped into the range $ -\pi$ to $ \pi$, yielding

$\displaystyle \Theta_3(\omega) \isdef \angle e^{-j3\omega T} = \left\{\begin{ar... ... -3\omega T - 2\pi, & -\pi \leq \omega T < -\frac{\pi}{3}. \end{array}\right. $

I.e., the phase of $ z^{-3}$ is wrapped to lie between $ -\pi$ and $ \pi$. This wrapped phase is shown in Fig.B.28.

Figure: Phase response of a $ 3$-sample delay line, computed as angle(exp(-j*3*omega*T)) in Matlab.
\includegraphics[width=3.5in,height=2.25in]{eps/phasewrapSp}

Figure B.28 depicts the phase spectrum of S = exp(-j*3*w*T), while Fig.B.25a shows the phase spectrum of S^1.5. Consider now what happens when each sample of the complex array S is raised to the power 1.5 in matlab by the expression S .^ 1.5. Let $ z = x + jy$ denote an arbitrary complex number. Then we need to know how matlab computes $ z^{1.5}$. This can be thought of as fractionally iterated multiplication. In general, we are asking how to compute $ z^\alpha$, where $ \alpha $ is a real exponent and $ z$ is any complex number.

One method for computing exponentiation by a real number, $ z^\alpha$, is to approximate $ \alpha $ by a rational number $ N/M \approx \alpha$ [243]. We try out this approach in Example 2 of the next section.

Another method for computing exponentiation $ z^\alpha$ is by means of logarithms. Expressing $ z$ in polar notation as $ z=r\,e^{j\theta}$, the natural logarithm of $ z^\alpha$ can be written as

$\displaystyle \ln\left(z^\alpha\right) = \alpha \ln(z) = \alpha \ln\left(r\,e^{j\theta}\right) = \alpha \ln(r) + j\alpha\theta. $

We may therefore compute $ z^\alpha$ as

$\displaystyle z^\alpha \isdef e^{\alpha \ln(r)} e^{j\alpha\theta}. \protect$ (B.6)

From Fig.B.25a, we can deduce that Matlab computes $ z^\alpha$ by the logarithm method. This is because Fig.B.25a is clearly

$\displaystyle \Theta_{4.5}(\omega) = \hbox{\sc wrap}_{\pm\pi}(1.5\cdot\hbox{\sc wrap}_{\pm\pi}(3\omega T)) \protect$ (B.7)

where $ \hbox{\sc wrap}_{\pm\pi}(x)$ denotes the operation wrapping $ x$ into the range $ [-\pi,\pi)$ by adding or subtracting multiples of $ 2\pi$. Thus, we see that Matlab effectively multiplies the wrapped phase $ \Theta_3(\omega) = -\hbox{\sc wrap}_{\pm\pi}(3\omega T)$ (shown in Fig.B.28) by 1.5, thereby stretching it vertically by 50%. Exponentiating the result, according to (B.6), effectively applies $ \hbox{\sc wrap}_{\pm\pi}$ again, so that a subsequent angle() call sees (B.7), as shown in Fig.B.25a.

Order a Hardcopy of Spectral Audio Signal Processing

Previous: FIR Fractional Delay Filter Design
Next: Example 2: Rational Approximation
written by Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.


Comments

No comments yet for this page


Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )