Matlab Demos

I was the teaching assistant for Julius Smith's digital signal processing courses, Music 320A & B in fall 2014 and winter 2015. Much of the class time was spent on some matlab demos that I had prepared. Here they are!

I have included explanations and wikipedia links for subjects that deserve a bit more explanation.

---

Vector Projections

This demo takes two vectors and projects one onto the other. This demo is in preparation of the DFT (Discrete Fourier Transform) derivation in which some signal is projected onto the sinusoidal basis functions. The code below has one moving vector projected onto a stationary vector. The projection is shown as well as the difference between the projection and the moving vector. The projected component is always in the direction of y, and the difference is always perpendicular to the projection.

% Vector projection
for theta = 0:0.01:2*pi % Go around the circle once
    x = 3*[cos(theta), sin(theta)]; 
    y = [1,0];
    % Project x onto y
    x_projy = dot(x,y)/dot(y,y)*y;
    hold off
    % The vector we are projecting (blue)
    plot([0 x(1)], [0 x(2)], 'b-','LineWidth',3);
    axis([-5 5 -5 5]); % Change the window size
    hold on
    % The line we are projecting onto (red)
    plot([0 y(1)], [0 y(2)], 'r--','LineWidth',3);
    
    % The removed component (x - x_projy) (black)
    plot([x(1) x_projy(1)], [x(2) x_projy(2)], 'k--'); 
    % The projected component (green)
    plot([0 x_projy(1)], [0 x_projy(2)], 'g-','LineWidth',2);
    pause(0.01)
end
        

Reconstructing a signal

Sinusoids are used to reconstruct some arbitrary signal. We start with DC and increase in frequency, filling in the spectrum. For any arbitrary signal, we see that progressively higher frequency sinusoids can be used to reconstruct. This is known as reconstruction using the Fourier series.

N = 512;     % The length of the signal
n = (0:N-1); % The indices of the signal
fs = 800;    % A sample rate
f = 1;       % A scalar for the frequency (try raising and lowering this)
x = sin(2*pi*f*0.05*n.^2/fs)+cos(2*pi*f*0.07*n.^2/fs); % Some random signal

