• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于解析插值離散時(shí)間傅里葉變換的精確頻率估計(jì)

    2022-04-08 05:42:28胡緒權(quán)付志紅
    電工技術(shù)學(xué)報(bào) 2022年6期
    關(guān)鍵詞:信號(hào)

    楊 超 李 波,2 胡緒權(quán),3 付志紅

    基于解析插值離散時(shí)間傅里葉變換的精確頻率估計(jì)

    楊 超1李 波1,2胡緒權(quán)1,3付志紅1

    (1. 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué)) 重慶 400044 2. 云南電網(wǎng)有限責(zé)任公司電力科學(xué)研究院 昆明 650217 3. 重慶璀陸探測(cè)技術(shù)有限公司 重慶 402660)

    短時(shí)波動(dòng)周期下,來自負(fù)頻率的頻譜泄漏是頻域插值算法估計(jì)誤差的主要來源。為消除頻譜泄漏影響以實(shí)現(xiàn)任意波動(dòng)周期下實(shí)正弦信號(hào)頻率的精確估計(jì),該文提出一種基于廣義離散時(shí)間傅里葉變換(DTFT)插值解析的頻率估計(jì)算法。該算法利用位于頻譜主瓣范圍內(nèi)的任意DTFT譜線構(gòu)建對(duì)應(yīng)實(shí)/虛部的線性方程組,進(jìn)而通過方程組求解得到基于解析插值DTFT的無偏頻率估計(jì)公式。此外,考慮實(shí)際中加性白噪聲的影響,根據(jù)頻率估計(jì)的統(tǒng)計(jì)特性分析給出一種簡(jiǎn)單且有效的迭代方法,從而確保頻率估計(jì)的方均誤差(MSE)在任意波動(dòng)周期下能夠一致逼近無偏頻率估計(jì)的克拉美羅下屆(CRLB)。仿真和實(shí)驗(yàn)結(jié)果驗(yàn)證了所提算法優(yōu)于現(xiàn)有的DFT插值算法。

    短時(shí)波動(dòng)周期 頻譜泄漏 解析插值DTFT 迭代方法 克拉美羅下屆(CRLB)

    0 引言

    正弦信號(hào)的頻率估計(jì)是一個(gè)經(jīng)典且重要的信號(hào)處理問題,在科學(xué)和工程領(lǐng)域具有極為廣泛的應(yīng)用。在電力系統(tǒng)中,頻率作為最重要的電能質(zhì)量參數(shù)之一,其變化反映了發(fā)電和負(fù)載之間的動(dòng)態(tài)平衡[1-3],通過頻率及其變化的準(zhǔn)確估計(jì)可以實(shí)現(xiàn)系統(tǒng)的動(dòng)態(tài)過程監(jiān)測(cè)。此外,在計(jì)量基準(zhǔn)測(cè)試和標(biāo)準(zhǔn)溯源等特定應(yīng)用中,對(duì)頻率估計(jì)精度往往具有更高的要求。

    頻率估計(jì)算法總體上可分為時(shí)域和頻域兩類,頻域算法主要基于由快速離散傅里葉變換(Fast Discrete Fourier Transform,FDFT)實(shí)現(xiàn)的插值離散傅里葉變換(Interpolated Discrete Fourier Transform, IpDFT),在計(jì)算復(fù)雜度方面相對(duì)于通常需要矩陣求逆的時(shí)域方法更具優(yōu)勢(shì)(時(shí)域?yàn)?3),頻域?yàn)?log2)[4]),因而在實(shí)時(shí)頻率估計(jì)應(yīng)用中表現(xiàn)出特別的吸引力。但是,過去大多數(shù)針對(duì)IpDFT頻率估計(jì)的研究都集中在單頻解析信號(hào)上[5-9],對(duì)于實(shí)際中應(yīng)用更為廣泛的實(shí)正弦信號(hào),受共軛負(fù)頻率的頻譜長(zhǎng)泄漏影響,尤其是在被測(cè)信號(hào)對(duì)應(yīng)波動(dòng)周期(Cycles in Record, CiR)較小時(shí),IpDFT存在嚴(yán)重的估計(jì)偏差。

    短時(shí)CiR頻譜泄漏問題已引起國內(nèi)外學(xué)者的高度關(guān)注并開展了大量研究,提出多種改進(jìn)IpDFT算法,總體上可分為直接方法和迭代方法。直接方法力求在計(jì)及負(fù)頻率頻譜長(zhǎng)泄漏基礎(chǔ)上直接得到頻率估計(jì)的解析解。文獻(xiàn)[10]提出了一種基于矩形窗的3點(diǎn)IpDFT(3IpDFT),利用變換構(gòu)建對(duì)應(yīng)譜線的線性方程組,進(jìn)而求解得到頻率估計(jì)的解析解,但仿真結(jié)果顯示,即便在最優(yōu)情況下算法方均誤差(Mean Square Error, MSE)仍然無法逼近無偏頻率估計(jì)的克拉美羅下界(Cramer-Rao Lower Bound, CRLB)。文獻(xiàn)[11-12]在文獻(xiàn)[13]解析泄漏補(bǔ)償方法基礎(chǔ)上,提出了利用信號(hào)補(bǔ)零和DFT系數(shù)變換的2點(diǎn)IpDFT(2IpDFT),但補(bǔ)零操作和繁瑣的譜線選擇判據(jù)使得算法計(jì)算復(fù)雜度相對(duì)較高。文獻(xiàn)[14-15]基于最大旁瓣衰減窗離散等間隔頻譜近似線性比例特性,以求解頻譜線性方程組的方式得到了3點(diǎn)IpDFT(3IpDFT),但鑒于窗函數(shù)的數(shù)學(xué)復(fù)雜性增加[16],事實(shí)上在高階情況下很難得到頻率估計(jì)的解析解,因而無法進(jìn)一步消除窗函數(shù)頻譜近似導(dǎo)致的估計(jì)偏差。此外,相對(duì)矩形窗而言,采用具有更高等效噪聲帶寬(Equivalent Noise Band Width, ENBW)[17]的最大旁瓣衰減窗,算法MSE在噪聲擾動(dòng)下將無法逼近CRLB。

    迭代方法則常以適用于單頻解析信號(hào)的IpDFT為基礎(chǔ),力求通過迭代不斷消除頻譜泄漏干擾,從而確保算法最終收斂于頻率估計(jì)的解析解。文獻(xiàn)[18]在最大旁瓣衰減窗頻譜近似基礎(chǔ)上,提出了一種基于頻譜系數(shù)加權(quán)的單步迭代IpDFT(3IpDFT),但以經(jīng)典IpDFT[19]作為其起始解析方式,將導(dǎo)致短時(shí)CiR下算法MSE中存在明顯的估計(jì)偏差。文獻(xiàn)[20]對(duì)AM算法[21]進(jìn)行了有效改進(jìn),提出了一種利用位于DFT譜線中間的離散時(shí)間傅里葉變化(Discrete Time Fourier Transform, DTFT)系數(shù)插值的多步迭代算法(AM)。其基本思想與2IpDFT類似,即盡量采用幅度最大的頻譜系數(shù)以降低噪聲干擾影響,進(jìn)而通過不斷重構(gòu)并消除頻譜中負(fù)頻率泄漏分量以逐漸降低頻率估計(jì)偏差,使得算法MSE最終一致收斂逼近于CRLB,是現(xiàn)有頻域估計(jì)算法中精度最高的[22]。但由于其核心兩點(diǎn)插值方法為有偏頻率估計(jì),在估計(jì)偏差占主導(dǎo)情況下需要多次計(jì)算DTFT系數(shù)并進(jìn)行迭代,這無疑大大增加了算法的計(jì)算成本。此外,采用固定0.5D(D為頻率分辨率)的DTFT系數(shù)還使得算法在CiR<1時(shí)性能將急劇下降。

    基于AM算法的基本思想,本文提出了一種基于廣義DTFT插值的精確頻率估計(jì)算法。一方面,通過準(zhǔn)確解析由任意位置的廣義DTFT譜線構(gòu)成的線性方程組,可以確保算法在任意CiR下得到無偏頻率估計(jì);另一方面,通過統(tǒng)計(jì)特性分析導(dǎo)出用以優(yōu)化算法性能的簡(jiǎn)單迭代方法,可以最大程度地抑制背景噪聲干擾以確保算法MSE能夠一致逼近CRLB,從而實(shí)現(xiàn)任意CiR下頻率的精確估計(jì)。

    1 頻率估計(jì)的解析解

    1.1 基本原理

    假設(shè)以采樣頻率s進(jìn)行等間隔均勻離散采樣,則點(diǎn)采樣序列()的數(shù)學(xué)模型可表示為

    式中,=0, 1,…,-1?N;、和分別為信號(hào)()的頻率、幅值和初始相位;()為具有零均值和方差2的加性高斯白噪聲;/s//為頻域歸一化頻率,在時(shí)域中對(duì)應(yīng)于被測(cè)信號(hào)的波動(dòng)周期個(gè)數(shù),在頻域中則表示在頻譜中所處位置。據(jù)奈奎斯特采樣定律可知,max<0.5,且無量綱,其單位通常采用bin標(biāo)示;=[]?N,符號(hào)[ ]表示四舍五入取整,由D=s/可知,非整余數(shù)?[-0.5, 0.5],則點(diǎn)離散信號(hào)()的廣義DTFT可表示為

    (2)

    式中,j為復(fù)數(shù)單位;(?)為寬帶噪聲的DTFT;第二項(xiàng)為負(fù)頻率分量引起的頻譜長(zhǎng)泄漏干擾。值得注意的是,此處參量可為任意[0,-1]范圍內(nèi)的實(shí)數(shù),當(dāng)其為整數(shù)時(shí),式(2)即為標(biāo)準(zhǔn)DFT,此時(shí)()的等間隔抽樣即可表示為(|?N)。

    便于推導(dǎo),首先忽略噪聲項(xiàng)(?),通過計(jì)算可將式(2)轉(zhuǎn)化為

    將式(2)改寫為式(3)所示方程具有兩個(gè)非常重要的作用:一方面,式(3)中含有3個(gè)未知參量、和,可知還需另外兩條DTFT譜線{(±)}以形成方程組進(jìn)行解析;更為重要的一方面在于式(3)左右側(cè)為復(fù)數(shù),自然地可分解為對(duì)應(yīng)實(shí)部和虛部的方程,這是求解方程組并導(dǎo)出解析插值頻率估計(jì)的關(guān)鍵。

    1.2 齊次線性方程組

    如式(3)左側(cè)所示,()e-jpk(1-N)/N與未知參量無關(guān),僅由()及其位置索引決定,簡(jiǎn)單起見采用 =()e-jpk(1-N)/N表示,則據(jù)式(3)可得,對(duì)應(yīng)實(shí)部的齊次線性方程組,以矩陣形式表示為

    方程組式(4)中系數(shù)矩陣定義如下

    其中

    同理,對(duì)應(yīng)虛部的齊次線性方程組為

    方程組式(9)中系數(shù)矩陣的相關(guān)向量定義如下

    在式(4)~式(12)中,需要明確指出的是?R+,其范圍控制在(0, 0.5],使得所用DTFT譜線全部位于矩形窗譜的主瓣范圍內(nèi),以確保具有高信噪比。

    1.3 頻率估計(jì)

    由于齊次線性方程組式(4)和式(9)必有非零解,則根據(jù)線性代數(shù)中克萊姆法則(Cramer’s Rule)可知方程組系數(shù)矩陣的行列式必然為零[22],即

    雖然據(jù)1.2節(jié)中各矩陣定義可知,上述方程都僅含有單一的未知變量,至此可以單獨(dú)求解式(13)或式(14)而直接得到的解析解。但是,當(dāng)所用DTFT譜線的實(shí)部或虛部接近零時(shí),噪聲干擾將導(dǎo)致估計(jì)極為不準(zhǔn)確。為防止這種情況的出現(xiàn),需結(jié)合式(13)和式(14)得到用于頻率估計(jì)的解析插值公式(具體求解過程見附錄第1節(jié)),結(jié)果如下

    其中

    2 頻率估計(jì)的統(tǒng)計(jì)特性分析

    第1節(jié)通過研究無噪聲情況,得到了基于廣義插值DTFT的頻率估計(jì)解析式,本節(jié)將分析加性高斯白噪聲背景下頻率估計(jì)算法的性能并得到估計(jì)的漸近性質(zhì),進(jìn)而通過簡(jiǎn)單迭代優(yōu)化算法性能。

    據(jù)式(23)可知,頻率估計(jì)由決定,因而首先分析噪聲對(duì)估計(jì)的影響。令=( )e-jpk(1-N)/N,則背景噪聲干擾下的估計(jì)應(yīng)準(zhǔn)確表示為

    其中

    由此可得噪聲干擾下的估計(jì)偏差D

    其中

    由此利用泰勒級(jí)數(shù)展開,式(23)所示的頻率估計(jì)解析式可改寫為

    則噪聲干擾下頻率估計(jì)的統(tǒng)計(jì)偏差為

    由白噪聲特性(1)=(2)=0可知,期望值(D)= 0,表明式(23)為頻率的無偏估計(jì)。

    進(jìn)一步考慮噪聲項(xiàng)之間的相關(guān)性,可得頻率估計(jì)的方差(詳細(xì)推導(dǎo)過程見附錄第2節(jié))為

    式中,[Re2(D )]如式(A20)所示。

    進(jìn)一步分析可知,任意頻點(diǎn)下式(31)的變化趨勢(shì)及其量級(jí)大小始終由|1-G2-G3|決定,因而通過分析|1-G2-G3|即可掌握頻率估計(jì)的統(tǒng)計(jì)方差變化。根據(jù)中對(duì)應(yīng)各定義,將式(2)和式(21)代入1-G2-G3即可簡(jiǎn)化為

    式中,=(/2)exp{j[+p(-1)/]};′=2+;()定義為

    通常情況下>1,此時(shí)中的負(fù)頻率泄漏分量可忽略不計(jì),在中心譜線位置搜索正確的前提下(即=),1-G2-G3可近似為

    根據(jù)式(34)可知:

    (1)>0。當(dāng)Re()=0和Im()=1時(shí),對(duì)應(yīng)|1-G2-G3|的最大值,此時(shí)單一頻點(diǎn)下式(31)僅包含Im(1)時(shí)接近其最小值(下界);當(dāng)Re()=1和Im()=0時(shí),對(duì)應(yīng)|1-G2-G3|的最小值,此時(shí)單一頻點(diǎn)下式(31)僅包含Re(1)時(shí)接近其最大值(上界)。

    (2)≤0。當(dāng)Re()=1和Im()=0時(shí),對(duì)應(yīng)|1-G2-G3|的最大值,此時(shí)單一頻點(diǎn)下式(31)僅包含Im(1)時(shí)接近其最小值(下界);當(dāng)Re()=0和Im()=1時(shí),對(duì)應(yīng)|1-G2-G3|的最小值,此時(shí)單一頻點(diǎn)下式(31)僅包含Re(1)時(shí)接近其最大值(上界)。

    實(shí)際中由于信號(hào)相位為隨機(jī)分布,因此取其中值作為平均估計(jì),即Re()=Im()=cos(p/2)。以= 1 024,=2,s=1 024Hz,SNR=40dB,?[2.55, 3.45]和=0.1為例,對(duì)應(yīng)每一個(gè)時(shí)頻率估計(jì)方差的理論分布如圖1所示。

    圖1 i=0.1時(shí)頻率估計(jì)方差的范圍

    圖1結(jié)果顯示,在=0.1且相位隨機(jī)分布前提下,所提頻率估計(jì)方法的方差在→0時(shí)接近最小值,且最小值十分接近無偏頻率估計(jì)的CRLB[11]。進(jìn)一步地,考慮相位在[-p,p]范圍內(nèi)隨機(jī)變化,值在0.1~0.5范圍內(nèi)逐漸增大,頻率估計(jì)方差理論均值的最小值及其與CRLB之比見表1。

    表1 頻率估計(jì)的方差最小值

    Tab.1 The minimum variance of the frequency estimates

    據(jù)表1結(jié)果可知,隨著逐漸增大,頻率估計(jì)方差最小值逐漸偏離CRLB,但差距十分微小。這表明,算法在≤0.5范圍內(nèi)性能十分穩(wěn)定且受值的影響極小,理論上,隨著→0頻率估計(jì)方差將完全逼近無偏頻率估計(jì)的CRLB。

    另一更為重要且關(guān)鍵的點(diǎn)在于:該結(jié)果表明算法具有可迭代執(zhí)行的能力,即通過頻移(不斷迭代),將≠0時(shí)的頻率估計(jì)轉(zhuǎn)化為→0時(shí)的頻率估計(jì),進(jìn)一步優(yōu)化提升算法性能,使得≠0時(shí)的頻率估計(jì)方差能一致逼近于CRLB,從而最大程度地降低噪聲對(duì)算法的干擾影響。迭代執(zhí)行程序的具體執(zhí)行步驟如下所示:

    (2)通過峰值搜索(κ)=argmax{|DFT()|}得到中心譜線,根據(jù)DTFT定義計(jì)算輔助譜線(κ±)。

    (4)=+1。

    (6)據(jù)DTFT定義計(jì)算(κ+1)及輔助譜線(κ+1)。

    (9)否則,重復(fù)步驟(4)~步驟(8)。

    值得注意的是,雖然理論上可取任意接近于零的數(shù)值,但當(dāng)極小時(shí)有限字長(zhǎng)效應(yīng)可能導(dǎo)致式(23)的分子和分母接近于零,使得頻率估計(jì)可能產(chǎn)生無意義的結(jié)果[23];且當(dāng)?shù)陀谀骋婚撝禃r(shí),受限于采樣點(diǎn)數(shù)和計(jì)算精度限制,頻率估計(jì)方差幾乎保持恒定而無法隨著減小而降低,具體結(jié)果見3.1節(jié)。

    3 仿真結(jié)果

    為全面驗(yàn)證算法的有效性和準(zhǔn)確性,分別進(jìn)行三個(gè)不同的仿真測(cè)試。第一個(gè)為分析所用譜線的間隔對(duì)算法在非迭代情況下性能的影響;第二個(gè)為驗(yàn)證算法在迭代情況下是否能夠一致逼近無偏頻率估計(jì)的CRLB;第三個(gè)為不同性能優(yōu)異的插值算法與本文算法的性能比對(duì)。

    3.1 譜線間隔i的影響分析

    噪聲干擾情況下,據(jù)附錄第2節(jié)可知,單一譜線及與其他譜線間的統(tǒng)計(jì)特性隨著譜線間隔的變化而變化,進(jìn)而導(dǎo)致不同值下頻率估計(jì)方差存在差異。設(shè)定信號(hào)采樣點(diǎn)數(shù)=1 024,幅度=2,采樣頻率s=1 024Hz,信號(hào)初始相位在[-p,p]范圍內(nèi)隨機(jī)變化,對(duì)應(yīng)信號(hào)波動(dòng)周期在[2.55, 3.45]以步長(zhǎng)0.1逐漸變化;噪聲干擾設(shè)置為加性高斯白噪聲,對(duì)應(yīng)信噪比為SNR=40dB,單一信號(hào)條件下重復(fù)= 10 000次獨(dú)立蒙特卡洛仿真測(cè)試,則頻率估計(jì)的MSE與CRLB之比隨的變化趨勢(shì)如圖2所示。

    圖2 頻率估計(jì)MSE

    圖2結(jié)果顯示:①頻率估計(jì)的MSE與理論方差結(jié)果一致吻合,表明噪聲擾動(dòng)下算法為嚴(yán)格無偏頻率估計(jì)。②當(dāng)≤0.5時(shí),頻率估計(jì)MSE在單一頻率條件下幾乎保持恒定,且在→0時(shí)十分接近于CRLB;反之,則隨著的增大而逐漸增大。其原因在于當(dāng)≤0.5時(shí),由于所用三條譜線都位于矩形窗主瓣范圍內(nèi),因而其信噪比相對(duì)較高;反之,隨著旁瓣衰減,輔助譜線的信噪比則很低,受噪聲擾動(dòng)影響較大。③當(dāng)≤0.5時(shí),頻率估計(jì)MSE不受頻率偏移符號(hào)影響,即在>0或<0范圍內(nèi)的MSE在量級(jí)上幾乎相等;當(dāng)>0.5時(shí)(特別是在→1時(shí)),在>0范圍內(nèi)的頻率估計(jì)MSE要小于<0范圍內(nèi),其原因在于當(dāng)<0時(shí),量級(jí)較小的輔助譜線受負(fù)頻率的頻譜泄漏影響將不可忽略。

    3.2 迭代有效性驗(yàn)證

    仍以3.1節(jié)中設(shè)定的信號(hào)為例,對(duì)應(yīng)信號(hào)波動(dòng)周期在[2.55, 3.45]則以步長(zhǎng)0.02逐漸變化。根據(jù)3.1節(jié)仿真結(jié)果,此處僅考慮≤0.5范圍內(nèi)取值時(shí)的頻率估計(jì)結(jié)果,即對(duì)比>0.5時(shí)保持相對(duì)最優(yōu);同時(shí),考慮實(shí)際應(yīng)用中數(shù)值穩(wěn)定性,此處僅示出=0.1時(shí)的頻率估計(jì)結(jié)果,如圖3所示。

    圖3 i=0.1時(shí)頻率估計(jì)的MSE

    圖3結(jié)果顯示,噪聲干擾下雖然所提無偏頻率估計(jì)算法在頻偏偏離0時(shí)無法一致逼近CRLB,但伴隨迭代程序的執(zhí)行,僅需1次迭代頻率估計(jì)的MSE在頻偏?[-0.5, 0.5]全域內(nèi)都一致趨近于最優(yōu)估計(jì)結(jié)果(→0時(shí)),即在任意頻偏情況下都十分接近于CRLB。在不顯著增加計(jì)算復(fù)雜度的前提下,迭代執(zhí)行程序有效降低了算法對(duì)頻率偏移的敏感度,極大地提高了算法的魯棒性。對(duì)于≤0.5范圍內(nèi)其他任意值,據(jù)圖2和第2節(jié)中算法統(tǒng)計(jì)特性分析可知,迭代程序仍有效且在噪聲干擾下頻偏估計(jì)結(jié)果將十分接近,不再贅述。

    3.3 性能比對(duì)

    為進(jìn)一步驗(yàn)證算法在不同波動(dòng)周期下的有效性和準(zhǔn)確性,設(shè)定信號(hào)波動(dòng)周期在[0.55, 4.45]范圍內(nèi)逐漸變化,其他信號(hào)參數(shù)保持與3.1節(jié)和3.2節(jié)一致。由于頻域類頻率估計(jì)算法計(jì)算復(fù)雜度通常要遠(yuǎn)遠(yuǎn)小于時(shí)域類參數(shù)化算法,分別為(log2)和(3),因此,僅選取目前為數(shù)不多的幾種考慮負(fù)頻率泄漏影響的頻域插值算法進(jìn)行性能比對(duì),包括基于頻譜系數(shù)加權(quán)的3IpDFT[18]、采用矩陣解析的3IpDFT[15]、基于補(bǔ)零和DFT系數(shù)變換的2IpDFT[11]、解析插值算法3IpDFT[10]以及改進(jìn)的迭代插值DTFT算法AM[20]。不同插值算法的頻率估計(jì)MSE結(jié)果如圖4所示。圖4a為SNR=40dB時(shí)?[0.55, 4.45]范圍內(nèi)的頻率估計(jì)結(jié)果;圖4b為本文算法單步迭代和AM多步迭代收斂的頻率估計(jì)結(jié)果比對(duì);圖4c為SNR?[-20, 100]dB范圍內(nèi)=1.25時(shí)的頻率估計(jì)結(jié)果。此外,圖4中以IpDTFT(,)表示本文算法,其中為譜線間隔,為迭代次數(shù);鑒于在≤0.5范圍內(nèi)算法性能十分接近,此處仍僅示出=0.1時(shí)的仿真結(jié)果。

    圖4 不同算法頻率估計(jì)性能比對(duì)

    圖4a結(jié)果顯示:①非迭代情況下,本文算法和2IpDFT整體上要遠(yuǎn)優(yōu)于其他參考比對(duì)算法。②當(dāng)信號(hào)波動(dòng)周期≤1時(shí),本文算法要優(yōu)于2IpDFT;當(dāng)信號(hào)波動(dòng)周期>1時(shí),在偏移||≤0.25范圍內(nèi)兩種算法性能十分接近;在0.25<||≤0.5范圍內(nèi)2IpDF要優(yōu)于本文算法。其原因在于,當(dāng)||≤0.25時(shí)兩種算法的解析方式極為接近,進(jìn)一步分析可知,當(dāng)且僅當(dāng)=0.5時(shí)兩者具有幾乎相同的解析插值公式,由3.1節(jié)結(jié)果可知,兩種算法性能相當(dāng);而當(dāng) 0.25<||≤0.5時(shí),由于2IpDFT算法通過補(bǔ)零改變了解析插值中所用譜線,使得其在該范圍內(nèi)的性能與||≤0.25時(shí)保持一致,因而優(yōu)于本文算法的非迭代結(jié)果。

    圖4b結(jié)果顯示:①當(dāng)信號(hào)波動(dòng)周期>1時(shí),通過多次不斷迭代收斂的AM算法與本文算法單次迭代結(jié)果幾乎保持一致,明顯受限于多次迭代計(jì)算,AM算法計(jì)算量將遠(yuǎn)遠(yuǎn)高于本文算法;②當(dāng)信號(hào)波動(dòng)周期≤1時(shí),受限于與輔助譜線間隔固定為=0.5,AM算法中所用譜線量級(jí)極小,受負(fù)頻率的頻譜泄漏影響遠(yuǎn)大于本文算法(=0.1),且因未采用信噪比較高的中心峰值譜線,導(dǎo)致算法MSE要遠(yuǎn)遠(yuǎn)高于本文算法。

    圖4c結(jié)果顯示:①在信號(hào)波動(dòng)周期=1.25時(shí),此時(shí)負(fù)頻率頻譜泄漏干擾相對(duì)占主導(dǎo),非迭代AM算法和3IpDFT算法存在明顯的估計(jì)偏差,因而性能上無法一致逼近CRLB;同理,迭代AM算法在信噪比SNR>80dB時(shí),因此處迭代次數(shù)設(shè)定限制無法完全消除估計(jì)偏差影響,導(dǎo)致其MSE偏離CRLB。②總體上本文算法(包括迭代)和其他算法對(duì)比,因準(zhǔn)確計(jì)及負(fù)頻率泄漏使得頻率估計(jì)結(jié)果幾乎一致逼近CRLB,在不同SNR情況下本文算法的單次迭代結(jié)果依然保持最優(yōu)。

    綜上可知,相較其他插值類算法,雖然本文算法在非迭代情況下未保持全域最優(yōu),但在不明顯增加計(jì)算量的前提下,僅通過單步迭代以優(yōu)化性能,即能達(dá)到最優(yōu)的頻率估計(jì)且全域一致逼近CRLB。

    表2 不同算法計(jì)算復(fù)雜度比對(duì)

    Tab.2 Comparison of the computational complexity of different algorithms

    4 實(shí)驗(yàn)結(jié)果

    實(shí)驗(yàn)驗(yàn)證測(cè)試平臺(tái)如圖5所示,以精度為0.05級(jí)的DK-56B1三相交直流指示儀表檢定裝置作為功率信號(hào)源,額定輸出信號(hào)頻率為50Hz,信號(hào)幅值為5V;以24bit的IOtech數(shù)據(jù)記錄儀作為信號(hào)采集裝置,采樣頻率為100kHz。以截取信號(hào)采樣點(diǎn)數(shù)的改變對(duì)應(yīng)于變化,對(duì)應(yīng)于每個(gè)通過單點(diǎn)移動(dòng)改變信號(hào)初始相位,并連續(xù)進(jìn)行10 000次分析。此外,鑒于白噪聲擾動(dòng)下IEEE 1057標(biāo)準(zhǔn)[24]建議的4參量正弦擬合算法(4PSF)提供了頻率的最大似然估計(jì),其估計(jì)方差非常接近CRLB,因此以該算法估計(jì)結(jié)果作為基準(zhǔn),其迭代次數(shù)設(shè)定為6。在?[1, 4]范圍內(nèi)不同算法的頻率估計(jì)MSE實(shí)驗(yàn)結(jié)果如圖6所示。

    圖5 實(shí)驗(yàn)平臺(tái)

    圖6 實(shí)驗(yàn)結(jié)果

    從圖6可看出,實(shí)驗(yàn)結(jié)果與圖4仿真結(jié)果基本完全吻合。此處,雖然實(shí)驗(yàn)中頻率分辨率隨著增大而增高,頻率估計(jì)的方差下界將逐漸降低,但本文算法的單步迭代結(jié)果在任意頻點(diǎn)上仍然與4PSF基準(zhǔn)結(jié)果基本保持一致,這也再次印證了本文算法的準(zhǔn)確性和迭代有效性。

    圖7 實(shí)測(cè)電力信號(hào)

    圖8 實(shí)測(cè)電壓信號(hào)的頻率估計(jì)

    如圖8a所示,當(dāng)采樣點(diǎn)數(shù)較小即≤175時(shí),實(shí)測(cè)電壓信號(hào)含有的一定量諧波泄漏干擾導(dǎo)致頻率估計(jì)波動(dòng)偏大;隨著采樣點(diǎn)數(shù)增多,諧波泄漏干擾逐漸降低直至忽略不計(jì),頻率估計(jì)趨于穩(wěn)定且能夠準(zhǔn)確反映負(fù)載運(yùn)行狀態(tài)變化。圖8b所示頻率估計(jì)MSE隨變化而變化的趨勢(shì)與圖6實(shí)驗(yàn)室測(cè)試結(jié)果基本一致,雖然在較小時(shí)受限于諧波泄漏干擾,導(dǎo)致與實(shí)驗(yàn)室測(cè)試結(jié)果出現(xiàn)一定偏差,但總體上本文算法的單步迭代頻率估計(jì)結(jié)果依然保持最優(yōu),而且極為接近于以4PSF算法為基準(zhǔn)近似的CRLB。

    5 結(jié)論

    針對(duì)短時(shí)波動(dòng)周期下實(shí)正弦信號(hào)頻率精確估計(jì)問題,首先利用頻譜實(shí)/虛部線性方程組解析得到了基于廣義DTFT插值的無偏頻率估計(jì)公式,不僅擺脫了對(duì)窗函數(shù)旁瓣衰減特性的嚴(yán)重依賴,而且打破了以往傳統(tǒng)頻域插值算法僅限于單頻解析信號(hào)和離散DFT譜線插值的局限性;進(jìn)而基于噪聲背景下頻率估計(jì)的統(tǒng)計(jì)性能分析導(dǎo)出了算法的迭代執(zhí)行方式,解決了算法在背景噪聲干擾下無法保持全域最優(yōu)且一致逼近CRLB的問題。

    仿真和實(shí)驗(yàn)結(jié)果顯示,一方面,得益于算法解析過程準(zhǔn)確計(jì)及來自負(fù)頻率的頻譜泄漏干擾,且采用間隔低至0.1D的DTFT輔助譜線,即便信號(hào)波動(dòng)周期≤1時(shí)算法在未迭代情況下性能仍保持穩(wěn)定;另一方面,對(duì)于噪聲干擾影響,算法通過簡(jiǎn)單有效的迭代執(zhí)行方式以優(yōu)化估計(jì)性能,在不顯著增加計(jì)算量和復(fù)雜度的前提下,僅需單步迭代即能達(dá)到最優(yōu)的頻率估計(jì)且全域一致逼近CRLB,在=0.1時(shí)理論上算法最小MSE與CRLB之比接近1.008 822。

    對(duì)于諧波擾動(dòng)的多頻信號(hào)頻率估計(jì),在主瓣頻譜不存在混疊的前提下,可首先采用本文算法得到單頻點(diǎn)的頻率粗估計(jì),進(jìn)一步通過單頻點(diǎn)的粗估計(jì)實(shí)現(xiàn)頻譜重構(gòu),逐步迭代消除諧波頻譜長(zhǎng)泄漏的干擾影響,從而再次采用本文算法實(shí)現(xiàn)精確頻率估計(jì)。

    這種針對(duì)多頻信號(hào)頻率估計(jì)的算法具體實(shí)現(xiàn)及相應(yīng)測(cè)試結(jié)果將在下一步研究工作中詳細(xì)探討。

    附 錄

    1.的具體求解過程

    便于推導(dǎo),首先將矩陣+/-和+/-中與位置索引相關(guān)的已知三角函數(shù)變量分別以ab表示為

    由此,式(13)和式(14)改寫為

    則根據(jù)行列式性質(zhì)將式(A2)左側(cè)行列式分別進(jìn)行簡(jiǎn)化,即

    再次利用行列式性質(zhì),將式(A6)虛數(shù)部分乘以虛數(shù)單位j后與式(A5)實(shí)數(shù)部分進(jìn)行合并,進(jìn)而經(jīng)過展開計(jì)算后可得的解析解為

    2. 頻率估計(jì)的統(tǒng)計(jì)方差推導(dǎo)

    考慮中噪聲項(xiàng)之間的相關(guān)性,首先將D展 開,即

    令=-,=-和=-,由于[Re(D)]=0,Re(D)的自相關(guān)與方差相等,即

    式中,由于各參量方差及互相關(guān)具有相似性,因此,此處僅分別以的方差以及和的互相關(guān)為例進(jìn)行推導(dǎo),即的方差為

    其中

    (A13)

    其中

    由此,據(jù)式(A10)~式(A18)可得Re(D)的方差為

    式中,P和N分別為

    [1] Sun J, Aboutanios E, Smith D B, et al. Robust frequency, phase, and amplitude estimation in power Systems considering harmonics[J]. IEEE Transactions on Power Delivery, 2020, 35(3): 1158-1168.

    [2] 張政, 溫和, 黎福海, 等. 多水平集單周期電力系統(tǒng)頻率測(cè)量方法及應(yīng)用[J]. 電工技術(shù)學(xué)報(bào), 2017, 32(7): 119-127.

    Zhang Zheng, Wen He, Li Fuhai, et al. Multi-level set single-cycle power system frequency measurement method and its application[J]. Transactions of China Electrotechnical Society, 2017, 32(7): 119-127.

    [3] 肖勇, 趙偉, 黃松嶺. 基于離散傅里葉級(jí)數(shù)的非同步采樣下諧波功率測(cè)量算法[J]. 電工技術(shù)學(xué)報(bào), 2018, 33(7): 1570-1578.

    Xiao Yong, Zhao Wei, Huang Songling. Harmonic power measurement algorithm based on discrete fourier series in asynchronous sampling[J]. Transa- ctions of China Electrotechnical Society, 2018, 33(7): 1570-1578.

    [4] Jafarpisheh B, Madani S M, Parvaresh F, et al. Power system frequency estimation using adaptive accelerated MUSIC[J]. IEEE Transactions on Instrumentation and Measurement, 2018, 67(11): 2592-2602.

    [5] 楊超, 張淮清, 王耀, 等. 計(jì)及全泄漏影響的多點(diǎn)插值離散傅里葉變換校正方法[J]. 電工技術(shù)學(xué)報(bào), 2020, 35(16): 3385-3395.

    Yang Chao, Zhang Huaiqing, Wang Yao, et al. Multipoint interpolated discrete fourier transform correction method considering total leakage effect[J]. Transactions of China Electrotechnical Society, 2020, 35(16): 3385-3395.

    [6] Shin D, Kwak C, Kim G. An efficient algorithm for frequency estimation from cosine-sum windowed DFT coefficients[J]. Signal Processing, 2020, 166: 107245.1-107245.9.

    [7] 孫仲民, 何正友, 臧天磊. 一種混合卷積窗及其在諧波分析中的應(yīng)用[J]. 電工技術(shù)學(xué)報(bào), 2016, 31(16): 207-214.

    Sun Zhongmin, He Zhengyou, Zang Tianlei. A kind of hybrid convolution window and its application in harmonic analysis[J]. Transactions of China Electro- technical Society, 2016, 31(16): 207-214.

    [8] Romano P, Paolone M. Enhanced interpolated-DFT for synchrophasor estimation in FPGAs: theory, implementation, and validation of a PMU prototype[J]. IEEE Transactions on Instrumentation and Measure- ment, 2014, 63(12): 2824-2836.

    [9] 翟曉軍, 周波. 一種改進(jìn)的插值FFT諧波分析算法[J]. 中國電機(jī)工程學(xué)報(bào), 2016, 36(11): 2952-2958.

    Zhai Xiaojun, Zhou Bo. An improved interpolated FFT algorithm for harmonic analysis[J]. Proceedings of the CSEE, 2016, 36(11): 2952-2958.

    [10] Wang Kai, Wen He, Li Guoqing. Accurate frequency estimation by using three-point interpolated discrete fourier transform based on rectangular window[J]. IEEE Transactions on Industrial Informatics, 2021, 17(1): 73-81.

    [11] Sun Yuxin, Zhuang Chungang, Xiong Zhenhua. A scale factor based interpolated DFT for chatter frequency estimation[J]. IEEE Transactions on Instrumentation and Measurement, 2015, 64(10): 2666-2678.

    [12] Sun Yuxin, Zhuang Chungang, Xiong Zhenhua. A switch-based inter-polated DFT for the small number of acquired sine wave cycles[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(4): 846-855.

    [13] Hlupic N, Basch D, Fertalj K. A note on optimality of analytical leakage compensation at boundary fre- quencies[J]. IEEE Transactions on Instrumentation and Measurement, 2010, 59(2): 301-308.

    [14] Borkowski J, Kania D, Mroczka J. Interpolated-DFT- based fast and accurate frequency estimation for the control of power[J]. IEEE Transactions on Industrial Electronics, 2014, 61(12): 7026-7034.

    [15] Chen Kuifu, Cao Xu, Li Yangfeng. Sine wave fitting to short records initialized with the frequency retrieved from Hanning windowed FFT spectrum[J]. Measurement, 2009, 42(1): 127-135.

    [16] Murakami T, Wang W. An analytical solution to Jacobsen estimator for windowed signals[C]//IEEE International Conference on Acoustics Speech and Signal Processing (ICASSP), Barcelona, Spain, 2020: 5950-5954.

    [17] Belega D, Petri D. Effect of noise and harmonics on sine-wave frequency estimation by interpolated DFT algorithms based on few observed cycles[J]. Signal Processing, 2017, 140: 207-218.

    [18] Belega D, Petri D, Dallet D. Frequency estimation of a sinusoidal signal via a three-point interpolated DFT method with high image component interference rejection capability[J]. Digital Signal Processing, 2014, 24: 162-169.

    [19] Candan C. A method for fine resolution frequency estimation from three DFT samples[J]. IEEE Signal Process, Lett, 2011, 18(6): 351-354.

    [20] Ye S, Sun J, Aboutanios E. On the estimation of the parameters of a real sinusoid in noise[J]. IEEE Signal Processing Letters, 2017, 24(5): 638-642.

    [21] Aboutanios E, Mulgrew B. Iterative frequency estimation by interpolation on Fourier coefficients[J]. IEEE Transactions on Signal Processing, 2005, 53(4): 1237-1242.

    [22] Horn R A, Johnson C R. Matrix analysis[M]. Cambridge: Cambridge University Pre, 1985.

    [23] Lei Fan, Guoqing Qi. Frequency estimator of sinusoid based on interpolation of three DFT spectral lines[J]. Signal Processing, 2018, 144: 52-60.

    [24] IEEE standard for digitizing waveform recorders[S]. New Jersey: Piscataway, 2007.

    Accurate Frequency Estimation Based on Analytical Interpolated Discrete Time Fourier Transform

    11,21,31

    (1. State Key Laboratory of Power Transmission Equipment & System Security and New Technology Chongqing University Chongqing 400044 China 2. Electric Power Research Institute of Yunnan Power Grid Co. Ltd Kunming 650217 China 3. Chongqing Triloop Prospecting Technology Co. Ltd Chongqing 402660 China)

    When the cycles in record (CiR) of the signal is small, the spectral leakage from the image component is the main source of estimation errors in the frequency-domain interpolation. To eliminate the influence of spectral leakage and accurately estimate the frequency of a real sinusoidal signal with any CiR, an analytic interpolation algorithm based on generalized DTFT is proposed. The algorithm uses any DTFT spectral lines located in the main lobe of the spectrum to construct a linear equation set corresponding to the real and imaginary parts, and then an unbiased frequency estimator based on analytic interpolated DTFT is obtained. Besides, considering the influence of additive white noise in practice, a simple and effective iterative method is presented based on the analysis of the statistical characteristics of the frequency estimator. The iterative method ensures that the mean square errors (MSE) of the frequency estimation can consistently achieve the Cramer-Rao lower bound (CRLB) of the unbiased frequency estimation under any CiR. Simulation and experimental results confirm that the proposed algorithm is superior to other existing interpolated DFT algorithms.

    Small cycles in record, spectral leakage, analytically interpolated discrete time fourier transform (DTFT), iterative method, Cramer-Rao lower bound (CRLB)

    10.19595/j.cnki.1000-6753.tces.201667

    TM935

    楊 超 男,1986年生,博士,研究方向?yàn)閿?shù)字信號(hào)處理及其在精密測(cè)量、電能計(jì)量與電能質(zhì)量分析中的深度應(yīng)用。E-mail: 328383369@qq.com(通信作者)

    李 波 男,1982年生,碩士,高級(jí)工程師,主要從事電氣測(cè)量以及智能電網(wǎng)測(cè)控研究工作。E-mail: libo2010@163.com

    2020-12-20

    2021-04-27

    國家重點(diǎn)研發(fā)計(jì)劃課題(2018YFC0406904)和重慶市重點(diǎn)產(chǎn)業(yè)共性關(guān)鍵技術(shù)創(chuàng)新專項(xiàng)(重點(diǎn)研發(fā)項(xiàng)目)(CSTC2017ZDCY-ZDYFX0045)資助。

    (編輯 陳 誠)

    猜你喜歡
    信號(hào)
    信號(hào)
    鴨綠江(2021年35期)2021-04-19 12:24:18
    完形填空二則
    7個(gè)信號(hào),警惕寶寶要感冒
    媽媽寶寶(2019年10期)2019-10-26 02:45:34
    孩子停止長(zhǎng)個(gè)的信號(hào)
    《鐵道通信信號(hào)》訂閱單
    基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
    電子制作(2018年11期)2018-08-04 03:25:42
    基于Arduino的聯(lián)鎖信號(hào)控制接口研究
    《鐵道通信信號(hào)》訂閱單
    基于LabVIEW的力加載信號(hào)采集與PID控制
    Kisspeptin/GPR54信號(hào)通路促使性早熟形成的作用觀察
    久久久久久久久大av| 在线观看人妻少妇| 女人久久www免费人成看片| 日韩av在线免费看完整版不卡| 最近最新中文字幕免费大全7| 少妇猛男粗大的猛烈进出视频| 久久99热6这里只有精品| 99久久精品一区二区三区| 最近中文字幕2019免费版| 男女下面进入的视频免费午夜| 国产高清国产精品国产三级 | 欧美激情国产日韩精品一区| 日本av免费视频播放| 精品一区在线观看国产| 青春草国产在线视频| 99久久综合免费| av女优亚洲男人天堂| 直男gayav资源| 精品酒店卫生间| av国产免费在线观看| 青青草视频在线视频观看| 少妇人妻精品综合一区二区| 久久亚洲国产成人精品v| 日韩人妻高清精品专区| 精品一区在线观看国产| 啦啦啦啦在线视频资源| 人妻制服诱惑在线中文字幕| 一区二区av电影网| 一级毛片我不卡| 国产爽快片一区二区三区| 亚洲最大成人中文| av国产免费在线观看| 黄色怎么调成土黄色| 国产视频内射| 国产在线视频一区二区| 美女主播在线视频| 久久午夜福利片| 国产久久久一区二区三区| 亚洲中文av在线| 五月伊人婷婷丁香| 亚洲国产最新在线播放| 国产乱人偷精品视频| 亚洲av二区三区四区| 精品人妻一区二区三区麻豆| 成人亚洲精品一区在线观看 | 久久久久性生活片| 欧美+日韩+精品| 在线观看免费日韩欧美大片 | 亚洲内射少妇av| 麻豆成人午夜福利视频| 在线播放无遮挡| freevideosex欧美| 亚洲欧美日韩卡通动漫| 国产乱人偷精品视频| 黑丝袜美女国产一区| 国产伦精品一区二区三区视频9| 啦啦啦视频在线资源免费观看| 国产精品蜜桃在线观看| 久久久久久久久久久免费av| 国产精品国产三级国产av玫瑰| 性色avwww在线观看| 亚洲国产精品国产精品| av福利片在线观看| 久久人人爽人人爽人人片va| av女优亚洲男人天堂| 精品一品国产午夜福利视频| 自拍偷自拍亚洲精品老妇| 一级二级三级毛片免费看| 天堂俺去俺来也www色官网| 日韩成人伦理影院| 少妇丰满av| 午夜福利视频精品| 国产精品麻豆人妻色哟哟久久| 国产精品久久久久久久电影| 亚洲欧洲日产国产| 久久久久国产精品人妻一区二区| 熟女av电影| 蜜臀久久99精品久久宅男| 久久影院123| 国产视频内射| 欧美人与善性xxx| 亚洲天堂av无毛| 成人漫画全彩无遮挡| 欧美丝袜亚洲另类| 少妇高潮的动态图| 爱豆传媒免费全集在线观看| 免费黄色在线免费观看| 成人高潮视频无遮挡免费网站| 日韩中字成人| 99久国产av精品国产电影| 国产免费又黄又爽又色| 亚洲av国产av综合av卡| 国产av国产精品国产| 纯流量卡能插随身wifi吗| 看十八女毛片水多多多| 又粗又硬又长又爽又黄的视频| 男人爽女人下面视频在线观看| 在线观看免费视频网站a站| 亚洲av国产av综合av卡| 毛片一级片免费看久久久久| 乱码一卡2卡4卡精品| 亚洲国产高清在线一区二区三| av又黄又爽大尺度在线免费看| 国产亚洲精品久久久com| 成人美女网站在线观看视频| 亚洲欧美精品专区久久| 插阴视频在线观看视频| 大又大粗又爽又黄少妇毛片口| 22中文网久久字幕| 成人18禁高潮啪啪吃奶动态图 | 一级毛片久久久久久久久女| 国产乱人偷精品视频| 亚洲欧美一区二区三区黑人 | 国产中年淑女户外野战色| 日韩成人伦理影院| 国产成人91sexporn| 日本一二三区视频观看| 午夜视频国产福利| 久久精品国产自在天天线| 亚洲婷婷狠狠爱综合网| 亚洲国产毛片av蜜桃av| 蜜臀久久99精品久久宅男| 最近中文字幕高清免费大全6| 欧美成人午夜免费资源| 日韩av不卡免费在线播放| .国产精品久久| 韩国av在线不卡| 欧美zozozo另类| 久久久久视频综合| 午夜精品国产一区二区电影| 一区二区三区四区激情视频| 成人无遮挡网站| 毛片一级片免费看久久久久| 婷婷色综合www| 青春草国产在线视频| 久热久热在线精品观看| 久久精品国产鲁丝片午夜精品| 久久久久久人妻| 黄色欧美视频在线观看| 午夜免费男女啪啪视频观看| 日韩,欧美,国产一区二区三区| 黑丝袜美女国产一区| 国产av精品麻豆| 久久久亚洲精品成人影院| 国产极品天堂在线| 精品人妻视频免费看| 亚洲欧美成人综合另类久久久| 久久久久久伊人网av| 国产精品人妻久久久影院| 伦理电影免费视频| 成人漫画全彩无遮挡| 国产精品爽爽va在线观看网站| 欧美一区二区亚洲| 国产精品久久久久久av不卡| 王馨瑶露胸无遮挡在线观看| 国产av国产精品国产| 日本-黄色视频高清免费观看| a 毛片基地| 亚洲国产精品成人久久小说| 亚洲欧美清纯卡通| 最后的刺客免费高清国语| 国产成人免费无遮挡视频| 久久久久久久精品精品| 欧美zozozo另类| 大香蕉97超碰在线| 日本欧美国产在线视频| 久久精品久久精品一区二区三区| 日本vs欧美在线观看视频 | 午夜福利在线观看免费完整高清在| 夫妻性生交免费视频一级片| 黑丝袜美女国产一区| 精品国产露脸久久av麻豆| 男的添女的下面高潮视频| 久久热精品热| 国产真实伦视频高清在线观看| 99久久精品热视频| 有码 亚洲区| 日本av手机在线免费观看| 久久99热这里只频精品6学生| videossex国产| 观看美女的网站| 美女内射精品一级片tv| 乱系列少妇在线播放| 91精品国产九色| 精品亚洲乱码少妇综合久久| 在线观看免费视频网站a站| 国产又色又爽无遮挡免| 日韩,欧美,国产一区二区三区| 性高湖久久久久久久久免费观看| www.av在线官网国产| 精品国产露脸久久av麻豆| 男的添女的下面高潮视频| 一区在线观看完整版| 十分钟在线观看高清视频www | 天堂中文最新版在线下载| 亚洲精品视频女| 五月伊人婷婷丁香| 国产欧美另类精品又又久久亚洲欧美| 免费久久久久久久精品成人欧美视频 | 少妇猛男粗大的猛烈进出视频| 大香蕉97超碰在线| 麻豆乱淫一区二区| 蜜桃久久精品国产亚洲av| 青春草视频在线免费观看| www.av在线官网国产| 22中文网久久字幕| 国产精品秋霞免费鲁丝片| 五月玫瑰六月丁香| 久久久久人妻精品一区果冻| 在线观看免费高清a一片| av不卡在线播放| 麻豆精品久久久久久蜜桃| 亚洲av不卡在线观看| 精品久久久噜噜| 国产毛片在线视频| 黑丝袜美女国产一区| 午夜日本视频在线| 国产成人91sexporn| 麻豆国产97在线/欧美| 亚洲高清免费不卡视频| 美女高潮的动态| 国产色爽女视频免费观看| 午夜免费男女啪啪视频观看| 国产男女超爽视频在线观看| 久久热精品热| 菩萨蛮人人尽说江南好唐韦庄| 大码成人一级视频| 国产av国产精品国产| 日韩 亚洲 欧美在线| 在线观看人妻少妇| 久久精品熟女亚洲av麻豆精品| 九草在线视频观看| 91午夜精品亚洲一区二区三区| 日韩不卡一区二区三区视频在线| 午夜激情久久久久久久| av国产精品久久久久影院| 亚洲av免费高清在线观看| 香蕉精品网在线| 在线观看一区二区三区激情| 国产精品99久久久久久久久| 一级毛片黄色毛片免费观看视频| 国产极品天堂在线| 亚洲国产av新网站| 久久影院123| 麻豆成人午夜福利视频| 国产在线一区二区三区精| 极品少妇高潮喷水抽搐| 高清视频免费观看一区二区| 你懂的网址亚洲精品在线观看| 2022亚洲国产成人精品| 亚洲怡红院男人天堂| 狂野欧美激情性xxxx在线观看| 日日摸夜夜添夜夜爱| 在现免费观看毛片| 狂野欧美激情性bbbbbb| 蜜桃亚洲精品一区二区三区| 精品一品国产午夜福利视频| 国产精品99久久久久久久久| 两个人的视频大全免费| 精品国产乱码久久久久久小说| 日韩中文字幕视频在线看片 | 久久婷婷青草| 看十八女毛片水多多多| 亚洲精品自拍成人| 亚洲av.av天堂| 亚洲国产精品一区三区| 99久国产av精品国产电影| 亚洲av综合色区一区| 久久精品熟女亚洲av麻豆精品| 舔av片在线| 性色avwww在线观看| 极品教师在线视频| 亚洲人成网站在线观看播放| 国产精品久久久久久久电影| 婷婷色综合大香蕉| 亚洲av电影在线观看一区二区三区| 国模一区二区三区四区视频| 国产精品国产三级国产专区5o| h视频一区二区三区| 亚洲经典国产精华液单| 蜜臀久久99精品久久宅男| 日韩大片免费观看网站| 女人十人毛片免费观看3o分钟| 狂野欧美白嫩少妇大欣赏| 国产日韩欧美在线精品| 在线看a的网站| 日韩制服骚丝袜av| 久久国内精品自在自线图片| 青春草视频在线免费观看| 男人舔奶头视频| 久久久久久久精品精品| 亚洲人成网站在线播| 亚洲经典国产精华液单| 色吧在线观看| 在线免费十八禁| 人人妻人人看人人澡| 只有这里有精品99| 在线观看免费高清a一片| 在线观看人妻少妇| 国内少妇人妻偷人精品xxx网站| 日韩成人av中文字幕在线观看| 久久久久久久久久久丰满| 亚洲电影在线观看av| 国产精品欧美亚洲77777| 国产视频内射| 久久久久国产网址| 国产伦精品一区二区三区视频9| 插阴视频在线观看视频| 国产成人精品一,二区| 男的添女的下面高潮视频| 国产高清国产精品国产三级 | 亚洲一区二区三区欧美精品| 亚洲国产精品一区三区| 午夜精品国产一区二区电影| av天堂中文字幕网| 国产亚洲欧美精品永久| 青春草国产在线视频| 免费大片黄手机在线观看| 久久99热这里只有精品18| 国产精品熟女久久久久浪| 国产精品国产av在线观看| 观看av在线不卡| 插阴视频在线观看视频| 在线观看一区二区三区激情| 亚洲一区二区三区欧美精品| 国产av国产精品国产| 最近最新中文字幕大全电影3| 婷婷色麻豆天堂久久| 久久97久久精品| 中文字幕人妻熟人妻熟丝袜美| 日本黄大片高清| 青春草亚洲视频在线观看| 亚洲精品国产av蜜桃| 亚洲欧美精品自产自拍| 亚洲精品国产成人久久av| 内地一区二区视频在线| 99视频精品全部免费 在线| 国产精品三级大全| 18禁动态无遮挡网站| 国产熟女欧美一区二区| 夫妻午夜视频| 如何舔出高潮| 亚洲国产高清在线一区二区三| 免费av中文字幕在线| 欧美成人精品欧美一级黄| av天堂中文字幕网| 色哟哟·www| 中文字幕制服av| 欧美区成人在线视频| 一级片'在线观看视频| 热re99久久精品国产66热6| 久久久久久久久久久免费av| 97超视频在线观看视频| 亚洲欧美中文字幕日韩二区| 99久久人妻综合| 美女脱内裤让男人舔精品视频| 亚洲国产av新网站| 99久久综合免费| 一级毛片aaaaaa免费看小| 精品久久久久久久久亚洲| 久久国内精品自在自线图片| 成人午夜精彩视频在线观看| 黄片wwwwww| 精品人妻熟女av久视频| 亚洲最大成人中文| 亚洲自偷自拍三级| 欧美丝袜亚洲另类| 精品酒店卫生间| 亚洲内射少妇av| 国产精品蜜桃在线观看| 大话2 男鬼变身卡| 天美传媒精品一区二区| 乱码一卡2卡4卡精品| 欧美亚洲 丝袜 人妻 在线| 大片免费播放器 马上看| 久久这里有精品视频免费| 亚洲熟女精品中文字幕| 2021少妇久久久久久久久久久| 国产成人一区二区在线| 久久99热这里只有精品18| 熟女人妻精品中文字幕| 亚洲色图av天堂| 一边亲一边摸免费视频| 老熟女久久久| 久久热精品热| 观看美女的网站| av在线老鸭窝| av播播在线观看一区| 大片免费播放器 马上看| 精品视频人人做人人爽| 黄片wwwwww| 麻豆成人午夜福利视频| 国产在线免费精品| 色哟哟·www| 亚洲精品乱码久久久v下载方式| 亚洲精品乱码久久久久久按摩| 午夜老司机福利剧场| 色5月婷婷丁香| 日韩中文字幕视频在线看片 | 国产精品成人在线| 欧美成人一区二区免费高清观看| 少妇的逼水好多| 十分钟在线观看高清视频www | 秋霞在线观看毛片| 高清黄色对白视频在线免费看 | 卡戴珊不雅视频在线播放| 成人一区二区视频在线观看| 国内揄拍国产精品人妻在线| 午夜激情福利司机影院| 成人亚洲欧美一区二区av| 亚洲精品国产成人久久av| 亚洲精品乱码久久久v下载方式| 亚洲精品久久午夜乱码| 欧美3d第一页| 午夜福利网站1000一区二区三区| 久久久久久久精品精品| 成年女人在线观看亚洲视频| 久久99热这里只有精品18| 妹子高潮喷水视频| 我的老师免费观看完整版| 91精品国产国语对白视频| 亚洲精华国产精华液的使用体验| 久久99热这里只有精品18| 韩国av在线不卡| 国产免费福利视频在线观看| 久久6这里有精品| 国产午夜精品一二区理论片| 成人国产麻豆网| 老司机影院成人| 黑丝袜美女国产一区| 一级a做视频免费观看| 91精品国产九色| 亚洲精品视频女| 精品久久久精品久久久| 亚洲在久久综合| 午夜激情久久久久久久| 直男gayav资源| 久久国产乱子免费精品| 国产精品一区www在线观看| 网址你懂的国产日韩在线| 我要看日韩黄色一级片| 国产精品麻豆人妻色哟哟久久| 国产一区二区三区av在线| 国产美女午夜福利| 国产av精品麻豆| 欧美zozozo另类| 色视频在线一区二区三区| 国产淫片久久久久久久久| av视频免费观看在线观看| 视频区图区小说| 美女高潮的动态| 亚洲精品自拍成人| 日韩亚洲欧美综合| 夫妻性生交免费视频一级片| 午夜免费鲁丝| 80岁老熟妇乱子伦牲交| 亚洲一级一片aⅴ在线观看| 男人狂女人下面高潮的视频| 国产女主播在线喷水免费视频网站| 王馨瑶露胸无遮挡在线观看| 看免费成人av毛片| 欧美一区二区亚洲| 久久这里有精品视频免费| 亚洲精品乱久久久久久| 国产v大片淫在线免费观看| 22中文网久久字幕| 日韩亚洲欧美综合| 欧美区成人在线视频| 国产成人免费无遮挡视频| 欧美日本视频| 精品久久国产蜜桃| 国产成人a区在线观看| 一区二区av电影网| 三级国产精品片| 男女啪啪激烈高潮av片| 五月天丁香电影| 日本免费在线观看一区| 欧美xxxx黑人xx丫x性爽| 国产在线免费精品| 一级毛片久久久久久久久女| 成人一区二区视频在线观看| 免费不卡的大黄色大毛片视频在线观看| 高清毛片免费看| 黄片wwwwww| 欧美3d第一页| 在线观看免费日韩欧美大片 | 欧美xxⅹ黑人| 一二三四中文在线观看免费高清| 18+在线观看网站| 亚洲国产欧美人成| 国产成人a区在线观看| 我的女老师完整版在线观看| 久久精品国产a三级三级三级| 91aial.com中文字幕在线观看| 久久久久网色| 最近中文字幕2019免费版| 美女内射精品一级片tv| 看免费成人av毛片| 免费观看性生交大片5| av视频免费观看在线观看| 777米奇影视久久| 老女人水多毛片| 亚洲美女视频黄频| 嘟嘟电影网在线观看| freevideosex欧美| 久久久久国产精品人妻一区二区| 1000部很黄的大片| 日韩亚洲欧美综合| 我的老师免费观看完整版| 国产亚洲精品久久久com| 国产高清不卡午夜福利| 欧美zozozo另类| 黄色配什么色好看| 国产av国产精品国产| 女人久久www免费人成看片| 国产成人a∨麻豆精品| 女人久久www免费人成看片| 91精品国产国语对白视频| 国产在线免费精品| 高清不卡的av网站| 日韩,欧美,国产一区二区三区| 国产高清有码在线观看视频| 欧美激情国产日韩精品一区| 在线观看国产h片| 97超碰精品成人国产| 精品一区在线观看国产| 国产精品国产av在线观看| 欧美另类一区| 少妇的逼水好多| 精品一区在线观看国产| 国产精品国产av在线观看| 成人特级av手机在线观看| 人人妻人人看人人澡| 国产精品欧美亚洲77777| 欧美日韩亚洲高清精品| 只有这里有精品99| 日本黄色片子视频| 3wmmmm亚洲av在线观看| 亚洲国产成人一精品久久久| 亚洲熟女精品中文字幕| 国产精品久久久久久久久免| 日本色播在线视频| 国产v大片淫在线免费观看| 最黄视频免费看| 内地一区二区视频在线| 日韩制服骚丝袜av| 日日啪夜夜爽| 最近最新中文字幕免费大全7| tube8黄色片| 高清av免费在线| 在线观看免费视频网站a站| 国产成人精品婷婷| 欧美97在线视频| 一级a做视频免费观看| 欧美变态另类bdsm刘玥| 国产高清国产精品国产三级 | 国国产精品蜜臀av免费| 小蜜桃在线观看免费完整版高清| 久久久久人妻精品一区果冻| 精品少妇黑人巨大在线播放| 亚洲aⅴ乱码一区二区在线播放| 精品酒店卫生间| 丰满少妇做爰视频| 在线看a的网站| 欧美xxxx性猛交bbbb| 这个男人来自地球电影免费观看 | 亚洲第一区二区三区不卡| 亚洲国产精品成人久久小说| 日韩av在线免费看完整版不卡| 国产视频内射| av免费观看日本| 亚洲图色成人| 国产欧美另类精品又又久久亚洲欧美| 大又大粗又爽又黄少妇毛片口| 国产国拍精品亚洲av在线观看| 久久毛片免费看一区二区三区| 精品一区在线观看国产| 91狼人影院| 啦啦啦视频在线资源免费观看| 欧美日韩综合久久久久久| 老司机影院成人| av在线观看视频网站免费| 如何舔出高潮| 少妇人妻久久综合中文| tube8黄色片| 欧美另类一区| 精品少妇黑人巨大在线播放| 国产免费视频播放在线视频| 熟女电影av网| 在现免费观看毛片| 日韩av在线免费看完整版不卡| 免费看av在线观看网站| 日本wwww免费看| 黄色欧美视频在线观看| 亚洲精品国产成人久久av| 18+在线观看网站| 精品少妇久久久久久888优播| 黄色日韩在线| 一本一本综合久久| 日日啪夜夜撸| 欧美日韩亚洲高清精品| 日韩大片免费观看网站| 日韩欧美一区视频在线观看 | 99久久综合免费| 国产精品久久久久久久电影| 狂野欧美激情性bbbbbb| 哪个播放器可以免费观看大片| 亚洲国产精品国产精品| 一级片'在线观看视频| 欧美老熟妇乱子伦牲交| 六月丁香七月| 成人毛片a级毛片在线播放| 免费看av在线观看网站| 少妇人妻久久综合中文| 成人综合一区亚洲|