burg演算法

『壹』 Burg演算法對任意平穩隨機過程是否仍有效

是會有效的

『貳』 Olderburg-Parker法

（一）計算單一密度界面的Olderbury-Parker法

R.L.Parker和P.W.Oldenburg（1972年，1974年）提出一種快速計算界面的正反演方法。目前，該方法應用較廣泛。

1.單一界面異常的波譜

如圖11-5所示，單一密度界面在水平面（xoy面）內的引力位為

圖11-5 單一界面快速反演法圖示

勘探重力學與地磁學

引力位的波譜為

勘探重力學與地磁學

Erdelyi積分表（1954年）可知，1/（x²＋y²＋ζ²）^1/2的波譜為2πe^-ωζ/ω，根據傅氏變換（FT）的位移性質可知，1/［（x-ξ）²＋（y-η）²＋ζ²］^1/2的波譜應為

，因此引力位的波譜為

勘探重力學與地磁學

ξ和ζ的積分限是-∞～∞，而ζ的積分限是（h₀＋Δh）～h₀，其中Δh為界面相對於平均深度h₀的起伏，界面在h₀以上Δh為負，在h₀以下Δh為正，而Δh又是ξ和ζ的函數，故

勘探重力學與地磁學

又已知重力異常波譜與引力位波譜之間的關系為S_Δg＝ωS_V，所以重力異常的波譜為

勘探重力學與地磁學

因為

！，而（-Δhw）⁰/（0！）＝1，所以

勘探重力學與地磁學

式中：（S_Δh）ⁿ為（Δh）ⁿ的波譜。上式即為單一密度界面的重力異常波譜公式。

2.計算界面深度的迭代法

將（11-30）式改寫成

勘探重力學與地磁學

從上式解出

勘探重力學與地磁學

在已知界面平均深度h₀和剩餘密度ρ的條件下，可用迭代法計算界面深度。先忽略式中（S_Δh）ⁿ的高次項，S_Δh的一級近似為

勘探重力學與地磁學

將

代入（11-30）式計算S_Δg的一級近似值：

勘探重力學與地磁學

計算S_Δh的修改量：

勘探重力學與地磁學

而S_Δh的二級近似值為

勘探重力學與地磁學

再將

代入（11-30）式求出S_Δg的二級近似值

，如此反復迭代下去，直到n次迭代後S_Δg與

的均方誤差小於給定的允許誤差為止。對S_Δh計算逆富氏變換求出界面起伏Δh，而界面深度為

勘探重力學與地磁學

在具體計算時，尚需加一適當的正則化因子，以防止高波數因子

的發散作用。

（二）計算單一磁性界面的Olderburg-Parker方法

根據磁位與引力位的關系——泊松公式的推導可知，均勻磁化體的磁性為

勘探重力學與地磁學

式中：V為引力位；l_s，m_s，n_s為磁化強度M的方向餘弦（即cosα_s，cosβ_s，cosv_s）。根據傅氏變換（FT）的線性和微分性質、磁位的波譜應為

勘探重力學與地磁學

式中：q_s＝j（l_scosθ＋m_ssinθ）＋n_s；w_x＝wcosθ；w_y＝wsinθ；w_x，w_y分別為與x軸、y軸對應的圓波數；θ為徑向圓波數w的極角。（11-38）式是物體引力位與磁位波譜間的關系式，也稱波數域的泊松公式。

由ΔT的物理意義和位函數的方向導數性質可知：

勘探重力學與地磁學

式中：∂T為地磁場T₀方向的長度元，l_t＝cosα_t，m_t＝cosβ_t，n_t＝cosγ_t為T₀的方向餘弦。依FT的線性和微分性質可知ΔT的波譜應為

勘探重力學與地磁學

勘探重力學與地磁學

式中：q_t＝j（l_tcosθ＋m_tsinθ）＋n_t。又根據前面介紹S_Δg＝wS_V，

（-w)^n-1/n！，所以

勘探重力學與地磁學

