CQT-采自http://read.pudn.com/downloads166/sourcecode/multimedia/audio/764229/genlgftkern.m__.htm

const Q trans

1 % genlgftkern 7 / 99 write this in the time domain
2   % write a genkern like genkern.c for a const Q trans (?? genkern was for transform from fft)
3   % %
4   % % Example: [kerncos, kernsin, freqs ] = genlgftkern( 174.6141 , 1.0293022366 , 11025 , 120 , 1102 );
5   %
6   % nfreqs and windsizmax optional but must input 0 if not a value
7   % will choose 100 ms max and nfreqs up to Nyq
8   % minfreq = 440 / ( 2 ^ ( 1 / 12 )) ^ 16 ; % F3 174.6141 G3 = 195.9977
9   % nfreqs = 60 ;
10   % freqrat = 2 ^ ( 1 / 12 ); % 2 ^ ( 1 / 24 ) = 1.0293022366
11 % SR = 11025 ;
12 % windsizmax = fix(. 1 * SR); % take a 100 ms max
13
14 function [kerncos, kernsin, freqs] = genlgftkern(minfreq, freqrat, SR, nfreqs, windsizmax, winnam, constwind);
15
16 if (exist( ' constwind ' ) & constwind ~= 0 ), fprintf( ' genlgftkern; Calculating kernels with constant windowsiz = %.0f\n ' , constwind);
17 elseif ~ exist( ' constwind ' ), constwind = 0;
18 end
19 % Put these steps back in if want to run this prog indep of logftS
20 % if (nargin < 5 | windsizmax == 0) % no input windsizmax
21 % if windsizmax == 0
22 % windsizmax = floor( . 1 * SR); % take a 100 ms max
23 % fprintf( ' No input maximum window size. Taking 100 ms max.\n ' );
24 % end
25 % if (nargin < 4 | nfreqs == 0) % no nfreqs either
26 % if nfreqs == 0
27 % nfreqs = fix( log(SR / ( 2 * minfreq)) / log(freqrat) ); % to give highest freq at the Nyq
28 % fprintf( ' No input number of freqs; taking freqs from minfreq to Nyquist = %.0f freq bins\n ' , nfreqs);
29 % end
30
31 % winnam = ' hamming ' ;
32 if ~ exist( ' winnam ' ), winnam = ' hamming ' ;
33 fprintf( ' Using default window Hamming\n ' );
34 else , fprintf( ' Input window %s \n ' , winnam);
35 end
36 if winnam( 1 : 4 ) == ' rect ' , winnam = ' boxcar ' ; end
37
38 Q = 1 / (freqrat - 1 );
39 TWOPI = 2 * pi;
40 mindigfreq = TWOPI * minfreq / SR;
41 freqs = minfreq * (freqrat . ^ [(0: 1 :nfreqs - 1 )]);
42 pos = find(freqs < SR / 2 );
43 freqs = freqs(pos);
44 nfreqs = length(freqs);
45 digfreq = freqs * TWOPI / SR;
46 % shouldn ' t need the following since fixed up freqs
47 if sum(find(digfreq > pi)) ~= 0, error( ' freq over Nyq ' ); end
48
49 flag = 1 ;
50 if constwind == 0
51 windsizOk = fix (TWOPI * Q . / digfreq); % period in samples time Q
52 % arg = (pi / 2 ) * ones(nfreqs, windsiz);
53 windsizmax;
54 windsizOk( 1 );
55 if (windsizmax > windsizOk( 1 )),
56 windsizmax = windsizOk( 1 );
57 flag = 0;
58 end
59 pos = find(windsizOk > windsizmax);
60 % if windsizmax < windsizOk( 1 ) so get some windows not as long as necess for that Q
61 if (flag)
62 fprintf( ' genlgftkern: Const window %.0f up to freq position %.0f and frequency %.0f out of %.0f frequencies. windsizMinfreq = %.0f Q=%.1f\n ' , ...
63 windsizmax, max(pos), digfreq(max(pos)) * SR / ( 2 * pi), length(digfreq), windsizOk( 1 ), Q);
64 else
65 fprintf( ' genlgftkern: No const window; windsizmax = %.0f = windsizMinfreq = %.0f Q=%.1f\n ' , ...
66 windsizmax, windsizOk( 1 ), Q);
67 end
68 fprintf( ' \n ' );
69 windsizOk(pos) = windsizmax;
70
71 kerncos = zeros(nfreqs, windsizOk( 1 ) );
72 kernsin = zeros(nfreqs, windsizOk( 1 ) );
73 numzeros = windsizOk( 1 ) - windsizOk;
74 numzerosO2 = round(numzeros / 2 );
75
76 else
77 kerncos = zeros(nfreqs, constwind );
78 kernsin = zeros(nfreqs, constwind );
79 end
80
81 % Get kaiser number if window is kaiser
82 if length(winnam) > 5
83 if winnam( 1 : 6 ) == ' kaiser '
84 if length(winnam) == 7 , kaiserno = winnam( 7 );
85 elseif length(winnam) == 8 , kaiserno = winnam( 7 : 8 );
86 else , kaiserno = ' 8 ' ; % default is 8 for no input kaiser number
87 end
88 winnam = ' kaiser ' ;
89 end
90 end
91
92 if constwind == 0,
93 for k = 1 :nfreqs
94 sz = windsizOk(k);
95 switch(winnam)
96 case ' kaiser '
97 winstr = [ winnam ' ( ' num2str(sz) ' , ' kaiserno ' ) ' ];
98 otherwise
99 winstr = [ winnam ' ( ' num2str(sz) ' ) ' ];
100 end
101 % fprintf( ' calc window %s \n ' , winstr);
102 wind = eval(winstr);
103 wind = wind ' ;
104
105 % wind = boxcar(windsizOk(k)) ' ;
106 numz = 1 ;
107 if numzerosO2(k) ~= 0, numz = numzerosO2(k); end;
108 kerncos(k, numz: numz + windsizOk(k) - 1 ) = ( 1 / windsizOk(k)) * ...
109 cos(digfreq(k) * ( - sz / 2 : sz / 2 - 1 )). * wind;
110 % cos(digfreq(k) * (0:windsizOk(k) - 1 )). * wind;
111
112 % cos(digfreq(k) * (0:windsizOk(k) - 1 )). * wind(( 1 :windsizOk(k)));
113 kernsin(k, numz: numz + windsizOk(k) - 1 ) = ( 1 / windsizOk(k)) * ...
114 sin(digfreq(k) * ( - sz / 2 : sz / 2 - 1 )). * wind;
115 % sin(digfreq(k) * (0:windsizOk(k) - 1 )). * wind;
116 % sin(digfreq(k) * (0:windsizOk(k) - 1 )). * wind(( 1 :windsizOk(k)));
117 end
118 else
119 for k = 1 :nfreqs
120 sz = constwind;
121 switch(winnam)
122 case ' kaiser '
123 winstr = [ winnam ' ( ' num2str(sz) ' , ' kaiserno ' ) ' ];
124 otherwise
125 winstr = [ winnam ' ( ' num2str(sz) ' ) ' ];
126 end
127 % fprintf( ' calc window %s \n ' , winstr);
128 wind = eval(winstr);
129 wind = wind ' ;
130
131 % wind = boxcar(windsizOk(k)) ' ;
132
133 kerncos(k, 1 : constwind) = (cos(digfreq(k) * ( - constwind / 2 : constwind / 2 - 1 ))). * wind;
134 % kerncos(k, 1 : constwind) = (cos(digfreq(k) * (0:constwind - 1 ))). * wind;
135
136 % cos(digfreq(k) * (0:windsizOk(k) - 1 )). * wind(( 1 :windsizOk(k)));
137 kernsin(k, 1 : constwind) = (sin(digfreq(k) * ( - constwind / 2 : constwind / 2 - 1 ))). * wind;
138 % kernsin(k, 1 :constwind) = (sin(digfreq(k) * (0:constwind - 1 ))). * wind;
139 end
140 end
141
142 % fprintf( ' nfreqs = %.0f ; windsizmax = %.0f \n ' ,nfreqs, windsizmax);
143 ...

