[Chapter 6][Table of Contents][Appendix B]

Appendix A - MATLAB code for test tone synthesis, selection, and all-pass filtering

%--------------------------------------------------
% MATLAB code for test tone synthesis and filtering
% by Daisuke Koya
%--------------------------------------------------

clear all;
h = get(0,'Children'); % get handle of current plots
set(h,'color',[1 1 1]); % set background of plot to white

% test tone generation
fs = 44100; % sampling frequency of test tone in Hz
time = (0:1/fs:1);% 1 sec. time duration of sine wave, T = 1/44100Hz

%--------
% impulse
%--------

imp = zeros(44100,1);
imp(1) = 1;
% f0 = 3500 is used for the AP filter f0
% (assumed to be most sensitive by the Fletcher-Munson curves)

%-------------------
% 70Hz sawtooth wave
%-------------------

%st70 = sawtooth(2*pi*70*time); % 70Hz sawtooth wave
%st70 = 0.35.*st70;
% reduce amplitude so that filtered signal
% is not clipped beyond -1, +1

%---------------------
% 3.5kHz sawtooth wave
%---------------------

%st3500 = sawtooth(2*pi*3500*time); % 3.5kHz sawtooth wave
%st3500 = 0.6.*st3500;
% reduce amplitude so that filtered signal
% is not clipped beyond -1, +1

%--------------------
% 10kHz sawtooth wave
%--------------------

%st10000 = sawtooth(2*pi*10000*time); % 10kHz sawtooth wave
%st10000 = 0.6.*st10000;
% reduce amplitude so that filtered signal
% is not clipped beyond -1, +1

% note: need to use f0=12100Hz for AP filter formulation
% since approximation has larger frequency errors at higher frequencies

% note: tmax approximation for 10kHz sawtooth waveform
% fails, therefore,
% for tmax = 0.001 mS, use tmax = 0.000285
% for tmax = 0.002 mS, use tmax = 0.000575
% for tmax = 0.004 mS, use tmax = 0.00115
% for tmax = 0.008 mS, use tmax = 0.00235

%----------------------------------------------------
% University of Miami Percussion Ensemble
% use f0 ~ 150Hz for APF formulation since spectogram
% has significant frequency content in this area
%----------------------------------------------------

%perc = wavread('perc.wav');
%perc = 0.8.*perc;
% reduce amplitude so that filtered signal
% is not clipped beyond -1, +1

%----------------------------------------------------
% M-Pact (Jazz vocal group)
% use f0 ~ 160Hz for APF formulation since spectogram
% has significant frequency content in this area
%----------------------------------------------------

%mpact = wavread('mpact.wav');

%--------------------
% test tone selection
%--------------------

x = imp; % select signal from above to be used

% !! check f0 of APF to make sure it corresponds to above signal !!

%------------------------------------------------------------
% Sampling frequency for 2nd-order allpass filter formulation
%------------------------------------------------------------

fs = 44100; % sampling frequency in Hz
T = 1/44100; % sampling period in seconds

%--------------------------
% Allpass filter parameters
%--------------------------

tmax = 0.004 % peak delay in mS of AP
f0 = 3500; % resonant frequency in Hz of AP

%---------------------------
% Allpass filter formulation
%---------------------------

Q = 0.5*tmax*pi*f0 % Q of AP transfer function
w0 = 2*pi*f0; % resonant frequency in radians
sn = 1/w0; % normalized frequency sn
numAP = [sn^2 -sn/Q 1]; % numerator of allpass filter
denAP = [sn^2 +sn/Q 1]; % denominator of allpass filter

f = (0:20:22050); % frequency vector, to 1/2 Nyquist
w = 2*pi*f; % frequency vector in radians
H = freqs(numAP,denAP,w); % complex frequency response
magH = 20*log10(abs(H));% magnitude of H

phaserad = angle(H); % phase of complex frequency response (rad.)
phasedeg = phaserad*180/pi; % convert to degrees
gd = -diff(unwrap(phaserad))./diff(w); % calculate group delay

%---------------------------------
% Analog to Digital Transformation
%---------------------------------

[numAPd,denAPd] = BILINEAR(numAP,denAP,fs);%Bilinear transformation
[Gd,F] = GRPDELAY(numAPd,denAPd,f,fs);
Gdsec = Gd*T;
[H,F] = freqz(numAPd,denAPd,f,fs);

%----------------------------------
% Filter signal with allpass filter
%----------------------------------

y = filter(numAPd, denAPd, x);

%------------
% Plot result
%------------

figure(1)
semilogx(f,Gdsec)
axis([0 22050 0 0.01])
grid
title('Digital Group Delay Response of AllPass Filter');
xlabel('frequency (Hz)'),ylabel('Group Delay (seconds)')

tsignal = ( 0:(length(x)-1) ).*T ; % convert # samples to time in seconds

figure(2)
plot(tsignal,x,'-') % snapshot of steady-state response
axis([0.000 0.05 -1 1]);
hold on
plot(tsignal,y,':')
hold off
title('Unfiltered (solid) and filtered (dotted) signals');
xlabel('time (seconds)'),ylabel('amplitude')


[Chapter 6][Table of Contents][Appendix B]