Шумный входной сигнал

Вход является звуковым сигналом, выбранным на уровне 22 050 Гц. Объект dsp.AudioFileReader читает этот сигнал в кадрах 512 выборок. Звуковой сигнал повреждается белым Гауссовым шумом, который имеет стандартное отклонение 0,05. Используйте объект audioDeviceWriter проигрывать шумный звуковой сигнал к аудио устройству вашего компьютера.

FrameLength = 512;
afr = dsp.AudioFileReader('speech_dft.wav',...
    'SamplesPerFrame',FrameLength);
adw = audioDeviceWriter('SampleRate',afr.SampleRate);

noiseStd = 0.05;
while ~isDone(afr)
    cleanAudio = afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    adw(noisyAudio);
end
reset(afr)

Инициализируйте кратковременные и обратные кратковременные объекты преобразования Фурье

Инициализируйте объекты dsp.ISTFT и dsp.STFT. Установите длину окна, равную входной длине кадра и длине транзитного участка к 16. Длина перекрытия является различием между длиной окна и длиной транзитного участка, OL = WL – HL. Установите длину БПФ на 1 024.

WindowLength = FrameLength;
HopLength = 16;
numHopsPerFrame = FrameLength / 16;
FFTLength = 1024;

Окно, используемое, чтобы вычислить STFT и ISTFT, является периодическим окном hamming с длиной 512. Флаг ConjugateSymmetricInput объекта istf установлен в true, указав, что вывод объекта istf является сопряжено-симметричным сигналом.

win = hamming(WindowLength,'periodic');
stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);

Получите средство оценки

Следующий шаг должен задать параматерей средства оценки усиления. Это усиление применяется к шумному спектру, чтобы ослабить поддиапазоны с более высокой шумовой степенью и подъемником поддиапазоны с меньшей шумовой степенью.

dec = 16;
alpha = 15;
stftNorm = (sum(win.*win) / dec).^2;

Спектральное затухание

Накормите звуковым сигналом stf одной длиной транзитного участка за один раз. Примените предполагаемое усиление к преобразованному сигналу. Восстановите denoised версию исходного речевого сигнала путем выполнения обратного преобразования Фурье на отдельных диапазонах частот. Проигрывайте denoised звуковой сигнал к аудио устройству компьютера.

while ~isDone(afr)
    cleanAudio =  afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    y = zeros(FrameLength,1); % y holds the denoised audio frame
    
    % Feed audio to stft one hop-length at a time
    for index = 1:numHopsPerFrame        
        X = stf(noisyAudio((index-1)*HopLength+1:index*HopLength));        
        % Gain estimator
        Z = abs(X).^2 / (noiseStd^2 * alpha) / stftNorm;
        Z(Z<=1) = 1;
        Z = 1 - 1./Z;
        Z = sign(Z) .* sqrt(abs(Z));
        X = X .* Z;        
        % Convert back to time-domain
        y((index-1)*HopLength+1:index*HopLength) = istf(X);        
    end    
    % Listen to denoised audio:
    adw(y);
end

Функция iscola

Функция iscola проверяет, что заданное окно и перекрытие удовлетворяют ограничение COLA, чтобы гарантировать, что обратное кратковременное преобразование Фурье (ISTFT) приводит к совершенной реконструкции для неизмененных спектров. Функция возвращает логический true, если комбинация входных параметров совместима КОЛОЙ и логический false если нет. Аргумент method функции установлен в 'ola' или 'wola' в зависимости от того, добавляет ли использование метода инверсии взвешенное перекрытие - (ВОЛЯ).

Проверяйте, является ли окном hann() длины 120 выборок и продолжительность перекрытия 60 выборок совместимый COLA.

winLen = 120;
overlapLen = 60;
win = hann(winLen,'periodic');
tf = iscola(win,overlapLen,'ola')

tf = logical
   1

Инициализация

Инициализируйте dsp.STFT и Системные объекты dsp.ISTFT с этим окном hann, которое является совместимым COLA. Установите длину БПФ равняться длине окна.

frameLen = winLen-overlapLen;
stf = dsp.STFT('Window',win,'OverlapLength',overlapLen,'FFTLength',winLen);
istf = dsp.ISTFT('Window',win,'OverlapLength',overlapLen,'WeightedOverlapAdd',0);

Восстановите данные

Вычислите STFT случайного сигнала. Установите длину входного сигнала равняться длине транзитного участка (длина окна – длина перекрытия). Поскольку окном является совместимый COLA, ISTFT, неизмененные спектры отлично восстанавливают исходный сигнал временного интервала.

Чтобы подтвердить, сравните вход, x к восстановленному выводу, y. Из-за задержки, введенной объектами, восстановленный вывод переключен вовремя по сравнению с входом. Поэтому, чтобы выдержать сравнение, возьмите норму различия между восстановленным выводом, y и предыдущим входом, xprev. Норма является очень маленькой, указывая, что выходной сигнал является отлично восстановленной версией входного сигнала.

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 2.6883e-15

ISTFT со взвешенным перекрытием - добавляют (ВОЛЯ)

В ВОЛИ второе окно вызвало окно синтеза, $\mathit{f}\left(\mathit{n}\right)$, применяется после операции IFFT и перед перекрытием - добавляют. Окна синтеза и анализа обычно идентичны и обычно получаются путем пущения квадратного корня из окон, удовлетворяющих COLA (таким образом, гарантирующий совершенную реконструкцию).

Функция iscola

Проверяйте, является ли окном sqrt(hann()) длины 120 выборок и продолжительность перекрытия 60 выборок совместимая ВОЛЯ. Установите аргумент method функции iscola к 'wola'. Вывод функции iscola является 1 указанием, что этим окном является совместимая ВОЛЯ.

winWOLA = sqrt(hann(winLen,'periodic'));
tfWOLA = iscola(winWOLA,overlapLen,'wola')

tfWOLA = logical
   1

Восстановите данные с ВОЛЕЙ

Выпустите dsp.STFT и Системные объекты dsp.ISTFT и установите окно на окно sqrt(hann(winLen,'periodic')). Чтобы использовать взвешенное перекрытие - прибавляют сторону ISTFT, устанавливают 'WeightedOverlapAdd' на true.

release(stf);
release(istf);
stf.Window = winWOLA;
istf.Window = winWOLA;
istf.WeightedOverlapAdd = true;

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 2.4499e-15

norm различия между входным сигналом и восстановленным сигналом является очень маленьким указанием, что сигнал был восстановлен отлично.