令q_s＝q_t＝1時，便得到垂直磁化垂直磁異常ΔZ_⊥的波譜公式：

勘探重力學與地磁學

下面介紹計算磁性界面深度的迭代法：

將（11-41）式改寫成

勘探重力學與地磁學

這樣

由（11-42）式出發，在已知界面的平均深度h₀和磁化強度M的條件下，可用迭代法計算界面深度。先忽略式中（S_Δh）ⁿ的高次項，S_Δh的一級近似值為

勘探重力學與地磁學

將

代入（11-41）式，計算

的一級近似值：

勘探重力學與地磁學

計算修正量：

勘探重力學與地磁學

而S_Δh的二級近似值為

勘探重力學與地磁學

再將

代入（11-41）求出

的二級近似值

。如此反復迭代下去，直到幾次迭代後

與

的均方誤差小於給定的允許誤差為止。對S_Δh計算逆FT，求出界面起伏Δh，而界面深度為

勘探重力學與地磁學

在具體計算時，尚需加一適當的正則化因子，以防高波數因子

的發散作用。

由以上所介紹的磁性界面的計算方法可見，利用的是垂直磁化垂直磁異常

或其波譜

。但是實際磁異常ΔZ或ΔT往往是斜磁化的，所以，在計算界面之前，須將斜磁化的ΔZ或ΔT轉換成垂直磁化垂直磁異常ΔZ_⊥或其波譜

。

『叄』 請問matlab中如何在功率譜密度圖中添加一條置信曲線

這么多,看都懶得看!~

『肆』 如何用matlab實現去除語音信號中雜訊

感覺應該用濾波器。
舉個例子，如果類似就可套用，如果不對就當借鑒了。
正弦信號x=sin50*pi*t加入了隨機雜訊，采樣頻率為1000Hz，濾除掉其中30Hz的雜訊。
程序：
clear
t=(0:0.001:0.5);
x=sin(50*pi*t)+rand(size(t));
[b,a]=butter(10,30/500);
y=filter(b,a,x);
plot(t,x);
figure,plot(t,y);
%end

『伍』 信號功率的求法

能量有限平均功率為零的信號是能量信號，如單位沖擊信號；
能量無限，平均功率有限的信號是功率信號。
題中信號有一個常量10，平均功率不為零且有限，能量無限，屬於"功率信號「
功率為：10^2+4^2+2^2=120

『陸』 從經典譜估計到現代譜估計

諧波分析最早可追溯到古代對時間的研究。18世紀伯努利（Bernoulli）、歐拉（Euler）和拉格朗日（Lagrange）等人對波動方程及其正弦解進行了研究，19世紀初葉，傅立葉（Fourier）證明了在有限時間段上定義的任何函數都可以用正弦和餘弦分量的無限諧波的總和來表示。1898年舒斯特（Schuster）以傅立葉分析為基礎來擬合待分析信號，研究太陽黑子數的周期變化，並提出了周期圖的概念。1930年維納（Wiener）發表了經典性論文《廣義諧波分析》，對平穩隨機過程的自相關函數和功率譜密度作了精確的定義，證明了二者之間存在著傅氏變換（以下簡稱傅氏變換）的關系，從而為功率譜分析奠定了堅實的統計學基礎。由於1934年辛欽（Khintchine）也獨立地證明了自相關函數和功率譜之間的傅氏變換關系，即維納-辛欽（Wiener-Khintchine）定理。根據這個定理（詳見第一章），平穩離散隨機信號x（n）的自相關函數r_xx（m）

r_xx（m）=E［x（m+n）x^*（n）］ （4-1）

與功率譜P_xx（e^jω）之間構成一傅氏變換對，即

地球物理信息處理基礎

若x（n）還是各態遍歷性的，則其自相關函數可由它的一個采樣時間序列用時間平均的方法求出，即

地球物理信息處理基礎

在大多數應用中x（n）是實信號，於是上式可寫成

地球物理信息處理基礎

實際上，一般只能在時域觀測到隨機信號的有限個采樣值（例如N個值），可表示為

