Сгенерируйте трехканальный сигнал, состоящий из трех различных щебетов, дискретизированных с частотой 1 кГц в течение 1 секунды.

Вычислите STFT многоканального сигнала с помощью периодического окна Хэмминга длины 256 и длины перекрытия 15 выборки.

fs = 1e3;
t = 0:1/fs:1-1/fs;
x = [chirp(t,100,1,300,'quadratic',45,'concave');
      chirp(t,200,1,600,'quadratic',[],'convex');
      chirp(t,300,1,500,'logarithmic')]'; 
  
[S,F,T] = stft(x,fs,'Window',hamming(256,'periodic'),'OverlapLength',15);

Постройте график оригинальной и реконструированной версий первого и второго каналов.

[ix,ti] = istft(S,fs,'Window',hamming(256,'periodic'),'OverlapLength',15);

plot(t,x(:,1)','LineWidth',1.5)
hold on
plot(ti,ix(:,1)','r--')
hold off
legend('Original Channel 1','Reconstructed Channel 1')

plot(t,x(:,2)','LineWidth',1.5)
hold on
plot(ti,ix(:,2)','r--')

legend('Original Channel 2','Reconstructed Channel 2')

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

Фазовый вокодер принимает STFT сигнала с окном анализа размера скачка ${\mathit{R}}_1$ а затем выполняет ISTFT с окном синтеза размера скачка ${\mathit{R}}_2$. Таким образом, вокодер использует в своих интересах метод WOLA. Чтобы растянуть сигнал во времени, в окне анализа используется большее количество выборок перекрытия, чем синтез. В результате на выходе больше выборок, чем на входе (${\mathit{N}}_{\mathit{S},\mathrm{Out}}>{\mathit{N}}_{\mathit{S},\mathrm{In}}$), хотя содержимое остается неизменным. Теперь можно поднять шкалу этого сигнала путем воспроизведения его с более высокой частотой дискретизации, которая выдает сигнал с исходной длительностью, но с более высоким тангажом.

Загрузите аудио файла, содержащее фрагмент Hallelujah Chorus Генделя, выбранный с частотой дискретизации 8192 Гц.

load handel

Спроектируйте окно корня-Ханна длиной 512. Установите длину перекрытия анализа равную 192, и длину перекрытия синтеза равную 166.

wlen = 512;
win = sqrt(hann(wlen,'periodic'));
noverlapA = 192;
noverlapS = 166;

Реализуйте фазовый вокодер с помощью окна анализа перекрытия 192 и окна синтеза перекрытия 166.

S = stft(y,Fs,'Window',win,'OverlapLength',noverlapA);
iy = istft(S,Fs,'Window',win,'OverlapLength',noverlapS);


%To hear, type soundsc(w,Fs), pause(10), soundsc(iw,Fs);

Если окна анализа и синтеза совпадают, но длина перекрытия изменена, будет дополнительная прибыль/потеря, которую вам нужно будет настроить. Это общий подход к реализации фазы вокодера.

Вычислите коэффициент скачка и используйте его, чтобы настроить усиление восстановленного сигнала. Также вычислите частоту сдвинутых по тангажу данных с помощью коэффициента скачка.

hopRatio = (wlen-noverlapS)/(wlen-noverlapA);
iyg = iy*hopRatio;
Fp = Fs*hopRatio;


%To hear, type soundsc(iwg,Fs), pause(15), soundsc(iwg,Fp);

Постройте график исходного сигнала и растянутого по времени сигнала с фиксированным усилением.

plot((0:length(iyg)-1)/Fs,iyg,(0:length(y)-1)/Fs,y)
xlabel('Time (s)')
xlim([0 (length(iyg)-1)/Fs])
legend('Time Stretched Signal with Fixed Gain','Original Signal','Location','best')

Сравните растянутый по времени сигнал и сдвинутый по тангажу сигнал на том же графике.

plot((0:length(iy)-1)/Fs,iy,(0:length(iy)-1)/Fp,iy)
xlabel('Time (s)')
xlim([0 (length(iyg)-1)/Fs])
legend('Time Stretched Signal','Pitch Shifted Signal','Location','best') 

Чтобы лучше понять эффект данных перемены тангажа, рассмотрите следующую синусоиду частоты Fs более 2 секунд.

t =  0:1/Fs:2;
x = sin(2*pi*10*t);

Вычислите короткое преобразование Фурье и обратное короткое преобразование Фурье с длинами перекрытия 192 и 166 соответственно.

Sx = stft(x,Fs,'Window',win,'OverlapLength',noverlapA);
ix = istft(Sx,Fs,'Window',win,'OverlapLength',noverlapS);

Постройте график исходного сигнала на одном графике и растянутого по времени и смещенного по тангажу сигнала на другом.

subplot(2,1,1)
plot((0:length(ix)-1)/Fs,ix,'LineWidth',2) 
xlabel('Time (s)')
ylabel('Signal Amplitude')
xlim([0 (length(ix)-1)/Fs])
legend('Time Stretched Signal') 



subplot(2,1,2)
hold on
plot((0:length(x)-1)/Fs,x)
plot((0:length(ix)-1)/Fp,ix,'--','LineWidth',2)
legend('Original Signal','Pitch Shifted Signal','Location','best')
hold off
xlabel('Time (s)')
ylabel('Signal Amplitude')
xlim([0 (length(ix)-1)/Fs])

Сгенерируйте комплексную синусоиду частоты 1 кГц и длительности 2 секунды.

fs = 1e3;
ts = 0:1/fs:2-1/fs;

x = exp(2j*pi*100*cos(2*pi*2*ts));

Проектируйте периодическое окно Ханна длиной 100 и установите количество выборок перекрытия равным 75. Проверьте длину окна и перекрытия для податливости COLA.

nwin = 100;
win = hann(nwin,'periodic');
noverlap = 75;

tf = iscola(win,noverlap)

tf = logical
   1

Обнулите сигнал, чтобы удалить эффекты края. Чтобы избежать усечения, дополните входной сигнал нулями так, чтобы

- целое число. Установите длину БПФ равной 128. Вычислите кратковременное преобразование Фурье комплексного сигнала.

nPad = 100;
xZero = [zeros(1,nPad) x zeros(1,nPad)];
fftlen = 128;
s = stft(xZero,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);

Вычислите обратное кратковременное преобразование Фурье и удалите нули для идеальной реконструкции.

[is,ti] = istft(s,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);
is(1:nPad) = [];
is(end-nPad+1:end) = [];
ti = ti(1:end-2*nPad);

Постройте график реальных частей исходных и восстановленных сигналов. Мнимая часть сигнала также прекрасно реконструирована.

plot(ts,real(x))
hold on
plot(ti,real(is),'--')
xlim([0 0.5])
xlabel('Time (s)')
ylabel('Amplitude (V)')
legend('Original Signal','Reconstructed Signal') 
hold off

Сгенерируйте синусоиду, отобранную с частотой дискретизации 2 кГц в течение 1 секунды.

fs = 2e3;
t = 0:1/fs:1-1/fs;
x = 5*sin(2*pi*10*t);

Проектируйте периодическое окно Хэмминга длины 120. Проверьте ограничение COLA для окна с перекрытием 80 выборок. Комбинация оконного перекрытия совместима с COLA.

win = hamming(120,'periodic');
noverlap = 80;
tf = iscola(win,noverlap)

tf = logical
   1

Установите длину БПФ равной 512. Вычислите кратковременное преобразование Фурье.

fftlen = 512;
s = stft(x,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);

Вычислите обратное кратковременное преобразование Фурье.

[X,T] = istft(s,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen,'Method','ola','ConjugateSymmetric',true);

Постройте график исходных и восстановленных сигналов.

plot(t,x,'b')
hold on
plot(T,X,'-.r')
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Original and Reconstructed Signal')
legend('Original Signal','Reconstructed Signal')
hold off