k = (0:N-1);             % The frequency axis (bins)
omega = k/N*2-1;         % Relabel in units of rads/sample
S = exp(-j*2*pi*n'*k/N); % The sinusoids that we are projecting onto
X = x*S;                 % Multiply by the DFT matrix (the projection)

components = zeros(1,N);
X_components = zeros(1,N);
    
for this_k = 1:N/2+1
    % Split into magnitude and phase
    mag = abs(X(this_k));
    phase = angle(X(this_k));

    % Grab a positive and negative pair of frequencies
    component_contribution = mag * cos(2*pi*n*(this_k-1)/N+phase)/N;
    X_components(this_k) = X(this_k); 
    
    % DC and half the sampling rate don't have conjugate pairs
    if this_k ~= 1 && this_k ~= N/2+1
        component_contribution = component_contribution * 2; % Add more sinusoids
        X_components(N-this_k+2) = X(N-this_k+2);            % Fill in spectrum
    end
    % Add the next component into the accumulating signal
    components = components + component_contribution;
    
    % Plot the reconstruction of the time domain signal
    subplot(121)
    plot(components,'b')
    axis([0 N -2.5 2.5])
    hold on
    plot(x,'r--')
    xlabel('Samples')
    hold off
    
    % Frequency info (fill in the spectrum as the signal is reconstructed)
    subplot(122)
    plot(omega, fftshift(abs(X)),'r');
    hold on
    area(omega, fftshift(abs(X_components)));
    hold off
    xlabel('Radian Frequency')
    ylabel('Magnitude')
    title('Spectrum')
    
    pause(exp(-this_k*0.05))
    
end
        

Filtering the Messiah

Similar to the previous demo, we reconstruct Handel's Messiah using sinusoids, just a few at a time. Because we are starting with low frequencies and working our way up to higher ones, this is equivalent to low pass filtering the signal (with progressively less agressive filters). Listen to each one!

load handel %Handel's Messiah

N = length(y);
n = (0:N-1);
x = y; % Some random signal

k = (0:N-1);             % The frequency axis (bins)
omega = k/N*2-1;         % Relabel in units of rads/sample
X = fft(y);              % Perform an FFT on the signal

components = zeros(1,N);
X_components = zeros(1,N);
    
for this_k = 1:floor(N/2+1)
    % Split into magnitude and phase
    mag = abs(X(this_k));
    phase = angle(X(this_k));
    
    % Grab a positive and negative pair of frequencies
    component_contribution = mag * cos(2*pi*n*(this_k-1)/N+phase)/N;
    X_components(this_k) = X(this_k); 
    
    % DC and half the sampling rate don't have conjugate pairs
    if this_k ~= 1 && this_k ~= N/2+1
        component_contribution = component_contribution * 2; % Add more sinusoids
        X_components(N-this_k+2) = X(N-this_k+2);            % Fill in spectrum
    end
    % Add the next component into the accumulating signal
    components = components + component_contribution;
 
    %Do it in chunks so it doesn't take all day
    if (mod(this_k, round(N/10))==0)
        % Frequency info
        plot(omega, fftshift(abs(X)),'r');
        hold on
        area(omega, fftshift(abs(X_components)));
        hold off
        xlabel('Radian Frequency')
        ylabel('Magnitude')
        title('Spectrum')
        pause
        soundsc(ifft(X_components))
        
   end
    
end
        

Phase-stripping Handel

What is the importance of phase? Don't we usually look at magnitude plots? This demo zeros out the phase for the entire signal. The consequence of this is that we get impulsive behavior at time zero. The reason is becuase our reconstructing DFT cosines at every frequency are at zero phase. This is the only sample where each cosine is at a maximum at the same time. Even more interesting than that is the rest of the signal--it sounds like the entire song is happening simultaneously.

load handel %Handel's Messiah
sound(y,Fs);

N = length(y);
Y = fft(y);

% The magnitude and phase of the unprocessed spectrum
omega = (0:N-1)/N*2-1;
subplot(121)
plot(omega,abs(Y))
xlabel('Radian Frequency')
ylabel('Magnitude')
subplot(122)
plot(omega,unwrap(angle(Y)))
xlabel('Radian Frequency')
ylabel('Magnitude')
title('Spectrum')
pause

% Make a copy and zero out the phase
new_Y = abs(Y); 
new_y = ifft(new_Y);
soundsc(new_y,Fs); % Play the sound
subplot(111)
plot(ifft(new_Y))
xlabel('Samples')
ylabel('Amplitude')
pause

% Cut out the impulse and play again (rescaled)
new_y = new_y(500:end-500);
soundsc(new_y,Fs);
subplot(111)
plot(new_y)
xlabel('Samples')
ylabel('Amplitude')
        

A Graphical Convolution

This demo shows the geometric approach to a convolution. One signal slides past the other in and the amount of overlap (more accurately, the product of the signals) is computed. The overlap as well as the sliding signals are shown in the animation. How could we convert this to a correlation by only changing one line of code?

% Two signals
x = ones(1,100);
y = ones(1,100);
% Normalize the area to 1 so the signal is a good size to plot
norm_area = sum(abs(y));

% Zero padding
N = max(length(x),length(y))+300;
x = [x, zeros(1, N - length(x))];
y = [y, zeros(1, N - length(y))];

x_conv_y = zeros(1,N);

for i = 1:N
   accumulate = 0;
    y_shift = zeros(1,N);
    % Perform the convolution
    for j = 1:N
        if i-j1
            y_shift(j) = y(i-j);
            accumulate = accumulate + x(j)*y_shift(j);
        end
    end
    x_conv_y(i) = accumulate;
    
    area(x.*y_shift,'FaceColor',[1 1 0])
    hold on
    plot(x)
    % For speed, comment the legend out
    legend('overlapping area','x','y','conv(x,y)')
    plot(y_shift,'g')
    plot(x_conv_y/norm_area,'r')
    axis([0 N -1 2])
    hold off
    pause(0.001)
end
 
Filtering a real signal:
clear
clc
x = wavread('soundclip.wav');
x = x(1:5*44100,1);
Fs = 44100;
N = 100; % Change this
f = ones(1,N)/N; % Our filter
y = conv(x,f); % An easy convolution

sound(y,Fs);

omega = linspace(-1,1,8192);
plot(omega', abs(fftshift(fft(f,8192))));
        

Profiling the FFT

The time domain convolution of two length N signals can be computed using N^2 multiplication operations. What if we did an FFT? Just how much faster is the Fast Fourier Fransform?

All_N = 2.^(1:15);
t_conv = zeros(length(All_N),1); % Arrays to store the times 
t_fft  = zeros(length(All_N),1); 

count = 1;

for N = All_N
    N
    x = rand(1,N);
    y = rand(1,N);
    
    t = tic; % Start the timer
    % Convolution
    for i = 1:100
        conv(x,y); % An O(n^2), time domain convolution
    end
    t_conv(count) = toc(t)/100; %Stop the timer
    
    t = tic;
    % FFT
    for i = 1:100
        ifft(fft(x).*fft(y)); % An O(nlogn), frequency domain convolution
    end
    t_fft(count) = toc(t)/100;
    
    count = count + 1
end
% Plot the time domain convolution time
plot(All_N, t_conv);
hold on
% Plot the FFT convolution time (in seconds) in green
plot(All_N, t_fft,'g');
        

A Graphical Convolution - Making a Triangle Wave

This is nearly identical to the previous example, which shows a convolution. This shows how convolving a train of square pulses with a square pulse gives a triangle wave.

% Convolve an impulse train with a square pulse
%x = mod((1:256),32)==1;
%y = ones(1,16);

% Convolve a square wave with a square pulse
x = [zeros(1,10),mod((1:255),32)<16];
y = ones(1,16)/16;

norm_area = sum(abs(y));

N = max(length(x),length(y))+30;

x = [x, zeros(1, N - length(x))];
y = [y, zeros(1, N - length(y))];

x_conv_y = zeros(1,N);

for i = 1:N
   accumulate = 0;
    y_shift = zeros(1,N);
        
    for j = 1:N
        index = i-j+2;% Fixes some off by one errors
        if index1
            y_shift(j) = y(index);
            accumulate = accumulate + x(j)*y_shift(j);
        end
    end
    x_conv_y(i) = accumulate;
    
    area(x.*y_shift,'FaceColor',[1 1 0])
    hold on
    plot(x)
    %legend('overlapping area','x','y','conv(x,y)')
    plot(y_shift,'g')
    plot(x_conv_y,'r')
    axis([0 N -1 2])
    hold off
    pause(0.001)
end
        

The Short Time Fourier Transform

This demo computes a convolution via a Short Time Fourier Transform. The current signals are set to produce the same result as the previous example. There is a line of code to replace the FFT convolution with a time domain convolution. Note the extra computation time when this is done. If we don't zero pad the time domain signals before taking their FFT, we get time-aliasing (the tail of the segment wraps back around to the beginning). The reason is because we must account for the length of the filtered segment, which is longer than M. Additionally, set L = M and see what happens. This corresponds to ignoring the increased length of the signal when we accumulate the result (failing to overlap-add).

% Convolve a square wave with a square pulse to get a triangle wave
x = [zeros(10,1);mod((1:255)',32)<16];
h = ones(16,1)/16;

N = length(x); % The length of the signal x
M = length(h); % The length of the filter h

N = M*(ceil(N/M)+3); % Make the filter fit in the signal an integer number of times
x = [x; zeros(N - length(x),1)]; % Zero pad up to length N

y = zeros(N,1);

for n = 1 : M : N-M
   segment = x(n:n+M-1);

   % Take the convolution of each section
   % conv_seg = conv(segment,h);
   % FFTs are faster!
   conv_seg = ifft(fft(segment, 2*M-1).*fft(h,2*M-1));
   
   % Add it back into the buffer. Don't forget, the signal 
   % is longer than M now. We must not cut off the tails of 
   % the signals. (Overlap-add method)
   L = length(conv_seg);
   y(n:n+L-1) = y(n:n+L-1) + conv_seg(1:1+L-1);

end
% Look at the signal
plot(y)
        

Finding a signal using a correlation

Using a correlation, we attempt to find when an audio signal is played in the middle of a white noise signal. We do this by looking for a peak in the correlation of the noisy audio signal and the original audio signal. A peak occurs at the point where the two signals match up best. This position will be the correct position of the signal as long as the SNR doesn't get too bad. We plot the probability of recovering the signal as a function of SNR.

x_orig = wavread('dont_speak-no_doubt.wav');
x_orig = x_orig(1265000:1500000,1);
Fs = 22000;
%soundsc(x_orig,Fs) % Listen to song!

x_orig = x_orig/max(abs(x_orig)); % Normalize

% Testing amplitudes. Start at 0dB and decrease to -50dB
test_points = 10.^(-(0:3:50)/20);
data_y = zeros(length(test_points),1);

% Loop variables
j = 1;
i = 0;
% Number of successful recoveries
count = 0;
% Number of data points (for collecting stats)
tests = 100;

hold on

for factor = test_points % Try a new test point
    while (i < tests)
        % Scale down the waveform
        x = x_orig*factor;

        % Add it somewhere in a noisy signal
        L = 500000;
        r = rand(L,1);
        start = round(rand*(length(r)-length(x)));
        r(start: start + length(x) - 1) = r(start: start + length(x) - 1) + x;

        % This is the slow way
        %cross = xcorr(r,x);
        % This is the fast way
        cross = fftshift(ifft(fft(r,2*L-1).*conj(fft(x,2*L-1))));
  
        % We filter (this may not be necessary, it was done to
        % account for the finite length of the correlation, the 
        % edge effects). Don't focus too much on this.
        [b,a] = butter(4, 0.01,'high');
        cross = filtfilt(b,a,cross);
        
        % See if the peak in the autocorrelation is where we expect it
        [val, loc] = max(abs(cross));
        if loc - (length(r) - 1) == start
            % We found the signal!
            count = count + 1;
        end
        
        i=i+1;
    end
    
    i = 0;
    data_y(j) = count/tests;
    if count == 0
        break;
    end
    count = 0;
    j = j + 1
    
    
end

plot(20*log10(test_points), data_y)
xlabel('SNR (dB)')
ylabel('Recovery Probability')
        

Downsampling and Upsampling

This demo goes through some examples of upsampling and downsampling. The imaging and aliasing can be heard if the proper filtering is not done. The final examples show an ideal filter being applied to the signal. Here are some sound files that can be used.

%Choose between this:
x = wavread('TomsDiner_full.wav');
x = x(1:265000,1);
x = x/max(abs(x));
Fs = 44100;
% and this
%t= 0:1/Fs:4;                    
%fo=20;f1=20000;                 
%x=chirp(t,fo,4,f1,'linear');

%Upsampling without filtering. Sounds bad.
y = upsample(x,3);
N = length(y);
omega = linspace(-1,1,N);
plot(omega, abs(fftshift(fft(y/N))))
% Because we upsampled, the sampling rate has changed.
sound(y, Fs)
pause

% Listen. Sounds good. Why? The images are out of 
% the audible range and probably filtered before reaching 
% your speakers
sound(y, Fs*3)
pause

% Downsampling without filtering. Sounds bad. Wrong rate again.
y = downsample(x,2);
N = length(y);
omega = linspace(-1,1,N);
plot(omega, abs(fftshift(fft(y/N))))
sound(y, Fs)
pause

% Sounds bad. Downsampling caused aliasing!
sound(y, Fs/2)
pause
 

% Fix Upsampling by using a filter
K=3 % Upsampling factor
y = upsample(x,3);
N = length(y);
omega = linspace(-1,1,N);
f = (1:N);
f = (fN-N/2/K); % building an ideal filer
Y = fft(y);
max_Y = max(abs(Y));
plot(omega, fftshift(abs((Y/max_Y))))
hold on 
plot(omega, fftshift(f,'r'))
Y_filt = Y.*f;
y_up = real(ifft(Y_filt));
$ Sounds good!
sound(y, Fs*K)
pause

hold off 

% Fix Downsampling by using a filter
L=2 % Downsampling factor
N = length(x);
omega = linspace(-1,1,N);
f = (1:N);
f = (fN-N/2/L); % building an ideal filer
X = fft(x);
X_filt = X.*f;
x_filt = real(ifft(X_filt));
y = downsample(x_filt,L);

max_X = max(abs(X));
plot(omega, fftshift(abs((X/max_X))))
hold on 
plot(omega, fftshift(f),'r')
sound(y, Fs/L)
pause



% Upsample and downsample using only one filter
L=2 
K=3 % Upsampling factor
y = upsample(x,3);
N = length(y);
Y = fft(y);
f = (1:N)';
% Replace two filtering operations
  % f = (fN-N/2/K);
  % Y_filt = Y.*f;
  % y_up = real(ifft(Y_filt));
  % f = (fN-N/2/L);
  % Y_up_filt = Y_up.*f;

% with one
f = (fN-N/2/max(L,K)); % building an ideal filer

% Filter and downsample
Y_up_filt = Y.*f;
y_up_filt = real(ifft(Y_up_filt));
z = downsample(y_up_filt,L);

N = length(z);
omega = linspace(-1,1,N);
plot(omega, abs(fftshift(fft(z/N))))

sound(z, round(Fs*K/L))
        

IIR vs. FIR

Linear filters can be classified according to the length of their impulse response. If this impulse response goes to t=infinity, we call this an Infinite Impulse Response (IIR). If the impulse response truncates before time infinity, we call this a Finite Impulse Response (FIR). This demo shows that IIR filters can easily be constructed using simple recursive filters. These recursive filters have very long impulse responses, but can be calculated very efficiently. They also have much higher rolloff rates than FIR filters. By changing the value of alpha in this demo, it is shown that difference equations can be easily converted from high pass to low pass.

clear 
clc

N = 50
impulse = (1:N==1);


x = impulse;
y = zeros(1,N);
z = zeros(1,N);

alpha = -0.8
% FIR vs IIR
for i = 1:N
    if i>1
        % A simple nonrecursive and recursive filter
        y(i) = x(i) + alpha * x(i-1);
        z(i) = x(i) + alpha * z(i-1);
    else
        y(i) = x(i);
        z(i) = x(i);
    end
    
end
% Look at the impulse responses
subplot(131)
plot(y)
subplot(132)
plot(z,'r')
unstable = 0;
if (max(abs(z)>1))
    unstable = 1;
    title('WARNING: UNSTABLE')
end

% Look at the frequency responses
subplot(133)
N_pad = 1024;
Y = fft(y,N_pad);
Z = fft(z,N_pad);
omega = linspace(-1, 1, N_pad);
plot(omega, abs(fftshift(Y)/max(abs(Y))))
axis([-1 1 0 1]);
hold on
if unstable == 0
    plot(omega, abs(fftshift(Z)/max(abs(Z))), 'r')
else
    plot(omega, abs(fftshift(Z)/max(abs(Z))), 'r--')    
end
hold off
xlabel('Normalized Frequency')
ylabel('Magnitude')


        

Moving Poles & Zeros in the Z Plane

This demo shows creates poles and zeros on the unit circle and gives them some trajectory (changing their distance and/or angle to the origin as a function of time). The frequency response and the phase response are animated showing these changes. The group delay can be displayed instead of the phase response.

circle      = linspace(0,2*pi,1000);

fs = 48000;
snapshots = 150;
zero_radius = linspace(2, 2, snapshots);
zero_theta  = linspace(0, pi,snapshots);

pole_radius = linspace(0.5, 0.5, snapshots);
pole_theta  = linspace(0, pi,snapshots);

for i = 1:snapshots
    P = [];
    Z = [];
    
    % A zero and it's conjugate
    zero = zero_radius(i).*exp(1j*zero_theta(i));
    Z(end+1) = real(zero)+1j*imag(zero);
    Z(end+1) = real(zero)-1j*imag(zero);
    
    % A pole and it's conjugate
    pole = pole_radius(i).*exp(1j*pole_theta(i));
    P(end+1) = real(pole)+1j*imag(pole);
    P(end+1) = real(pole)-1j*imag(pole);
    
    
    % The unit circle
    subplot(131)
    plot(cos(circle), sin(circle), 'k-', 'LineWidth', 1)
    hold on
    plot(real(P),imag(P),'x', 'MarkerSize',12,'LineWidth',2)
    plot(real(Z),imag(Z),'o', 'MarkerSize',12,'LineWidth',2)
    axis([-2 2 -2 2])
    xlabel('real(Z)')
    ylabel('imag(Z)')
    hold off
    
    [B, A]     = zp2tf(Z',P',1);
    [H, omega] = freqz(B,A);
    
    % Plot the magnitude response
    subplot(132)
    plot(fs*omega/2/pi,20*log10(abs(H)))
    axis([0 fs/2 -30 30])
    xlabel('Radians/Sample')
    ylabel('Magnitude')
    
    % Plot the phase response
    subplot(133)
    plot(fs*omega/2/pi,unwrap(atan2(imag(H),real(H))))
    axis([0 fs/2 -10 10])
    xlabel('Radians/Sample')
    ylabel('Phase (Rads)')
    
    % Plot the group delay instead of phase response
    % subplot(133)
    % plot(fs*omega/2/pi, [-diff(unwrap(atan2(imag(H),real(H))));0]/(2*pi/length(omega)))
    % axis([0 fs/2 -10 10])
    % xlabel('Frequency')
    % ylabel('Group Delay  (samples)')

    pause(0.01)
end


        

Digitizing an Analog RC filter

A simple RC circuit can be digitized once we have obtained its continuous time transform function, (1/(1+sRC)), in this case. Using the bilinear transform, we can create a digital model of this filter. The transfer function magnitude responses are shown for the analog and digital case. The digital case shows the frequency response that is nearly identical to the analog frequencies for the passband region. In the higher frequencies, we see that the digital filter rolls off faster. This is due to the frequency warping properties of the bilinear transform.

R = 1000; % Ohms
C = 1e-6; % Farads

fs = 48000;

B = [1];
A = [R*C, 1];

subplot(121)
F = 1:24000
H_analog = freqs(B,A,2*pi*F);
loglog(F,abs(H_analog));
axis([10 20000 .00001 1.1])
title('Transfer function of a simple RC circuit (Analog)');
xlabel('Frequency')
ylabel('Magnitude Response')

% Bilinear likes descending powers of s
[b_z, a_z] = bilinear(B,A,fs);

w = linspace(0, pi, 10000);
H = freqz(b_z,a_z,w);
subplot(122)
f = w*fs/2/pi;
loglog(f,abs(H));
axis([10 20000 .00001 1.1])
title('Transfer function of a simple RC circuit (Digital)');
xlabel('Frequency')
ylabel('Magnitude Response')
        

Bessel, Butterworth, Chebyshev, and Elliptical Filter Zero/Pole Plots

This demo shows four filter types as a pole zero plot. They are shown to be arranged in arcs in the s-plane.

order = 6;

% Show the zeros and poles of the Bessel Filter    
[Z,P,k] = besselap(order);
plot(real(P),imag(P),'x', 'MarkerSize',12,'LineWidth',2)
title('S-Plane')
xlabel('Real')
ylabel('Imag')
axis([-2 2 -2 2])
hold on
plot(real(Z),imag(Z),'o', 'MarkerSize',12,'LineWidth',2)

% Show the zeros and poles of the Chebyshev Filter
[Z,P,k] = cheb1ap(order, 1); % 1dB max of passband error
plot(real(P),imag(P),'rx', 'MarkerSize',12,'LineWidth',2)
plot(real(Z),imag(Z),'ro', 'MarkerSize',12,'LineWidth',2)

% Show the zeros and poles of the Butterworth Filter
[Z,P,k] = buttap(order);
plot(real(P),imag(P),'gx', 'MarkerSize',12,'LineWidth',2)
plot(real(Z),imag(Z),'go', 'MarkerSize',12,'LineWidth',2)

% Show the zeros and poles of the Elliptical Filter
[Z,P,k] = ellipap(order,1,60); % 1dB max of passband error
plot(real(P),imag(P),'kx', 'MarkerSize',12,'LineWidth',2)
plot(real(Z),imag(Z),'ko', 'MarkerSize',12,'LineWidth',2)

grid
        

Bessel, Butterworth, Chebyshev, and Elliptical Filters: Frequency Response and Group Delay

This demo shows the advantages and disadvantages of four popular filter types. The Bessel filter exhibits a high degree of passband phase linearity at the expense of a sharp cutoff and high rolloff. The Butterworth filter is maximally flat in the pass band. The Chebyshev and Elliptical filters have very steep rolloff, but have ripple in their passband. There are a couple of commented out lines that will demonstrate an increase in order of the filters.

order = 4
% for order = 4:12

    w = linspace(.1,100,6000);
    cols = ['b','r','g','k'];
    for i = 1:4
        if i == 1
            % Bessel Filter
            [Z,P,k] = besselap(order);
        elseif i == 2
            % Chebyshev Filter
            [Z,P,k] = cheb1ap(order,1);
        elseif i == 3
            % Butterworth Filter
            [Z,P,k] = buttap(order);
        else
            % Elliptical Filter
            [Z,P,k] = ellipap(order,1,60);
        end

    [b,a] = zp2tf(Z,P,k);
    H = freqs(b,a,w);
    subplot(131)
    loglog(w,abs(H),cols(i));
    title('Magnitude Response')
    xlabel('Frequency')
    ylabel('Magnitude')
    axis([.1, 100, .000001, 10])
    hold on
legend('Bessel','Chebyshev','Butterworth','Elliptical')

    subplot(132)
    plot(w,unwrap(atan2(imag(H),real(H))),cols(i));
    title('Phase Response')
    xlabel('Frequency')
    ylabel('Phase Delay (rads)')
    axis([.1, 2, -8, 0])
    hold on

    subplot(133)
    plot(w(1:end-1),-diff(unwrap(atan2(imag(H),real(H))))/(2*pi/length(w)),cols(i));
    title('Group Delay')
    xlabel('Frequency')
    ylabel('Delay (samples)')
    axis([.1, 2, 0, 260])
    hold on

    end
    legend('Bessel','Chebyshev','Butterworth','Elliptical')

%     pause(2)
%     subplot(131);hold off
%     subplot(132);hold off
%     subplot(133);hold off
% end
        

Bilinear Transform Frequency Warping

This demo shows the transfer characteristic between analog and digital frequencies. The high frequencies in the s-plane get compressed to fit into the finite range of the z-plane

fs = 48000;
T = 1/fs;
f = 0:500000;
plot(f,f,'b-')
hold on

alpha = 2/T
plot(f, alpha/(2*pi)*atan(f/alpha*(2*pi)),'g')
plot(f, fs/2*ones(1,length(f)),'r--')

xlabel('Frequency (from the s plane)')
ylabel('Mapped Frequency')

legend('Original Frequency', 'Bilinear Mapped Frequency', 'Nyquist')
        

Butterworth Filter with Frequency Warping

Depending on whether the filter is properly pre-warped or not, we will find that we do not get the same filter once we have digitized the coefficients that we expected.

fs = 48000;
cutoff = 0.25*fs; % Hz

[B,A] = butter(10, 2*pi*cutoff,'s');
% Look at it in the s plane!
% [H,w] = freqs(B,A)
% semilogx(w,20*log10(abs(H)))
% semilogx(w,20*unwrap(phase(H)))

% No prewarping when making analog filter (no correction)
[b_z, a_z] = bilinear(B,A,fs);
w = linspace(0,pi,2048);
[H,w] = freqz(b_z,a_z,w);
plot(w/pi/2*fs,20*log10(abs(H)),'b')
hold on

% Butter digital corrects for bilinear warping
[b_z, a_z] = butter(10, cutoff/fs*2);
[H,w] = freqz(b_z,a_z,w);
plot(w/pi/2*fs,20*log10(abs(H)),'r--')

xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
axis([10 fs/2 -100 20])

legend('Unwarped','Warped')
        

Resonant Lowpass Frequency Warping

When the filter in question has an audible resonance, it is very important that the peak get mapped from the analog plane to the digital plane properly. This demo shows the deviation in peak placement as the frequency of the peak approaches the high frequencies.

fs = 48000;

for f = linspace(100,10000,1000)
    hold off
    omega = 2*pi*f;
    Q = 10;
    B = [1]
    A = [1/omega^2 1/omega/Q 1]

    [b_z,a_z] = bilinear(B,A,fs);
    w = linspace(0,pi,2000);
    [H,w] = freqz(b_z,a_z,w);
    plot(w/pi/2*fs,20*log10(abs(H)),'b')
    hold on

    [b_z,a_z] = bilinear(B,A,fs,f);
    w = linspace(0,pi,2000);
    [H,w] = freqz(b_z,a_z,w);
    plot(w/pi/2*fs,20*log10(abs(H)),'r')

    xlabel('Frequency (Hz)')
    ylabel('Magnitude (dB)')
    axis([10 fs/2 -100 50])
    pause(0.01)
    
end