x_N（n）=｛x（0），x（1），…，x（N-1）｝=｛x（n），n=0，1，…，N-1｝

其自相關函數只能由這N個采樣數據進行估計，常用有偏估計

地球物理信息處理基礎

這是一種漸近一致估計，稱之為采樣自相關函數。

用采樣自相關函數的傅氏變換作為功率譜的估計，這種方法是布萊克曼（R.Blackman）和杜基（J.Tukey）在1958年提出來的，稱為功率譜估計的自相關法（簡稱BT法）。此方法需要先求出有限個觀測數據估計自相關函數，然後再根據式（4-2）計算出功率譜。在快速傅氏變換（FFT）演算法提出之前，這是一種最流行的功率譜估計方法。

1965年庫利（Cooley）和杜基（Tukey）完善了著名的FFT演算法，使計算傅氏變換的速度提高了兩個數量級，運算量顯著降低，這樣DFT變換很快在各領域，特別是在工程實踐中得到了廣泛應用。由式（4-5）知，

為x（n）與x（-n）的卷積運算，因為

地球物理信息處理基礎

若x（n）的傅氏變換為X（e^jω），則x（-n）的傅氏變換等於X^*（e^jω）。對式（4-5）兩端取傅氏變換，得到

地球物理信息處理基礎

這表明：通過對隨機數據直接進行離散傅氏變換，然後取其幅值的平方，再對多樣本進行此種運算並取平均值作為功率譜的估計，即舒斯特的周期圖，這種譜估計受到了人們的普遍重視，因為它不需要計算自相關函數，而直接計算功率譜。

周期圖和自相關法以及它們的改進方法稱為功率譜估計的經典方法，周期圖和自相關法是經典功率譜估計的兩個基本方法。由於FFT的出現，周期圖和自相關法往往被結合起來使用，其步驟如下：

（1）對x_N（n）補N個零，求

；

（2）由

作傅氏變換，得

，這時｜m｜≤M=N-1；

（3）對

加窗函數v（m），這時｜m｜≤M＜＜N-1，得

；

（4）利用

，求

的傅氏變換，即

地球物理信息處理基礎

由周期圖法得到的功率譜

，其估計方差並不隨樣本長度的增加而趨於零，

，出人意料的是，不管數據記錄有多長，周期圖和自相關法得到的估計都不是功率譜的良好估計。事實上，隨著記錄長度增加，這兩種估計的隨機起伏反而會更加嚴重！此外，它們存在著以下兩個難以克服的固有缺點。

（1）頻率解析度（區分兩個鄰近頻率分量的能力）不高。因為它們的頻率解析度（赫茲）反比於數據記錄長度（秒）（即Δf=k/T_p=k/NT，k為常數，T_p=NT為數據的記錄長度，T為采樣周期），而實際應用中一般不可能獲得很長的數據記錄，即觀察到的數據只能是有限個，而觀察不到的數據被認為是0。這樣，如果只有N個觀測數據，而對於N以外的數據，信號仍有較強的相關性，那麼估計出的功率譜就會出現很大的偏差。

（2）對於有限的觀測數據，相當於將信號在時域內乘以矩形窗函數，因而在頻域內則相當於使真正的功率譜與sinc函數進行卷積，由於sinc函數不同於δ 函數，它有主瓣和旁瓣，這樣使卷積後的功率譜不同於真正的功率譜。sinc函數的主瓣不是無限窄的，引起功率譜向附近頻域擴展，造成譜的模糊，降低譜的解析度；同時，由於sinc函數的旁瓣存在，導致能量向旁瓣中「泄漏」（稱之為旁瓣泄漏），即引起頻譜間的干擾，信號強的功率譜旁瓣影響信號弱的功率譜檢測，嚴重時，會使主瓣產生很大失真，檢測不出弱信號，或者把旁瓣誤認為是信號，造成假信號。為了對經典功率譜估計進行改進，可以採用各種不同的窗函數，但其結果都是以增加主瓣寬度來換取旁瓣的壓低，因此功率譜解析度低是經典功率譜估計的致命缺點。

