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







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







posted @ 2011-05-28 23:33  ttschina  阅读(824)  评论(0)    收藏  举报