Загрузите сигнал и получите преобразование константы-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 ячеек на октаву по умолчанию, отношение должно быть равно $${\mathrm{}}^{\mathrm{}\mathrm{}\mathrm{21/12}}$$. Поскольку частоты постоянного тока и Найквиста не являются членами геометрической последовательности центральных частот, но включены в частотный вектор, исключите их из графика.

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']})

В преобразовании константы-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

Загрузите звуковой сигнал. Постройте график преобразования константы-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])