為了克服以上缺點，人們曾做過長期努力，提出了平均、加窗平滑等辦法，在一定程度上改善了經典功率譜估計的性能。實踐證明，對於長數據記錄來說，以傅氏變換為基礎的經典功率譜估計方法，的確是比較實用的。但是，經典方法始終無法根本解決頻率解析度和功率譜估計穩定性之間的矛盾，特別是在數據記錄很短的情況下，這一矛盾顯得尤為突出。這就促進了現代功率譜估計方法研究的展開。

現代功率譜估計方法主要是以隨機過程（Stochastic Process）的參數模型（Parameter Model）為基礎的，稱之為參數模型方法。雖然說現代功率譜估計技術的研究和應用主要起始於60年代，但實際上，時間序列模型在非工程領域早已被採用，如Yule在1927年、Walker在1931年都曾使用過自回歸模型預測描述經濟的時間序列的發展趨勢，而Prony則早在1795年就曾採用指數模型去擬合在氣體化學實驗中獲得的數據。在統計學和數值分析領域中，人們也曾採用過模型方法。

現代功率譜估計的提出主要是針對經典功率譜估計（周期圖和自相關法）的解析度和方差性能不好的問題而提出的。1967年Burg在地震學研究中受到線性預測濾波的啟發，提出了最大熵譜估計方法，在提高解析度方面作了最有意義的探索。1968年Parzen正式提出了自回歸譜估計方法。1971年Van der Bos證明了一維最大熵譜估計與自回歸譜估計等效。1972年出現的譜估計的Prony方法在數學上與自回歸方法有某些類似。目前以自回歸滑動平均模型為基礎的譜估計已經比自回歸模型譜估計具有更高的頻率解析度和更好的性能。1973年Pisarenko提出的諧波分解方法提供了可靠的頻率估計方法。1981年Schmidt提出了譜估計的多信號分類（MUSIC）演算法等。因此，現代功率譜分析主要有ARMA譜分析、最大似然、熵譜估計和特徵分解四種方法。ARMA譜分析是一種建模方法，即通過平穩線性信號過程建立模型來估計功率譜密度；熵譜估計包括最大熵譜和最小交叉法；特徵分解也叫特徵構造法和子空間法，包括Pisarenko諧波分解法、Prony法、MUSIC法和ESPRIT法（用旋轉不變技術估計參數方法）。

現代功率譜估計研究仍側重於一維功率譜分析，而且大部分是建立在二階矩（相關函數、方差、功率譜密度）基礎上的。但由於功率譜密度是頻率的實函數，缺少相位信息，因此，建立在高階譜基礎上的譜估計方法正引起人們的注意，特別是雙譜估計和三譜估計的研究受到了高度的重視。其它如多維譜估計、多通道譜估計等的研究也正在發展中。人們希望這些新方法能更多地在提取信息、估計相位和描述非線性等方面獲得應用。

『柒』 輸入字元不是 MATLAB 語句或表達式中的有效字元.

分號有問題 你的是中文下的分號，要英文下的分號

『捌』 利用Burg演算法計算AR模型的參數：arburg函數

（1）a=arburg（x，order）；返回利用Burg演算法計算得到的AR模型的參數：輸入參數or-der為AR模型的階數。

（2）［a，e］=arburg（…）；當輸入信號為白雜訊時，同時返回最後誤差e。

（3）［a，e，k］=arburg（…）：同時返回反射系數k。

［例4-5］用Burg演算法計算AR 模型的參數。

rand（『seed』，0）；a=［1 0.1 0.2 0.3 0.4 0.5 0.6］；

x=impz（1，a，20）+randn（20，1）/20；

［a，e，k］=arburg（x，5）。

『玖』 short-sequence是什麼意思

short-sequence
一個短序列；一個短串列
例句篩選
1.
Research on short sequence power spectrum estimates based onthe Burg algorithm
基於Burg演算法的短序列譜估計研究

