Загрузите сигнал и получите преобразование константы Q.

load noisdopp
cfs = cqt(noisdopp);

Загрузите реальный сигнал и получите преобразование константы Q. Верните приблизительные частоты центра полосы пропускания.

load handel
[cfs,f] = cqt(y);

Постройте график на логарифмической шкале частот центра полосы пропускания через частоту Найквиста.

lfreq = length(f);
nyquistBin = floor(lfreq/2)+1;
plot(f(1:nyquistBin))
title('Bandpass Center Frequencies')
grid on
set(gca,'yscale','log')

Чтобы подтвердить, что отношения последовательных пар частот постоянны, постройте графики коэффициентов. Начиная с cqt использует 12 интервалов на октаву по умолчанию, отношение должно равняться $$2^{1/12}$$. Поскольку частоты DC и Nyquist не являются представителями геометрической последовательности центральных частот, но включены в вектор частот, исключить их из графика.

figure
plot(f(3:nyquistBin-1)./f(2:nyquistBin-2))
grid on
title(['Ratio: ',num2str(2^(1/12))])

Получите минимально избыточное преобразование Q аудиосигнала. Используйте окно Блэкмана-Харриса в качестве функции прототипа для систем координат Габора.

load handel
df = Fs/numel(y);
[cfs,f,g,fshifts,fintervals,bw] = cqt(y,'SamplingFrequency',Fs,'TransformType',"sparse",'Window',"blackmanharris");

cfs является массивом ячеек, где каждый элемент массива соответствует центральной частоте полосы пропускания и системе координат Габора. Постройте график системы координат Габора, сопоставленный с частотой Найквиста.

lf = length(f);
ind = floor(lf/2)+1;
gFrame = fftshift(g{ind});
fvec = f(ind-1):df:f(ind+1)-df;
plot(fvec,gFrame)
xlabel('Frequency (Hz)')
grid on
title({['Gabor Frame - Freq: ',num2str(f(ind)),' Hz'];['Bandwidth ',num2str(bw(ind)*Fs/numel(y)),' Hz']})

В преобразовании constant-Q системы координат Габора применяются к дискретному преобразованию Фурье входного сигнала, и выполняется обратное дискретное преобразование Фурье. k-й системы координат Габора применяется к k-ому частотному интервалу, указанному в fintervals. Возьмите дискретное преобразование Фурье сигнала и постройте график его величины спектра. Использование fintervals чтобы указать, над какими коэффициентами Фурье применяются системы координат Габора, сопоставленные с частотой Найквиста.

yDFT = fft(y);
lyDFT = length(yDFT);
plot(Fs*(0:lyDFT-1)/lyDFT,abs(yDFT))
grid on
fIntervalGabor = fintervals{ind};
mx = max(abs(yDFT));
hold on
plot([df*fIntervalGabor(1) df*fIntervalGabor(1)],[0 mx],'r-','LineWidth',2)
plot([df*fIntervalGabor(end) df*fIntervalGabor(end)],[0 mx],'r-','LineWidth',2)
str = sprintf('Gabor Frame Interval (Hz): [%3.2f, %3.2f]',df*fIntervalGabor(1),df*fIntervalGabor(end));
title(str)

Отобразите коэффициенты Фурье в интервале с системой координат Габора и примите обратное дискретное преобразование Фурье. Нормализуйте результат и сравните с вычисленными коэффициентами константы Q и подтвердите, что они равны.

lGframe = length(gFrame);
indx = 1:lGframe;
indx = fftshift(indx);
winDFT(indx) = yDFT(fIntervalGabor).*fftshift(gFrame(indx));
cqCoefs = ifft(winDFT);
cqCoefs = (2*lGframe/length(y))*cqCoefs;
max(abs(cqCoefs(:)-cfs{ind}(:)))

ans = 0

Загрузите аудиосигнал. Постройте график преобразования constant-Q (CQT) с помощью максимально избыточной версии преобразования и с использованием 12 интервалов на октаву.

load handel
cqt(y,'SamplingFrequency',Fs)

Выполните CQT того же сигнала, используя 48 интервалов на октаву. Установите частотную область значений, в котором CQT имеет логарифмическую частотную характеристику, как минимально допустимую частоту, равную 2 кГц.

minFreq = Fs/length(y);
maxFreq = 2000;
figure
cqt(y,'SamplingFrequency',Fs,'BinsPerOctave',48,'FrequencyLimits',[minFreq maxFreq])