%%%% Filtering and Spectral Analysis %%%%
% designing a lowpass filter
f = [0 0.5 0.6 1]; a = [1 1 0 0];
order = 30;
b = firpm(order,f,a);
[h,w] = freqz(b,1,512);
plot(f,a,w/pi,abs(h))
legend('Ideal','firpm Design')
% do the same with order = 5;
% bandpass
f = [0 0.3 0.4 0.6 0.7 1]; a = [0 0 1 1 0 0];
b = firpm(17,f,a);
[h,w] = freqz(b,1,512);
plot(f,a,w/pi,abs(h))
legend('Ideal','firpm Design')
% filter a signal
y = rand(1,8000)*2-1; % white noise
yfilt = filter(b,1,y); % filter it
Hs = spectrum.periodogram; % plot the spectrum
psd(Hs,yfilt,'Fs',8000);
% Exercise: make a harmonic tone with 5 harmonics and suppress the middle harmonic by filtering with a stopband filter
% Better spectrograms for looking at speech signals
%MATLAB has a built-in function specgram() for spectrogram calculation.
%This function divides a long signal into windows and performs a fourier transform on each window, storing complex amplitudes in a table in which the columns represent time and the rows represent frequency.
%Call the specgram function on a signal x with a call:
[B,f,t] = specgram(y,nfft,fs,window,overlap)
%Here, nfft is size of the fourier transform, use this to control the number of frequency samples on the vertical axis of the spectrogram; fs is the sampling rate of the input signal; window is the number of input samples per vertical slice, use this to control the analysis bandwidth; overlap is the number of samples by which adjacent windows overlap, use this to control the number of vertical slices per second. The output B is a table of complex amplitudes; f is a vector of frequency values used to label the rows, and t is a vector of time values used to label the columns.
load mtlb
% calculate the table of amplitudes
[B,f,t]=specgram(mtlb,1024,Fs,256,192);
%
% calculate amplitude in dB
powerdB=20*log10(abs(B)+eps);
%
% plot top 10dB as image
imagesc(t,f,(max(powerdB,-10)*-1));
%
% label plot
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
colormap(gray);
% Windowing a signal
% Often we want to analyse a long signal in overlapping short sections called “windows”. For example we may want to calculate an average spectrum, or to calculate a spectrogram. Unfortunately we cannot simply chop the signal into short pieces because this will cause sharp discontinuities at the edges of each section. Instead it is preferable to have smooth joins between sections. Raised cosine windows are a popular shape for the joins:
% You can use the MATLAB function hamming() to design smooth windows of a given length, and then you can use code such as the following to divide the signal into sections:
% divide up a signal into windows
nx = length(x); % size of signal
w = hamming(128); % hamming window
nw = length(w); % size of window
pos=1;
while (pos+nw <= nx) % while enough signal left
y = x(pos:pos+nw-1).*w; % make window y
%%%% process window y %%%%
pos = pos + nw/2; % next window
end
% divide up a signal into windows
x = mtlb;
nx = length(x); % size of signal
w = hamming(128); % hamming window
nw = length(w); % size of window
nsteps = floor(nx/(nw/2))-2;
spec = zeros(512,nsteps);
f = linspace(0,Fs/2,512);
pos=1;
for I = 1:nsteps
y = x(pos:pos+nw-1).*w; % make window y
%%%% process window y %%%%
Y = fft(y,1024);
magnitude = abs(Y);
powerdB = 20*log10(magnitude+eps);
spec(:,I) = powerdB(1:512);
pos = pos + nw/2; % next window
end
% plot top 10dB as image
imagesc(t,f,(max(spec,-10)*-1));
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
colormap(gray);
% Now do something more fun - estimate pitch or formants in each window
% record 1s audio (say /a/)
recObj = audiorecorder;
recordblocking(recObj,1);
vowel = getaudiodata(recObj);
FS = 8000;
t = (0:length(a)-1)/FS;
plot(t,a)
sound(a,8000)
specgram(a);
% spectrum
nyq = FS/2;
Y = fft(a);
magnitude = abs(Y(1:nyq));
powerdB = 20*log10(magnitude+eps);
f = 0:nyq-1;
semilogx(f,powerdB)
plot(f,powerdB)
[jwd,idx] = max(powerdB)
f(idx)
f(b)*[1:15]
% the cepstrum
C=fft(log(magnitude));
figure; plot(abs(C(8:2000)))
% estimating pitch by autocorrelation
ms20=FS/50; % minimum speech pitch at 50Hz = 20ms period (could also do 0.02/ST)
ms2 = FS/500; maximum at 500Hz = 2ms peris
r = xcorr(x,ms20,'coeff');
d=(-ms20:ms20)/FS;
plot(d*1000,r)
r_speech = r((ms20+ms2):end);
d_speech = d((ms20+ms2):end);
[jwd,idx] = max(r_speech);
1/d_speech(idx)
% Homework:
% Put the above code to estimate pitch using autocorrelation into the code that devides a signal into windows so that you get a pitch estimate for each window. Record a few seconds of sound with changing pitch (try singing) and plot a curve (pitch frequency vs. time) showing how pitch changes over time in that sound.