2.
Even different strains of the same kind of bacteria will have uniqueidentifying patterns of a given short sequence.
甚至同種類型細菌的不同菌株也將有一個既定的短序列的唯一的識別模式。

『拾』 matlab提示 function [psdviaBurg, f, p] = myBurg(x, Fs, varargin) | Error: Function definitions

1、將下面函數復制保存為myBurg.m文件。就可以計算。預測誤差E 、系數a 、誤差功率psdviaBurg

function [psdviaBurg, f, p,a,E] = myBurg(x, Fs, varargin)
%MYBURG 根據burg演算法實現的AR模型功率譜計算
% psdviaBurg 根據burg演算法求出的功率譜值
% f 頻率軸參數
% p 模型階次
% a AR模型參數
%E AR模型誤差
% x 輸出信號
% Fs 采樣率
% varargin 若為數值型，則為AR模型階次
% 若為字元串，則為定階准則，AR模型階次由程序確定
%
% $Author: lskyp
% $Date: 2010.6.26
% 解析輸入參數內容
error(nargchk(3, 3, nargin)); % 該函數的輸入必須為三個個
if strcmp(class(varargin{1}), 'double')
p = varargin{1};
elseif ischar(varargin{1})
criterion = varargin{1};
else
error('參數2必須為數值型或者字元串');
end
x = x(:);
N = length(x);
% 模型參數求解
if exist('p', 'var') % p變數是否存在，存在則不需要定階，直接使用p階
[a, E] = computeARpara(x, p);
else % p不存在，需要定階，定階准則即criterion
p = ceil(N/3); % 階次一般不超過信號長度的1/3

% 計算1到p階的誤差
[a, E] = computeARpara(x, p);

% 根據誤差求解目標函數最小值
kc = 1:p + 1;
switch criterion
case 'FPE'
goalF = E.*(N + (kc + 1))./(N - (kc + 1));
case 'AIC'
goalF = N.*log(E) + 2.*kc;
end
[minF, p] = min(goalF); % p就是目標函數最小的位置，也即定階准則給出的階次

% 使用p階重新求解AR模型參數
[a, E] = computeARpara(x, p);
end
% 計算功率譜密度
[h, f] = freqz(1, a, [], Fs);
psdviaBurg = E(end)*abs(h).^2./Fs;

function [a, E] = computeARpara(x, p)
% 根據信號序列x和階次p計算AR模型參數和誤差
N = length(x);
% 初始值
ef = x; % 前向預測誤差
eb = x; % 後向預測誤差
a = 1; % 初始模型參數
E = x'*x/N; % 初始誤差
k = zeros(1, p); % 為反射系數預分配空間，提高循環速度
E = [E k]; % 為誤差預分配空間，提高速度
for m = 1:p
% 根據burg演算法步驟，首先計算m階的反射系數
efm = ef(2:end); % 前一階次的前向預測誤差
ebm = eb(1:end - 1); % 前一階次的後向預測誤差
num = -2.*ebm'*efm; % 反射系數的分子項
den = efm'*efm + ebm'*ebm; % 反射系數的分母項
k(m) = num./den; % 當前階次的反射系數

% 更新前後向預測誤差
ef = efm + k(m)*ebm;
eb = ebm + conj(k(m))*efm;

% 更新模型系數a
a = [a; 0] + k(m)*[0; conj(flipud(a))];

% 當前階次的誤差功率
E(m + 1) = (1 - conj(k(m))*k(m))*E(m);
end

2、例如：
參考論壇中的例子。
randn('state', 1);
x = randn(100, 1);
y = filter(1, [1 1/2 1/3 1/4 1/5], x);
pburg(y, 4, [], 1000);
[psd, f, p,a,E] = myBurg(y, 1000, 'FPE');
figure;
a,E
plot(f, 10*log10(psd));

結果：圖還是一樣，就補貼出來了。
a =
1.0000
0.3116
0.3647
0.2086
0.2088
0.0425

E =
1.0224 0.9859 0.9167 0.8989 0.8644 0.8629