１／ｆ噪聲的仿真產(chǎn)生和驗(yàn)證

摘 要：本文介紹了1/f噪聲的含義和原理，依照這些原理，通過白噪聲加簡(jiǎn)化得到的有理1/f數(shù)字濾波器的方法，使用Matlab編程，仿真產(chǎn)生了真實(shí)1/f噪聲。并根據(jù)常見的計(jì)算功率譜密度的算法，驗(yàn)證了仿真產(chǎn)生的1/f噪聲的性能。

關(guān)鍵詞:噪聲；1/f噪聲；仿真；功率譜密度；

中圖法分類號(hào):TB533 文獻(xiàn)標(biāo)志碼:B

Simulation and Verification of 1/f Noise

WU Yuan-zhen，WU MU-zhi

Suzhou Testing Instrument Factory Suzhou Jiangsu 215129 China

ABSTRACT:In this paper， the principles of the 1/f noise are discussed. In accordance with these principles， using the white noise to be added to simplify the right 1/f digital filter method; the author wrote programs based on Matlab and gets a better simulation of 1/f noise. And through the common computing power spectral density way of the algorithm， the verification is generated by the 1/f noise performance， the 1/f noise characteristics are shown.

Key words:Matlab;noise;1/f noise;simulation;PSD

近年來，隨著在物理、化學(xué)、生物醫(yī)學(xué)等基礎(chǔ)學(xué)科中對(duì)自然現(xiàn)象及規(guī)律的深入探索，在通信、測(cè)量、導(dǎo)航、醫(yī)療儀器等工業(yè)領(lǐng)域中人們必須測(cè)量及接受越來越小的電信號(hào)。測(cè)量或者接受裝置中的內(nèi)部噪聲影響越來越為明顯，這種噪聲已經(jīng)成為制約電子儀器和裝備進(jìn)一步提高信號(hào)檢測(cè)靈敏度及改進(jìn)接受信號(hào)質(zhì)量的關(guān)鍵因素。仿真產(chǎn)生噪聲信號(hào)可以幫助研究人員研究克服噪聲的方法，幫助技術(shù)人員研制出高抗干擾的儀器設(shè)備或指出改進(jìn)產(chǎn)品的方向。

1/f噪聲廣泛存在于電子元器件和電子線路中，從1/f噪聲測(cè)量分析中可以得到很多有用的信息。本文主要討論1/f噪聲的特點(diǎn)和仿真產(chǎn)生的方法。

Matlab的信號(hào)處理工具箱是信號(hào)算法文件的集合，它處理的對(duì)象是信號(hào)和系統(tǒng)，信號(hào)處理工具箱提供了很多命令行函數(shù)來進(jìn)行譜分析，包括經(jīng)典的無界技術(shù)以及特征值技術(shù)。另外，也添加了很多對(duì)象，用以提高可用性。因此我們使用該工具箱來仿真產(chǎn)生1/f噪聲源，并且用其來分析驗(yàn)證仿真產(chǎn)生出的噪聲可信度。

1 1/f噪聲的特性

1.1 1/f噪聲（低頻率噪聲）的特性

電子器件中的1/f 噪聲具有兩個(gè)基本特性：

1)1/f 噪聲在一個(gè)相當(dāng)寬的頻帶內(nèi)，功率譜密度與頻率f 成反比， 且頻帶上、下限都是有限值，上限頻率視1/f 噪聲與背景白噪聲的相對(duì)大小而定；下限頻率接近直流時(shí)功率密度仍呈很好的1/f 特性。

2)1/f 噪聲電壓和電流的功率譜密度近似與流經(jīng)器件的電流的平方成正比，這意味著1/f噪聲起源于電阻的漲落。設(shè)流過電阻R的電流保持恒定，電阻存在著漲落，引起電阻兩端電壓漲落，則有:

δV=1*δR(1-1)

Sv(f)=I2SR(f)I2(1-2)

假設(shè)功率譜密度是Pnf=Sf/f，相應(yīng)的電壓密度是e2nf=Af/f。從頻率f1到f2之間的總噪聲是

v2nf=f2f1Affdf=Aflog(f2f1)(1-3)

所以，1/f譜圖取決于上下的截止頻率，而不是絕對(duì)帶寬。

1.2 1/f噪聲的自相關(guān)函數(shù)

LTV(線性時(shí)不變)系統(tǒng)的輸出可以用卷積積分來表示：

(1-4)

