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

% 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

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

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

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')

% 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))));

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');

% 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

% 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)

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')

%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 = (f N-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 = (f N-N/2/K); % Y_filt = Y.*f; % y_up = real(ifft(Y_filt)); % f = (f N-N/2/L); % Y_up_filt = Y_up.*f; % with one f = (f N-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))

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')

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

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')

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

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

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')

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')

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