DCT变换的原理及算法.docx
《DCT变换的原理及算法.docx》由会员分享,可在线阅读,更多相关《DCT变换的原理及算法.docx(31页珍藏版)》请在冰豆网上搜索。
DCT变换的原理及算法
DCT变换的原理及算法
离散傅立叶变换(DiscreteFourierTransform)
离散傅立叶变换概述
傅立叶分析以法国数学家和物理学家JeanBaptisteJosephFourier命名,是一种将信号分解为谐波的方法。
如下三图所示,一个包含16个点的离散信号可以用9个余弦和9个正弦波来表示。
在表达任意一个离散信号时,这些三角波的周期是一定的,不同的只是振幅(amplitude)。
图1-1离散信号与对应的三角波
信号可以是连续的或离散的,同时也可以是周期性的或非周期性的,根据信号的这两个特点,傅立叶变换可以分为四种类型:
傅立叶变换(FourierTransform),处理非周期性的连续信号(Aperiodic-
Continuous)。
傅立叶序列(FourierSeries),处理周期性的连续信号(Periodic-
Continuous)。
离散时间域傅立叶变换(DiscreteTimeFourierTransform),处理非周期性的离散信号(Aperiodic-Discrete)。
离散傅立叶变换(DiscreteFourierTransform),处理周期性的离散信号
(Periodic-Discrete)。
计算机只能处理离散的和有限长度的信号,因此只有离散傅立叶变换(DFT)能在计算
机中以算法实现。
图1-2四种不同类型的傅立叶分析
实数离散傅立叶变换(RealDFT)的格式和表示
如图1-3所示,离散傅立叶变换将包含N个点的输入波转为两个包含N/2+1个点的输出波。
输入波常被称作时间域,因为信号的波形基本上都是随时间变化,输出波常被称作频率域。
时间域与频率域中存储的信息是一样的,只是表现方式不一样。
将时间域转为频率域的过程叫离散傅立叶变换(DFT),将频率域转换为时间域的过程叫反变换(IDFT)。
频率域可以分为两部分,实数部分ReX[]和虚数部分ImX[],分别存放余弦函数(Cosine)的振幅和正弦函数(Sine)的振幅。
图1-3实数傅立叶变换示意图
DFT基函数
DFT中使用的正弦和余弦函数统称为基函数(BasisFunction),这些三角函数的周
期是固定的,变化的只是振幅。
DFT基函数的表达式:
Ck[i]=cos(2πki/N)
Sk[i]=sin(2πki/N)
公式1-1
其中,Ck[i]和Sk[i]表示由N个点组成的离散正弦曲线,i的取值范围是张倒N-1。
k
决定了曲线的周期,取值范围是0到N/2。
多余的系数
完成DFT后,系数由原来的N个变为N+2个,似乎产生了两个多余的系数。
在频率域中,的确有两个系数是多余的,它们是ImX[0]和ImX[n/2]。
它们的存在使得频率域中的其他系数相互独立,并且它们的值永远为0,因此不会影响反变换。
反变换的计算(IDFT)
公式1-2
在上面的公式中,振幅使用的是
和
,而不是ImX[k]和ReX[k]。
两者的关系可以用下面的公式来表示:
公式1-3
反变换的算法实现
下面是反变换算法的伪代码实现:
100'THEINVERSEDISCRETEFOURIERTRANSFORM
110'Thetimedomainsignal,heldinXX[],iscalculatedfromthefrequencydomainsignals,
120'heldinREX[]andIMX[].130'
140DIMXX[511]'XX[]holdsthetimedomainsignal
150DIMREX[256]'REX[]holdstherealpartofthefrequencydomain160DIMIMX[256]'IMX[]holdstheimaginarypartofthefrequency
domain
170'
180PI=3.14159265'Settheconstant,PI
190N%=512'N%isthenumberofpointsinXX[]200'
210GOSUBXXXX'MythicalsubroutinetoloaddataintoREX[]and
IMX[]
220'
230
240'
'FindthecosineandsinewaveamplitudesusingEq.1-3250FORK%=0TO256
260REX[K%]=REX[K%]/(N%/2)
270IMX[K%]=-IMX[K%]/(N%/2)
280NEXTK%
290'
300REX[0]=REX[0]/2
310REX[256]=REX[256]/2
320'
330'
340FORI%=0TO511
'ZeroXX[]soitcanbeusedasanaccumulator350XX[I%]=0
360NEXTI%
370'
380'Eq.1-2SYNTHESISMETHOD#1.Loopthrougheach
390'frequencygeneratingtheentirelengthofthesineandcosine400'waves,andaddthemtotheaccumulatorsignal,XX[]
410'
420FORK%=0TO256'K%loopsthrougheachsampleinREX[]and
IMX[]
430FORI%=0TO511'I%loopsthrougheachsampleinXX[]440'
450XX[I%]=XX[I%]+REX[K%]*COS(2*PI*K%*I%/N%)
460XX[I%]=XX[I%]+IMX[K%]*SIN(2*PI*K%*I%/N%)
470'
480NEXTI%
490NEXTK%
500'
510END
图1-4IDCT示意图
图1-4解释了IDCT的过程以及频率域与反变换时振幅的不同。
图1-4a中是我们需要进行转换的时间域曲线,在0点坐标处振幅为32的一条曲线;图1-4b是进行DFT变换后的曲线,实数部分(ReX)的值为32,虚数部分的值全部为0,因此这里没有画虚数部分的曲线;公式1-3将频率域信号(图1-4b所示)转换为余弦函数的振幅(图1-4c所示)。
虚数部分(用正弦函数表示)的系数全部为0,因此这里没有显示。
频率域与三角函数振幅之所以不同,是因为频率域被定义为谱密度(spectraldensity)。
图1-4显示的是一个由32个点组成的信号在频率域中的实数部分,由编号从0到16的17个采样点组成。
谱密度是指单位带宽所能表达的信号量(振幅),其计算方法是使用三角函数的振幅除以相应的带宽(bandwidth)。
图1-4中的虚线解释了带宽的计算,即在采样点之间进行平均分割。
采样点0和16的带宽为1/N,其他采样点(1-15)的带宽为2/N。
这就是公式1-3中ReX[0]和ReX[N/2]与其他ReX不同的原因。
图1-5频率域的带宽(bandwidth)
DFT的分析与计算
DFT的计算有三种方法:
第一种是解线性方程祖,这种方法简单但计算量很大,实际应用中很少使用;第二种是关联法(correlation),基于已知的另一条的曲线;第三中方法是快速傅立叶变换(FFT),FFT算法将对一条含N个点曲线的计算转为对N条含1个点的曲线的计算,从而大大降低了计算量。
线性方程组求解DFT
使用这种方法来计算DFT是很自然的,从N个方程求解N个未知数。
但这种方法计算量很大,实际应用中很少使用。
关联法(correlation)求解DFT
通过correlation求解DFT是求解DFT的标准方法,下面使用一个例子来说明这种方法。
求解一个包含64个点的信号的DFT,意味着我们要计算频率域中实数部分的33个系数和虚数部分的33个系数。
在这个例子中,我们只解一个系数,ImX[3]。
ImX[3]是一条包含三个完整周期的正弦函数的振幅,这个正弦函数曲线分布在点0到点63之间,如图1-6a所示。
在图1-6中,a和b是两个时间域的示例信号,在这里分别称之为x1[]和x2[]。
曲线x1[]是一个包含3个完整周期、分布在0到63之间的正弦函数曲线;x2[]由多条正弦函数曲线和余弦函数曲线混合而成,在组成x2[]的三角函数曲线中,没有任何一条在0到63之间有完整的三个周期。
通过这两条曲线可以解释算法要实现的功能:
当输入函数为x1[]时,算法的计算结果应该是32,也就是信号中正弦曲线的振幅;而当输入函数为x2[]时,算法的计算结果应该为0,因为x2[]所表示的曲线并不在这个信号中。
之所以非x1[]的曲线与x1[]相乘结果为0,是因为除x1[]外的任意一个函数与正弦
函数相乘时,在0到63(3个完整周期)上的积分都为0。
在图1-6中,图e是图a和图c相乘的结果,将各个点的值相加即可得到32;图f是图b和图d相乘的结果,将各个点的值相加得到的结果为0。
上面的计算过程可以用公式1-4来表示。
公式1-4
图1-6关联法(correlation)求DFT
下面是关联法计算DFT的算法伪代码:
100'THEDISCRETEFOURIERTRANSFORM
110'Thefrequencydomainsignals,heldinREX[]andIMX[],arecalculatedfrom
120'thetimedomainsignal,heldinXX[].130'
140DIMXX[511]'XX[]holdsthetimedomainsignal
150DIMREX[256]'REX[]holdstherealpartofthefrequencydomain160DIMIMX[256]'IMX[]holdstheimaginarypartofthefrequency
domain
170'
180PI=3.14159265'Settheconstant,PI
190N%=512'N%isthenumberofpointsinXX[]200'
210GOSUBXXXX'MythicalsubroutinetoloaddataintoXX[]220'
230'
240FORK%=0TO256'ZeroREX[]&IMX[]sotheycanbeusedas
accumulators
250REX[K%]=0
260IMX[K%]=0
270NEXTK%
280'
290''CorrelateXX[]withthecosineandsinewaves,Eq.8-4300'
310FORK%=0TO256'K%loopsthrougheachsampleinREX[]and
IMX[]
320FORI%=0TO511'I%loopsthrougheachsampleinXX[]330'
340REX[K%]=REX[K%]+XX[I%]*COS(2*PI*K%*I%/N%)
350IMX[K%]=IMX[K%]-XX[I%]*SIN(2*PI*K%*I%/N%)
360'
370NEXTI%
380NEXTK%
390'
400END
二元性(Duality)
DFT和IDFT变换公式很类似。
从一个域到另一个域,都是用已有的值乘以基函数,然后将相应的值相加。
实际上,DFT和IDFT的区别仅仅是时间域中包含N个点,而频率域中包含N/2+1个点。
在复数傅里叶变换中,时间域和频率域中的信号都包含N个点,这使得两个域具有一种对称的性质,而在两个域之间的转换公式也就几乎一样了。
时间域和频率域的这种对称被称之为对称性(Duality),二元性可以产生很多有趣的性质。
例如:
在频率域中的一个单点对应于时间域中的一条三角函数曲线,由于Duality,这种性质反过来也成立,在时间域中的一个点对应于频率域中的一条三角函数曲线。
快速傅立叶变换(FastFourierTransform)
FFT算法是由J.W.Cooley和J.W.Tukey在论文”AnalgorithmforthemachinecalculationofcomplexFourierSeries”中提出的。
FFT是基于ComplexDFT来实现的。
通过ComplexDFT来计算RealDFT
尽管FFT算法是基于ComplexDFT实现的,但我们仍可以用其来计算RealDFT,因为RealDFT可以方便地转换为ComplexDFT。
从图2-1中可以看出RealDFT和ComplexDFT的区别。
在RealDFT中,时间域是一个包含N个点的信号,频率域则包括实数部分和虚数部分两个长度为N/2+1的信号;在ComplexDFT中,时间域也有两个部分,分别是实数部分和虚数部分,长度为N。
频率域的实数部分和虚数部分则长度增至N。
如图2-1所示,RealDFT和ComplexDFT的区别仅在于后者在时间域增加了一
个虚数部分,频率域长度的变化正是由这个虚数部分引起的。
图2-1RealDFT和ComplexDFT的区别
产生这个区别的原因是实数的虚数部分为0,因此将实数表达为虚数很简单,加上一个系数为0的虚数部分即可。
例如在图2-1中的ComplexDFT,若将时间域的虚数部分设为0,频率域中多出的部分也置为0,那么图2-1中的RealDFT和ComplexDFT就相等了。
当包含负频率时,DFT的频率域会具有周期性。
在ComplexDFT中,
频率域中0到N/2为正频率,N/2+1到N-1为负频率。
与使用ComplexDFT计算RealDFT相比,使用ComplexInverseDFT计算RealInverseDFT更为困难。
这是因为频率域中N/2+1到N-1部分的系数需要计算。
其计算过程也不复杂,系数N/2+1对应系数N/2-1的相反数,N/2+2对应N/2-2的相反数,即:
系数(N/2+1)=—系数(N/2-1)系数(N/2+2)=—系数(N/2-2)
注意,0与N/2没有相应的点与之对应。
进行RealInverseDFT计算时,首先将0到N/2复制到complexDFT的系数0到N/2,然后使用一个子过程来计算系数N/2+1到N-1。
这个子过程的伪代码实现如下:
6000'NEGATIVEFREQUENCYGENERATION
6010'Thissubroutinecreatesthecomplexfrequencydomainfromtherealfrequencydomain.
6020'Uponentrytothissubroutine,N%containsthenumberofpointsinthesignals,and
6030'REX[]andIMX[]containtherealfrequencydomaininsamples0toN/2.
6040'Onreturn,REX[]andIMX[]containthecomplexfrequencydomaininsamples0toN-1.
6050'
6060FORK%=(N%/2+1)TO(N%-1)
6070REX[K%]=REX[N%-K%]
6080IMX[K%]=-IMX[N%-K%]
6090NEXTK%
6100'
6110RETURN
FFT的实现原理
FFT算法很复杂,本文不讨论细节,只描述其实现原理。
在虚数域中,时间域和频率域表达的都是由N个虚数点组成的信号,每个虚数点都由实数部分和虚数部分的两个数字来表达。
例如虚数点X[6],就是由实数部分ReX[6]和虚数部分ImX[6]组成。
FFT算法的核心思想是将时间域中一个包含N个点的信号分解为N个包含一个点的信号。
然后分别计算这N个信号的频率域对应值,最后将这N个频率域的信号综合为频率域中的一个信号。
图2-2描述了一个包含12个点的示例信号在FFT中的分解过程。
图2-2FFT中的分解(decompose)过程
图2-2中的过程看似复杂,实际上可以通过如图2-3所示的位反转算法(bitreversalsorting)来实现。
算法将各点的二进制位反转为对称的形式,即可完成N个点的信号到N个单点信号的分解过程。
图2-3位反转排序
FFT算法的下一步是分别求出这N个单点信号在频率域的振幅。
这是算法中最容易
的一步,单点的振幅等于它自己本身的值,这意味着在这一步什么也不必做。
算法的最后一步是将这N个频率域的点按在时间域分解时的反序结合(combine)起
来,这里不能使用位反转算法,这一步是算法中最复杂的部分。
图2-4展现了两个长度为4的频率域信号组合为一个长度为8的频率域信号的过
程。
组合(synthesis)的顺序必须与在时间域中分解(decompose)的过程完全相逆。
以时间域的信号abcd和信号efgh为例,要将其整合为一个包含8个点的信号需要经过这两步:
首先将这两个信号进行稀释(dilute),即用0填充为长度为8的信号,然后两者相加即可得到新的信号。
如abcd稀释后得到a0b0c0d0,efgh稀释后得到0e0f0g0h,两者相加可得abcdefgh。
图2-4FFT组合(synthesis)
当时间域用0稀释时,对应的频率域会复制自己。
当时间域先移位再用0填充时,对应的频率域会乘以一个三角函数,然后再复制自己。
abcd与efgh的稀释方法并不相同,abcd稀释为a0b0c0d0,其偶数位为0;efgh稀释为0e0f0g0h,其奇数位为0。
也就是说efgh向右移动了一位,这个在时间域的移位对应于频率域乘以一个三角函数。
图2-5展示了在两个频率域中长度为4的信号组合的过程。
左侧的Odd-FourPointFrequencySpectrum指的是对应奇数位为0的时间域信号的频率域信号,如EFGH;右侧的Even-FourPointFrequencySpectrum指的是对应偶数位为0的时间域信号的频率域信号,如ABCD。
为了更清楚地表达这个过程,图2-6将其中的两个点拿出来,因为这个图形很想一只张开翅膀的蝴蝶,因此人们也将这个图所代表的过程称之为butterfly。
图2-5FFT组合过程
图2-6Butterfly
图2-7显示了FFT变换的流程图,包含了FFT变换的三个部分。
1.时间域的分解过程可以通过位变换算法来实现。
2.将时间域分解后得到的N个单点转换为频率域并不需要任何计算,因为对于单点而言,在频率域的振幅等于时间域的振幅。
3.第三部分是整个算法的核心,是图中重点要表达的部分。
图2-7FFT流程图
在图2-7中,最外面的循环表示要在lgN个层次上进行组合(synthesis),中间那层循环指在每一层上的组合过程,最内部的循环表示butterfly过程。
下面是FFT算法的一段Basic代码:
1000'THEFASTFOURIERTRANSFORM
1010'Uponentry,N%containsthenumberofpointsintheDFT,REX[]and1020'IMX[]containtherealandimaginarypartsoftheinput.Uponreturn,1030'REX[]andIMX[]containtheDFToutput.Allsignalsrunfrom0toN%-1.1040'
1050PI=3.14159265'Setconstants1060NM1%=N%-1
1070ND2%=N%/2
1080M%=CINT(LOG(N%)/LOG
(2))
1090J%=ND2%
1100'
1110FORI%=1TON%-2'Bitreversalsorting1120IFI%>=J%THENGOTO1190
1130TR=REX[J%]
1140TI=IMX[J%]
1150REX[J%]=REX[I%]
1160IMX[J%]=IMX[I%]
1170REX[I%]=TR
1180IMX[I%]=TI
1190K%=ND2%
1200IFK%>J%THENGOTO1240
1210J%=J%-K%
1220K%=K%/2
1230GOTO1200
1240J%=J%+K%
1250NEXTI%
1260'
1270FORL%=1TOM%'Loopforeachstage1280LE%=CINT(2^L%)
1290LE2%=LE%/2
1300UR=1
1310UI=0
1320SR=COS(PI/LE2%)'Calc