diff --git a/func/fft_.m b/func/fft_.m new file mode 100644 index 0000000..ed209aa --- /dev/null +++ b/func/fft_.m @@ -0,0 +1,14 @@ +function [frequencies, values] = fft_(signal, sample_frequency) + +L=length(signal); + +Y = fft(signal); +P2 = abs(Y); % two-sided spectrum +% P2 = abs(Y/L); % two-sided spectrum +P1 = P2(1:L/2+1); % single-sided spectrum +P1(2:end-1) = 2*P1(2:end-1); +frequencies = sample_frequency*(0:(L/2))/L; +values = P1; + +end + diff --git a/func/spectro.m b/func/spectro.m new file mode 100644 index 0000000..8741cb7 --- /dev/null +++ b/func/spectro.m @@ -0,0 +1,14 @@ +function spectro(signal, sample_frequency, windows, overlap_interval) + +sample_overlap = (overlap_interval / 1000) * sample_frequency; + +sample_size = size(signal); +%window_size = round(sample_size(1) / ((windows + 1)/2)) + +% Turn windows into window width in samples, take into account overlap +window_size = round((sample_size(1) + (windows + 1) * sample_overlap) / (windows+1)); + +spectrogram(signal, window_size, sample_overlap, [], sample_frequency, 'yaxis'); + +end + diff --git a/lpss.m b/lpss.m new file mode 100644 index 0000000..9e4d8ea --- /dev/null +++ b/lpss.m @@ -0,0 +1,83 @@ +%% lpss.m +%% +%% Coursework script + +close all;clear all;clc; + +LPC_ORDER = 8; +DISPLAY_SAMPLES = 1000; +WINDOW_NUMBER = 10; +WINDOW_OVERLAP = 5; % ms + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% READ SIGNAL +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +[y, Fs] = audioread('samples/hood_m.wav'); +L = length(y) % number of samples +DISPLAY_SAMPLES = min([DISPLAY_SAMPLES L]); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% LPC +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +a = aryule(y, LPC_ORDER) % signal, filter order +est_y = filter(0.02, a, y); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% PREDICTION ERROR +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +e = y - est_y; +[acs, lags] = xcorr(e,'coeff'); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% COMPARE TWO SIGNALS +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% x = 1:DISPLAY_SAMPLES; +% figure(1) +% plot(x, y(end-DISPLAY_SAMPLES+1:end), x, est_y(end-DISPLAY_SAMPLES+1:end), '--') + +%plot(x, y(end-DISPLAY_SAMPLES+1:end)) +%plot(x, est_y(end-DISPLAY_SAMPLES+1:end)) + +% grid +% xlabel('Sample Number') +% ylabel('Amplitude') +% legend('Original signal','LPC estimate') + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% PLOT AUTOCORRELATION +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% figure(2) +% plot(lags, acs) +% grid +% xlabel('Lags') +% ylabel('Normalized Autocorrelation') +%ylim([-0.1 1.1]) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% PLOT FREQUENCY RESPONSE +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% CALCULATE FFT +[freq_dom_freqs, freq_dom_vals] = fft_(y, Fs); + +% GET FILTER RESPONSE +[h, filter_freqs] = freqz(1, a, length(freq_dom_freqs), Fs); + +figure(3) +plot(freq_dom_freqs, 20*log10(freq_dom_vals), 'r', filter_freqs, 20*log10(abs(h)), 'b') +% plot(w/pi, 20*log10(abs(h))) +grid +xlabel('Frequency (Hz)') +ylabel('Magnitude (dB)') +legend('Original Signal', 'LPC Filter') + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% PLOT ORIGINAL SPECTROGRAM +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% figure(4) +% spectro(y, Fs, WINDOW_NUMBER, WINDOW_OVERLAP); +% colormap bone + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% PLAY +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%sound(y, Fs); diff --git a/lpss_autocorr.m b/lpss_autocorr.m index 0a5b564..6d8d776 100644 --- a/lpss_autocorr.m +++ b/lpss_autocorr.m @@ -3,12 +3,10 @@ %% Load wav and run through autocorr from econometrics toolbox %% Used max range of m lags calculated from sample frequency -close all; -clear all; +close all;clear all;clc; [y, Fs] = audioread('samples/hood_m.wav'); -sample_size = size(y); -L = sample_size(1); % number of samples +L = length(y); % number of samples f0 = 60; % low-pitched male speech %f0 = 600; % children diff --git a/lpss_cepstrum.m b/lpss_cepstrum.m new file mode 100644 index 0000000..18e4231 --- /dev/null +++ b/lpss_cepstrum.m @@ -0,0 +1,14 @@ +%% lpss_cepstrum.m +%% +%% Load wav and calculate cepstrum + +close all;clear all;clc; + +% READ SIGNAL +[y, Fs] = audioread('samples/hood_m.wav'); + +t = (0:length(y) - 1); +c = rceps(y); +% c = cceps(y); +plot(t, c) +xlabel('Quefrency') \ No newline at end of file diff --git a/lpss_fft.m b/lpss_fft.m index 385e38e..99a40f5 100644 --- a/lpss_fft.m +++ b/lpss_fft.m @@ -3,12 +3,10 @@ %% Load wav and plot in frequency domain using fft %% Construct proper one-sided function for display -close all; -clear all; +close all;clear all;clc; [y, Fs] = audioread('samples/hood_m.wav'); -sample_size = size(y); -L = sample_size(1); % number of samples +L = length(y); % number of samples % PLOT TIME DOMAIN figure(1); diff --git a/lpss_lpc.m b/lpss_lpc.m index 37c6346..4442ee4 100644 --- a/lpss_lpc.m +++ b/lpss_lpc.m @@ -2,45 +2,48 @@ %% %% Load wav and calculate LPC coeffs -close all; -clear all; +close all;clear all;clc; -ORDER = 40; +ORDER = 50; DISPLAY_SAMPLES = 1000; % READ SIGNAL [y, Fs] = audioread('samples/hood_m.wav'); -sample_size = size(y); -L = sample_size(1); % number of samples +L = length(y) % number of samples DISPLAY_SAMPLES = min([DISPLAY_SAMPLES L]); -for ITER=1:ORDER +% CALCULATE FFT +[freq_dom_freqs, freq_dom_vals] = fft_(y, Fs); + +for ITER=1:5:ORDER % LPC a = lpc(y,ITER); % signal, filter order - est_y = filter(0.02, a, y); - - % PREDICTION ERROR - e = y - est_y; - [acs, lags] = xcorr(e,'coeff'); - - % COMPARE TWO SIGNALS - x = 1:DISPLAY_SAMPLES; - figure(2) - plot(x, y(end-DISPLAY_SAMPLES+1:end), x, est_y(end-DISPLAY_SAMPLES+1:end), '--') - %plot(x, y(end-DISPLAY_SAMPLES+1:end)) - %plot(x, est_y(end-DISPLAY_SAMPLES+1:end)) + + % COMPARE FREQ DOMAIN + [h, filter_freqs] = freqz(1, a, length(freq_dom_freqs), Fs); + + figure(1) + plot(freq_dom_freqs, 20*log10(freq_dom_vals), 'r', filter_freqs, 20*log10(abs(h)), 'b') + % plot(w/pi, 20*log10(abs(h))) grid - xlabel('Sample Number') - ylabel('Amplitude') - legend('Original signal','LPC estimate') + xlabel('Frequency (Hz)') + ylabel('Magnitude (dB)') + legend('Original Signal', 'LPC Filter') - % PLOT AUTOCORRELATION + % COMPARE TWO SIGNALS TIME DOMAIN +% est_y = filter(0.02, a, y); +% x = 1:DISPLAY_SAMPLES; % figure(2) -% plot(lags, acs) +% plot(x, y(end-DISPLAY_SAMPLES+1:end), x, est_y(end-DISPLAY_SAMPLES+1:end), '--') +% %plot(x, y(end-DISPLAY_SAMPLES+1:end)) +% %plot(x, est_y(end-DISPLAY_SAMPLES+1:end)) % grid -% xlabel('Lags') -% ylabel('Normalized Autocorrelation') -% ylim([-0.1 1.1]) +% xlabel('Sample Number') +% ylabel('Amplitude') +% legend('Original signal','LPC estimate') + + + pause(0.5) end \ No newline at end of file diff --git a/lpss_play.m b/lpss_play.m index 6eeb431..263cc9c 100644 --- a/lpss_play.m +++ b/lpss_play.m @@ -2,8 +2,7 @@ %% %% Load wav and play -close all; -clear all; +close all;clear all;clc; [y, Fs] = audioread('samples/ee.wav'); size(y) diff --git a/lpss_spectrogram.m b/lpss_spectrogram.m index ef0171f..dbaa59b 100644 --- a/lpss_spectrogram.m +++ b/lpss_spectrogram.m @@ -2,15 +2,13 @@ %% %% Load wav and plot spectrogram -close all; -clear all; +close all;clear all;clc; WINDOW_NUMBER = 20; % matrix by audio channel and sample frequency [y, Fs] = audioread('samples/hood_m.wav'); -sample_size = size(y); -window_size = round(sample_size(1) / ((WINDOW_NUMBER + 1)/2)); +window_size = round(length(y) / ((WINDOW_NUMBER + 1)/2)); freqRes = Fs / window_size timeRes = window_size / Fs diff --git a/lpss_spectrogram_iter.m b/lpss_spectrogram_iter.m index 219fa1f..542194e 100644 --- a/lpss_spectrogram_iter.m +++ b/lpss_spectrogram_iter.m @@ -2,18 +2,16 @@ %% %% Load wav and plot spectrogram for iterating window sizes -close all; -clear all; +close all;clear all;clc; %WINDOW_NUMBER = 20; % matrix by audio channel and sample frequency [y, Fs] = audioread('samples/hood_m.wav'); -sample_size = size(y); for WINDOW_NUMBER=1:50 - window_size = round(sample_size(1) / ((WINDOW_NUMBER + 1)/2)); + window_size = round(length(y) / ((WINDOW_NUMBER + 1)/2)); freqRes = Fs / window_size; timeRes = window_size / Fs;