式中t是輸出的觀察時(shí)間。

通常我們更關(guān)心的是計(jì)算系統(tǒng)輸出的均方根值，假設(shè)輸入是由一個(gè)靜態(tài)噪聲源驅(qū)動(dòng)的，為了得到這個(gè)結(jié)果，第一步是計(jì)算輸出均方根值

x~2=∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v）dudv(1-5)

這里f(t)被一個(gè)噪聲信號(hào)n(t)取代。均方值通過計(jì)算整個(gè)時(shí)間段的平均值可以得到：

x-2=limk→∞1k∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v)dudv

x-2=∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)［limk→∞1kni(t-u)ni(t-v)］dudv(1-7)

ni(t)表示的是i次抽樣函數(shù)。

式(1-7)中間括號(hào)的內(nèi)容是隨機(jī)輸入噪聲信號(hào)的自相關(guān)函數(shù)，所以可以簡(jiǎn)化成為

x-2∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)［Rnn(u-v)］dudv（1-8）

上面這個(gè)式子指出了自相關(guān)函數(shù)在時(shí)域噪聲分析中的重要性。在噪聲處理中最常見的兩種噪聲是白噪聲和1/f噪聲的自相關(guān)函數(shù)。理想的1/f噪聲是一個(gè)多樣的過程，它的自相關(guān)函數(shù)無法用閉合的形式來簡(jiǎn)單的表達(dá)出來。實(shí)際應(yīng)用中進(jìn)行了一些近似。

白噪聲的雙邊功率譜密度(double-sidedPSD)是

v2Δf=N02（1-9）

N0是單邊噪聲功率譜密度，單位是V2/HZ。自相關(guān)函數(shù)就是簡(jiǎn)單的功率譜密度的富利葉反變換。所以，白噪聲的自相關(guān)函數(shù)由沖擊響應(yīng)函數(shù)給出

Rnn(τ)=N02δ(τ)(1-10)

1/f噪聲，可以看成是白噪聲通過一個(gè)噪聲濾波器得到的結(jié)果。白噪聲驅(qū)動(dòng)的通過線性時(shí)間濾波器的輸出的自相關(guān)函數(shù)是

Ryy(τ)=N02∞-∞h(λ)h(τ+λ)dλ(1-11)

h(t)是成型濾波器的沖擊響應(yīng)。理想的1/f成型噪聲濾波器得沖擊響應(yīng)是

h(t)2fct，t＞0(1-11)

fc是轉(zhuǎn)角頻率。1/f噪聲的自相關(guān)函數(shù)從下式計(jì)算得到：

Ryy(τ)=N0fc∞01λ1λ+τdλ(1-13)

可以得到

Ryy=2N0fc1n［λ+λ+τ］∞0(1-14)

上面這個(gè)式子形式是非常簡(jiǎn)單的，可惜的是，式(1-14)是發(fā)散的，所以并不可能去精確的計(jì)算得到1/f噪聲的自相關(guān)函數(shù)。我們可以簡(jiǎn)化得到有理形式的濾波器（或者一串濾波器），解決方法是定義一個(gè)總的操作時(shí)間，并使用恰當(dāng)?shù)慕?，?1-14)可以轉(zhuǎn)變?yōu)?/p>

Ryy(top，τ)=2N0fc1n［λ+λ+τ］top-τδ(1-15)

假設(shè)topτδ，可以得到

