//Spring 2010.
//Project 06.
//Copyright @ 2010 antonio081014 ;
//All codes are in Matlab.
Report Link:
Task:
1. Speech enhancement by using spectral subtraction;
%% Proj06;
clear all;
clc;
close all;
fd = fopen('chj.16', 'rb', 'b');
[spch, K] = fread(fd, 'int16');
fclose(fd);
% Remove the explosive in the background;
spch(28080:28100) = 0;
sspch = (spch - mean(spch));
y = sspch /max(sspch);
fs = 8000;
FL = 20 * fs / 1000; % 20ms;
NF = 1024;
frame = zeros(4, FL);
frame(1,:) = y(2241:2400);
frame(2,:) = y(6001:6160);
frame(3,:) = y(14912:15071);
frame(4,:) = y(33226:33385);
STFT = zeros(4, NF); % short-time spectrum;
STauto = zeros(4, 2 * FL -1); % short-time autocorrelation;
STPow = zeros(4, NF); % short-time
for i=1:4
STFT(i, :) = fft(frame(i,:).*hamming(FL)', NF);
STauto(i, :) = xcorr(frame(i,:));
STPow(i, :) = abs(fft(frame(i,:).*hamming(FL)', NF)).^2; %Using hamming window to reduce the error.
% STPow(i, :) = abs(fft(frame(i,:), NF)).^2;
end
AvergSTPow = sum(STPow, 1) / 4;
Newspch = zeros(1, length(y));
Fstart = 1;
Fend = Fstart + FL - 1;
overlap = FL / 2;
Fcnt = 0;
while(Fend <= length(y))
data = y(Fstart:Fend).*hamming(FL);
FFTdata = fft(data, NF);
Re = max(abs((FFTdata)').^2 - AvergSTPow, 0); %Get the Power of the spectral differences;
phase = angle((FFTdata)'); % Keep the phase of the original speech;
Newdata = sqrt(Re) .* (cos(phase) + i * sin(phase));
% Newdata = sqrt(Re) .* exp(i * phase); % This way doesn't work. I don't know why.
tmp = real(ifft(Newdata));
Newspch(1, Fstart: Fend) = Newspch(1, Fstart: Fend) + tmp(1, 1:160);
% don't use specified length when using ifft.
% Newspch(1, Fstart: Fend) = Newspch(1, Fstart: Fend) + real(ifft(Newdata, FL)); % Reconstruct the speech with 50% overlap.
Fstart = Fstart + overlap;
Fend = Fstart + FL - 1;
Fcnt = Fcnt + 1;
end
% soundsc(y, fs);
% soundsc(Newspch, fs);
figure; hold on;
plot(y, 'r'); plot(Newspch, 'g');