Audio & DSP

fs = 2000; tt = 0:1/fs:1;
a = sin(2*pi*40*tt);
b = sin(2*pi*44*tt);
plot(tt, a + b, 'LineWidth', 0.8), hold on
plot(tt, 2*abs(cos(2*pi*2*tt)), 'r--', 'LineWidth', 1.5)
plot(tt, -2*abs(cos(2*pi*2*tt)), 'r--', 'LineWidth', 1.5)
hold off
title('40 Hz + 44 Hz → a 4 Hz beat')
xlabel('time (s)'), legend('sum', 'envelope')

% hear it: 440 Hz + 444 Hz wobbles 4 times a second
fs2 = 22050; ts = 0:1/fs2:3;
soundsc(sin(2*pi*440*ts) + sin(2*pi*444*ts), fs2)

% Short-time spectra of a sweep from 100 Hz to 1700 Hz, stacked over time
rng(1)
fs = 4000; tt = 0:1/fs:2;
x = sin(2*pi*(100*tt + 400*tt.^2)) + 0.3*randn(size(tt));
win = 256; hop = 64;
w = 0.5 - 0.5*cos(2*pi*(0:win-1)/(win-1));      % Hann window
S = []; centers = [];
for k = 1:hop:numel(x) - win
  [f, A] = spectrum(x(k:k+win-1) .* w, fs);
  S = [S, 20*log10(A' + 1e-6)];
  centers = [centers, tt(k + win/2)];
end
imagesc(centers, f, S)
axis xy
colormap magma
clim([-60 0])
colorbar
xlabel('time (s)'), ylabel('frequency (Hz)')
title('Spectrogram (dB)')
soundsc(x, fs)                                  % listen to the sweep

rng(9)
fs = 1000; tt = 0:1/fs:0.5;
clean = sin(2*pi*8*tt);
noisy = clean + 0.5*randn(size(tt));
N = 25; h = ones(1, N) / N;               % 25 equal taps
y = conv(noisy, h); y = y(1:numel(tt));

subplot(2, 1, 1)
plot(tt, noisy, 'Color', [0.72 0.72 0.72]), hold on
plot(tt, y, 'LineWidth', 2), plot(tt, clean, 'k--'), hold off
legend('noisy', 'filtered', 'clean')
title('25-tap moving average (note the delay of 12 samples)')

subplot(2, 1, 2)
f = linspace(0.5, 500, 800);
H = abs(sin(pi*f*N/fs) ./ (N*sin(pi*f/fs)));
plot(f, 20*log10(H + 1e-9), 'LineWidth', 1.5)
ylim([-60 5])
xlabel('frequency (Hz)'), ylabel('gain (dB)'), title('frequency response')

% listen: one second of white noise, then the same noise through the filter (duller)
fs2 = 22050; n = 0.25*randn(1, fs2);
sound([n, 4*filter(ones(1, N)/N, 1, n)], fs2)

fs = 8000; tt = 0:1/fs:1.2;
A = 0.02; D = 0.15; S = 0.5; R = 0.4;     % attack, decay (s), sustain level, release (s)
env = zeros(size(tt));
on = tt < A;                env(on) = tt(on) / A;
on = tt >= A & tt < A + D;  env(on) = 1 - (1 - S) * (tt(on) - A) / D;
on = tt >= A + D & tt < 0.8; env(on) = S;
on = tt >= 0.8;             env(on) = S * exp(-(tt(on) - 0.8) / (R/4));

tone = zeros(size(tt));
for k = 1:6, tone = tone + sin(2*pi*220*k*tt) / k^1.5; end   % harmonics of 220 Hz
x = env .* tone;

subplot(2, 1, 1)
plot(tt, x, 'LineWidth', 0.5), hold on
plot(tt, env * max(abs(tone)), 'r', 'LineWidth', 2), hold off
title('220 Hz tone shaped by attack–decay–sustain–release'), xlabel('time (s)')
subplot(2, 1, 2)
[f, P] = spectrum(x, fs);
plot(f, P, 'LineWidth', 1.2), xlim([0 1600])
title('harmonic spectrum'), xlabel('frequency (Hz)')
sound(0.8 * x / max(abs(x)), fs)

% A burst of noise circulating in a delay line, averaged each pass, sounds like a plucked string
rng(3)
fs = 22050;
notes = [196 247 294 392];                % G3 B3 D4 G4: a G major chord, strummed
y = zeros(1, round(2.5*fs));
for k = 1:numel(notes)
  N = round(fs / notes(k));               % delay length sets the pitch
  burst = [rand(1, N) - 0.5, zeros(1, 2*fs - N)];
  % y(n) = x(n) + 0.996 * (y(n-N) + y(n-N-1)) / 2 : a delay line that averages and decays
  out = filter(1, [1, zeros(1, N - 1), -0.498, -0.498], burst);
  start = round((k - 1) * 0.06 * fs);
  y(start + 1:start + numel(out)) = y(start + 1:start + numel(out)) + out;
end
soundsc(y, fs)

tt = (0:numel(y) - 1) / fs;
subplot(2, 1, 1)
plot(tt, y, 'LineWidth', 0.4), xlabel('time (s)'), title('Strummed G major chord')
subplot(2, 1, 2)
[f, P] = spectrum(y(1:fs/2), fs);
plot(f, 20*log10(P + 1e-6)), xlim([0 2000])
xlabel('frequency (Hz)'), ylabel('dB'), title('spectrum of the first half second')

% Every key on a phone pad is the sum of one low and one high tone
keys = '1234567890';
low  = [697 697 697 770 770 770 852 852 852 941];
high = [1209 1336 1477 1209 1336 1477 1209 1336 1477 1336];
fs = 8000; tone = 0:1/fs:0.18; gap = zeros(1, round(0.07*fs));
number = '0612345678';
y = [];
for ch = number
  k = find(keys == ch);
  y = [y, 0.5*sin(2*pi*low(k)*tone) + 0.5*sin(2*pi*high(k)*tone), gap];
end
sound(0.9 * y, fs)

% decode the first digit again from its spectrum
seg = y(1:numel(tone));
[f, P] = spectrum(seg .* hann(numel(seg))', fs);
plot(f, P, 'LineWidth', 1.2), xlim([500 1700])
xline(low(10), ':', '941 Hz'), xline(high(10), ':', '1336 Hz')
title(['Dialling ' number ' — the first key is "0"'])
xlabel('frequency (Hz)'), ylabel('amplitude')

% Two columns = left and right channel. Equal-power panning keeps the loudness steady.
fs = 22050; tt = 0:1/fs:4;
src = sin(2*pi*330*tt) .* (0.6 + 0.4*sin(2*pi*3*tt));   % a pulsing tone
pan = (1 - cos(pi*tt/4)) / 2;                         % 0 = left … 1 = right
L = cos(pan*pi/2) .* src;
R = sin(pan*pi/2) .* src;
sound(0.7 * [L', R'], fs)                            % N×2 matrix: stereo
plot(tt, cos(pan*pi/2), tt, sin(pan*pi/2), 'LineWidth', 2)
legend('left gain', 'right gain'), xlabel('time (s)'), title('Equal-power pan')