Ryy(top，τ)≌N0fc［1n(4)+1n(topfc）-1n(τfc) ］ （1-16）

注意，top和δ，是等同于時(shí)域內(nèi)的兩個(gè)頻率限fmin(等同于/top)和(fc等同于1/δ)，所以，整個(gè)1/f噪聲過程包括了N個(gè)10倍頻段，(1-12)可以寫成

N=log10(topfc)（1-17）

Ryy(N，τ)≌N0fc［1n(4)+2.3N-1n(τfc)］(1-18)

在1/f噪聲仿真的實(shí)際應(yīng)用中，式(1-14)比式(1-19)更實(shí)用。

1.3 1/f噪聲的產(chǎn)生方法

產(chǎn)生1/f噪聲的最直觀的方法就是將白噪聲通過噪聲濾波器過濾之后得到。但是，理想的1/f濾波器并不是有理的，所以需要開發(fā)一個(gè)可以近似得到1/f噪聲的濾波器。需解決的主要問題是求這一近似的濾波器的算法。

使用下圖所示簡(jiǎn)化的噪聲濾波器。

圖1 仿真1/f噪聲成型濾波器的示意圖

一系列的hn(t)進(jìn)行疊加，可以得到一系列的濾波器，產(chǎn)生了完整的1/f噪聲的輸出，傳輸函數(shù)如下給出

H(s)=1+∑N+1n=010-0.5n22πfcs+2π10-nfc(1-19)

N定義了總頻率帶寬的10倍頻數(shù)量，fc是轉(zhuǎn)角頻率。1/f噪聲的數(shù)字化產(chǎn)生需要將模擬濾波器變成數(shù)字化濾波器，z域傳輸函數(shù)很容易計(jì)算得到，計(jì)算的結(jié)果如下式：

H(z)=1+∑Nn=010-0.5n2πfcTs1-z-1exp(-2π10-nfcTs)(1-20)

現(xiàn)在可以使用Matlab內(nèi)置的filter函數(shù)來進(jìn)行處理了。在仿真產(chǎn)生1/f噪聲時(shí)，比較合理N的值大約為10。這意味著一個(gè)帶寬為1MHz的系統(tǒng)，最低頻率在10-4左右。這樣，總的仿真時(shí)間長(zhǎng)度需要在s。一般來說104s，仿真N個(gè)10倍頻段需要大約10N個(gè)點(diǎn)。如果N>6，考慮到對(duì)計(jì)算機(jī)內(nèi)存和資源的要求，是不現(xiàn)實(shí)的。在實(shí)際應(yīng)用中，通常選定采樣頻率為107，模擬時(shí)間長(zhǎng)度為10-3?？紤]到所需要的精度大概需要104-105個(gè)點(diǎn)。為了解決這個(gè)問題我們定義一個(gè)參量Neff

Neff=1+log10(tsimTs)(1-21)

tsim是總的仿真時(shí)長(zhǎng)。然后整個(gè)數(shù)字濾波器被分為兩個(gè)部分，如圖2所示，所有的nNeff子濾波器，按照式(1-20)來計(jì)算，所有的n＞Neff部分，輸出在整個(gè)仿真過程中幾乎都是常量。意味著高階的子濾波器可以進(jìn)行簡(jiǎn)單的將數(shù)值（偏移量）相加。每一部分所得均方根等于穩(wěn)態(tài)時(shí)子濾波器的輸出均方根。

圖2 數(shù)字濾波器可以被Neff分為兩個(gè)部分

1.4 產(chǎn)生白噪聲并計(jì)算其功率譜密度

1.4.1 使用matlab仿真產(chǎn)生白噪聲

如上節(jié)分析，因?yàn)榉抡娈a(chǎn)生1/f噪聲的方法是白噪聲通過一個(gè)濾波器，所以先需要仿真得到一個(gè)白噪聲。(程序較簡(jiǎn)單略)

在Matlab中產(chǎn)生白噪聲，可以使用randn函數(shù)的得到一個(gè)具有高斯正態(tài)分布的白噪聲

noise=rms×randn(1，npts)

產(chǎn)生得到一個(gè)均值為0，均方根值為rms，總點(diǎn)數(shù)為npts的噪聲，均方根值由白噪聲底階和采樣頻率共同決定，如下式得到：

rms=N02Ts(1-22)

N0是白噪聲的單邊功率譜密度，Ts是采樣頻率。采樣頻率最好取為10的因子并且略大于系統(tǒng)帶寬，這樣在整個(gè)帶寬內(nèi)白噪聲都有一個(gè)平坦的功率譜。

程序輸入輸出量說明：

function noise=shot_noise(Ts，N0，npts)

Ts采樣時(shí)間，fs=1/Ts采樣頻率，N0白噪聲的單邊功率譜密度，npts仿真所取的點(diǎn)，noise為所得到的噪聲數(shù)值序列

圖3 仿真得到的白噪聲時(shí)域波形

1.4.2 計(jì)算功率譜密度的三種不同的方法

函數(shù)中使用direct關(guān)鍵詞則使用了直接法進(jìn)行計(jì)算，又稱周期圖法，它是把隨機(jī)序列x(n)的N個(gè)觀測(cè)數(shù)據(jù)視為一能量有限的序列，直接計(jì)算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實(shí)功率譜的估計(jì)。

使用indirect關(guān)鍵詞，使用間接法先由序列x(n)估計(jì)出自相關(guān)函數(shù)R(n)，然后對(duì)R(n)進(jìn)行傅立葉變換，得到x(n)的功率譜估計(jì)。