http://www.elec.qmul.ac.uk/people/anssik/cqt/smc2010.pdf

http://www.elec.qmul.ac.uk/people/anssik/cqt/

CQT

CQT把时域信号转换到时频域，使得频率下标的中心频率是成几何间隔排列的，并且Q因子相等。这意味着在低频域有这很好的频率分辨率，在高频域有好的时间分辨率。

CQT本质上是一种小波变换，但是我们考虑相对高Q因子的变换，等效于每八度音阶12-96个下标。很多传统的小波变换不满足这些要求，例如基于迭代滤波器组的方法要求对输入信号进行几百次滤波。

西方音乐
十二平均律

F0s   



 
 
 
 
 
 
 
 
 
 
 
 

从听觉的角度看，从20khz到500hz，外围听觉系统的频率分辨率是近似常Q的，低于这个范围Q值减小。从感知音频编码的角度看，最短的变换窗长必须是3ms为了保证高的质量，然而，在低频率编码需要更高的频率分辨率。

传统的DFT是线性间隔的频率下标，因此在可听见频率的范围内，不能满足变化的时间和频率分辨率的要求。

三个关键问题，使得目前在音频信号处理中，CQT没有广泛地代替DFT

1. 
   高的计算代价
2. 
   没有逆变换能够从变换系数，完美地重建原始信号
3. 
   在连续的时间帧里很难处理时频矩阵。不同的频率下标的抽样是不同步的。

Brown and Puckette提出一个计算有效的带有高品质因子的常Q变换，基于FFT和一个频域核。CQT实现的缺点是没有逆变换。

FitzGerald提出一个好的近似逆变换的，如果反转的信号在DFT域有稀疏的表达。但，一般地说，这对乐音信号是不准确的。

           k=1,2,…，K

是下标k的中心频率

指示了采样率

是连续窗函数，在

范围外是0.

最低的中心频率，B每八度音阶的下标数目，是CQT最重要的参数，确定了时频分辨率

是原子

频率响应的-3db带宽，并且

是窗函数频谱主瓣-3db带宽，

。

理想的Q尽可能大，为了使带宽

尽可能小，因此引入了最小频率拖尾。

是一个尺度因子，典型地q=1.

q值越小，则时间分辨率越好，但频率分辨率降低。

q=0.5和B=48与q=1和B=24获得一样的时频分辨率，但是前者每八度音阶包含两倍的频率下表。q<1可以堪称频率轴的过采样，类似于DFT时的补零。例如q=0.5对应于过采样因子2，有效的频率分辨率等于每八度音阶B/2个下标，尽管计算了B个下标。

与

无关

为了达到合理的信号重建，

3.1 Algorithm of Brown and Puckette

X(j)是x(n)的DFT,并且A(j)是a(n)的DFT.Parsecal’s
theorem

A_k(j)是n点a_k(n)变换基的DFT，ak(n)是N/2变换帧的中心。A_k(j)是谱核，是稀疏的，作为A的一列。 a_k(n)是时域核。

转载于:https://www.cnblogs.com/ttschina/archive/2011/05/28/2058504.html