其中complex關(guān)鍵詞，綜合方法會(huì)將加入三種窗函數(shù)(Hamming，Kaiser，Chebyshev)的圖譜估計(jì)直接顯示在一幅圖上，比較直觀和明顯。

1.4.3 白噪聲的功率譜密度

為了驗(yàn)證產(chǎn)生的白噪聲的性能，計(jì)算它的PSD進(jìn)行驗(yàn)證

function psd_s(Ts，N0，npts，method)

% PSD of shot noise

除了method外與產(chǎn)生程序的變量意義相同，method為產(chǎn)生功率譜密度算法參數(shù)。

圖4 仿真得到的白噪聲功率譜

1.5 仿真產(chǎn)生1/f噪聲

function noise=f1_noise(tsim，Ts，fc，N0，N，npts)

% noise=f1_noise(tsim，Ts，fc，N0，N，npts)

% variable declaration

% Ts: 采樣時(shí)間長(zhǎng)度

% N:10倍頻數(shù)輛

% tsim:總仿真時(shí)長(zhǎng)

% fs:采樣周期

% fc:轉(zhuǎn)角頻率

具體源程序代碼見2.1.1節(jié)

圖5 仿真得到的1/f噪聲時(shí)域波形

tsim，仿真持續(xù)時(shí)間，Ts采樣時(shí)間，N10倍頻段，fs采樣頻率，npts總共仿真的點(diǎn)數(shù)。

得到了一個(gè)雜亂無章的1/f噪聲時(shí)域波型，需要驗(yàn)證它就是我們所需要的。

1.6 1/f噪聲的功率譜密度

function psd_f(tsim，Ts，fc，N0，N，npts，method)

% pad_f(tsim，Ts，fc，N0，N，npts，method)

% variable declaration

% Ts: Sampling Period

% N:frequecncy decades

% tsim:total simulated time interval

% fs:Sampling Frequency

% fc:corner frequncy

具體的源程序代碼見2.1.2節(jié)。

除了method外與產(chǎn)生程序的變量意義相同，method為產(chǎn)生功率譜密度算法參數(shù)。

direct用直接法進(jìn)行計(jì)算，直接法又稱周期圖法，它是把隨機(jī)序列x(n)的N個(gè)觀測(cè)數(shù)據(jù)視為一能量有限的序列，直接計(jì)算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實(shí)功率譜的估計(jì)。

ndirect，間接法由序列x(n)估計(jì)出自相關(guān)函數(shù)R(n)，對(duì)R(n)進(jìn)行傅立葉變換，得到x(n)的功率譜估計(jì)。

其中complex方法，會(huì)將加入三種窗函數(shù)的圖譜估計(jì)直接顯示在一幅圖上。

圖6 仿真得到的1/f噪聲功率譜密度

可以看到，PSD近似與lnf成一個(gè)反比的關(guān)系，隨著f的增大而減小。對(duì)圖形進(jìn)行處理，橫軸以ln（1/f）為坐標(biāo)，下降的關(guān)系是成直線的，在低頻率范圍內(nèi)出現(xiàn)的扭曲，是濾波器傳遞函數(shù)在N＞Neff時(shí)的近似造成的。

圖7 對(duì)數(shù)坐標(biāo)-1/f功率譜密度

取對(duì)數(shù)坐標(biāo)作圖(見圖7)，舍棄超過轉(zhuǎn)角頻率的點(diǎn)，進(jìn)行線性擬合得到：b=-1.4486 a=-2.2333r=-3.9922

2 仿真過程的實(shí)現(xiàn)

2.1 仿真產(chǎn)生1/f噪聲

2.1.1 仿真產(chǎn)生1/f噪聲的程序

function noise=f1_noise(tsim，Ts，fc，N0，N，npts)

% noise=f1_noise(tsim，Ts，fc，N0，N，npts)

% variable declaration

% Ts: Sampling Period

% N:frequecncy decades

% tsim:total simulated time interval

% fs:Sampling Frequency

% fc:corner frequncy

Neff=1+log10(tsim/Ts);%parameter Neff

noise_tmp=shot_noise(Ts，N0，npts);

fs=1/Ts;

% Generation white noise using the shot_noise function

fn=0;

for n=0:N

if n<=Neff

b=［0 ((10^(-0.5*n))*sqrt(2)*pi*fc/fs)］;

a=［1-1*exp(-2*pi*(10^(-n))*fc/fs)］;

fn=fn+filter(b，a，noise_tmp);

else

fn=fn+sqrt(N0*pi*fc/4)*randn(1，1);

end

end

plot(fn);pause;

noise=fn;

2.1.2 1/f噪聲功率譜密度計(jì)算程序

function psd_f(tsim，Ts，fc，N0，N，npts，method)

% pad_f(tsim，Ts，fc，N0，N，npts，method)

% variable declaration

%Ts: Sampling Period

%N:frequecncy decades

%tsim: total simulated time interval

%fs: Sampling Frequency

%fc: corner frequncy

% methods

% 'direct'

% 'indirect'

% 'periodogram'

% 'Kaiser'

% 'Chebyshev'

% 'complex'

noise=f1_noise(tsim，Ts，fc，N0，N，npts);

Fs=1/Ts;

n=N0;

xn=noise;

%direct method

switch lower(method)

case{'direct'}

window= boxcar(length(noise));

nfft=1024;

［Pxx，f］=periodogram(noise，window，nfft，F(xiàn)s);

plot(f，10*log10(Pxx))

%indirect method

case{'indirect'}

nfft=1024;

window=boxcar(length(n));

noverlap=0;

p=0.9;

［Pxx，Pxxc］=psd(xn，nfft，F(xiàn)s，window，noverlap，p);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k，plot_Pxx);

pause;

figure(2)

plot(k，［plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc］);

case {'periodogram'}

% Instantiate spectrum object and call its PSD method.

h=spectrum.periodogram('rectangular');

hopts=psdopts(h，xn);% Default options based on the signal x

set(hopts，'Fs'，F(xiàn)s，'SpectrumType'，'twosided'，'CenterDC'，true);

psd(h，xn，hopts)

case{'kaiser'}

h=spectrum.welch('hamming'，64);

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

W=hpsd1.Frequencies;

h.WindowName='Kaiser';

hpsd2=psd(h，xn，'Fs'，F(xiàn)s);

Pxx2=hpsd2.Data;

hpsd=dspdata.psd(Pxx2，W，'Fs'，F(xiàn)s);

plot(hpsd);

legend('kaiser');

case{'chebyshev'}

h=spectrum.welch('hamming'，64);

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

W=hpsd1.Frequencies;

h.WindowName='Chebyshev';

hpsd3=psd(h，xn，'Fs'，F(xiàn)s);

Pxx3=hpsd3.Data;

hpsd=dspdata.psd(Pxx3，W，'Fs'，F(xiàn)s);

plot(hpsd);

legend('Chebyshev');

case{'complex'}

h= spectrum.welch('hamming'，64);

% Create three power spectral density estimates.

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

Pxx1=hpsd1.Data;

W=hpsd1.Frequencies;

h.WindowName='Kaiser';

hpsd2=psd(h，xn，'Fs'，F(xiàn)s);

Pxx2=hpsd2.Data;

h.WindowName='Chebyshev';

hpsd3=psd(h，xn，'Fs'，F(xiàn)s);

Pxx3=hpsd3.Data;

% Instantiate a PSD data object and store the three different

% estimates since they all share the same frequency information.

hpsd=dspdata.psd(［Pxx1， Pxx2， Pxx3］，W，'Fs'，F(xiàn)s);

plot(hpsd);

legend('Hamming'，'kaiser'，'Chebyshev');

end

3 結(jié)論

用產(chǎn)生白噪聲加濾波器的方法比較容易仿真產(chǎn)生1/f噪聲，并且仿真得到的1/f噪聲較真實(shí)。

參考文獻(xiàn)

［1］Tae-Hoon Lee， Gyuseong Cho Chap1. Monte Carlo based time-domain Hspice noise simulation for CSA-CRRC cricuit［J］. Nuclear Instruments and Methods in Physics Research Section A， 2003，505:328-333.

［2］Terry，S.C. Blalock， B.J. Rochelle， J.M. Ericson， M.N. Caylor， S.D. Time-domain noise analysis of linear time-Invariant and linear time-variant systems using MATLAB and HSPICE［C］. Nuclear Science， IEEE Transactions on 2005，52:805-812.

［3］王愛萍，王惠南.基于小波分析的1/f噪聲降噪 ［J］.數(shù)據(jù)采集與處理，2006.

［4］文 俊，陳良棟.1/f噪聲源以及MOS器件中的1/f噪聲［J］.電子器件，1996.