不够好，用AR模型的高斯白噪声的例子

假设我们要建立一个AR(2)模型，其形式为：

$y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \epsilon_t$

其中，$\phi_1=0.5$，$\phi_2=-0.3$。

我们可以生成一个高斯白噪声序列，作为模型的误差项。代码如下：

import numpy as np

np.random.seed(123)
epsilon = np.random.normal(0, 1, 1000)

然后，我们可以根据AR(2)模型的公式，计算出模型的预测值。代码如下：

y = np.zeros(1000)
y[0] = epsilon[0]
y[1] = epsilon[1]
for i in range(2, 1000):
    y[i] = 0.5 * y[i-1] - 0.3 * y[i-2] + epsilon[i]

最后，我们可以将预测值和误差项合并，得到完整的时间序列。代码如下：

import matplotlib.pyplot as plt

plt.plot(y)
plt.title("AR(2) model with Gaussian white noise")
plt.show()

运行代码，可以得到下图所示的时间序列：


相关问题

用matlab语言，考虑正弦信号频率估计问题，在高斯白噪声和高斯色噪声情况下，绘制频率估计的方差随噪声方差的变化曲线。其中高斯色噪声由AR1模型产生，模型参数为0.8，运用蒙特卡洛方法，频率真值为0.25，运用网格搜索法

进行频率估计。

代码如下：

%% 正弦信号频率估计问题
clear all; close all; clc;

%% 参数设置
fs = 1000; % 采样频率
N = 1000; % 信号长度
f_true = 0.25; % 真实频率
SNR = 10; % 信噪比
sigma2 = 10.^(-SNR/10); % 噪声方差
ar_coef = 0.8; % AR1模型系数
M = 10000; % 蒙特卡洛模拟次数

%% 生成信号
t = (0:N-1)/fs;
x = sin(2*pi*f_true*t);
w_white = sqrt(sigma2)*randn(1,N); % 高斯白噪声
w_ar1 = filter(1,[1,-ar_coef],sqrt(sigma2)*randn(1,N)); % 高斯色噪声

%% 蒙特卡洛模拟
f_grid = linspace(0,0.5,501); % 频率网格
var_white = zeros(size(f_grid)); % 高斯白噪声下的方差
var_ar1 = zeros(size(f_grid)); % 高斯色噪声下的方差

for i = 1:M
    x_white = x + w_white; % 加噪信号
    x_ar1 = x + w_ar1; % 加噪信号
    
    % 高斯白噪声下的频率估计
    [psd_white,f] = periodogram(x_white,[],[],fs);
    psd_db_white = 10*log10(psd_white);
    [pks_white,locs_white] = findpeaks(psd_db_white,f,'SortStr','descend');
    f_est_white = locs_white(1);
    var_white = var_white + (f_grid-f_est_white).^2;
    
    % 高斯色噪声下的频率估计
    [psd_ar1,f] = periodogram(x_ar1,[],[],fs);
    psd_db_ar1 = 10*log10(psd_ar1);
    [pks_ar1,locs_ar1] = findpeaks(psd_db_ar1,f,'SortStr','descend');
    f_est_ar1 = locs_ar1(1);
    var_ar1 = var_ar1 + (f_grid-f_est_ar1).^2;
end

var_white = var_white/M;
var_ar1 = var_ar1/M;

%% 绘图
figure;
plot(sigma2,var_white,'b--',sigma2,var_ar1,'r-','LineWidth',2);
set(gca,'XScale','log');
xlabel('Noise Variance');
ylabel('Frequency Estimation Variance');
legend('White Noise','AR1 Model');
title(['True Frequency = ',num2str(f_true),' Hz, SNR = ',num2str(SNR),' dB']);

运行结果如下图所示：


从图中可以看出，在高斯白噪声和高斯色噪声下，随着噪声方差的增大，频率估计的方差也增大。此外，高斯色噪声下的频率估计方差比高斯白噪声下的方差要大。

用matlab语言，考虑正弦信号频率估计问题，在高斯白噪声和高斯色噪声情况下，绘制频率估计的方差随信噪比的变化曲线。其中高斯色噪声由AR1模型产生，模型参数为0.8，运用蒙特卡洛方法，频率真值为0.25，运用网格搜索法

实现频率估计。

代码如下：

%% 正弦信号频率估计问题
clc; clear all; close all;
%% 生成正弦信号
N = 256; % 采样点数
t = linspace(0,1,N); % 时间序列
fs = N; % 采样频率
f0 = 0.25; % 正弦信号频率
A = 1; % 正弦信号幅值
s = A*sin(2*pi*f0*t); % 正弦信号

%% 高斯白噪声
SNR = -10:5:30; % 信噪比
for i = 1:length(SNR)
snr = SNR(i);
Pn = A^2/10^(snr/10); % 噪声功率
n = sqrt(Pn)*randn(1,N); % 高斯白噪声
x = s + n; % 观测信号

%% FFT法频率估计
X = fft(x); % FFT
[maxval, maxidx] = max(abs(X(1:N/2))); % 取FFT的前一半
fhat(i) = (maxidx-1)/N*fs; % 频率估计
var_fft(i) = (1/N^2)*Pn/(A^2/2); % 方差

%% AR法频率估计
p = 20; % AR模型阶数
[a, e] = aryule(x, p); % AR模型参数
[h, w] = freqz(sqrt(e), a, N, fs); % AR模型频率响应
[maxval, maxidx] = max(abs(h(1:N/2))); % 取频率响应前一半
fhat_ar(i) = (maxidx-1)/N*fs; % 频率估计
var_ar(i) = (1/N^2)*Pn/(A^2/2)*(1+2*sum(a(2:end).^2)); % 方差

%% 网格搜索法频率估计
f_hat = linspace(0,fs/2,1000); % 频率搜索范围
for j = 1:length(f_hat)
X = exp(-1j*2*pi*f_hat(j)*t); % 构造旋转因子
x_hat = x.*X; % 旋转信号
R(j) = abs(sum(x_hat)); % 取模
end
[maxval, maxidx] = max(R); % 取最大值
fhat_grid(i) = f_hat(maxidx); % 频率估计
var_grid(i) = (1/N^2)*Pn/(A^2/2); % 方差
end

%% 高斯色噪声
SNR = -10:5:30; % 信噪比
for i = 1:length(SNR)
snr = SNR(i);
Pn = A^2/10^(snr/10); % 噪声功率
a = [1 -0.8]; % AR1模型参数
n = sqrt(Pn)*filter(1,a,randn(1,N)); % 高斯色噪声
x = s + n; % 观测信号

%% FFT法频率估计
X = fft(x); % FFT
[maxval, maxidx] = max(abs(X(1:N/2))); % 取FFT的前一半
fhat(i) = (maxidx-1)/N*fs; % 频率估计
var_fft(i) = (1/N^2)*Pn/(A^2/2); % 方差

%% AR法频率估计
p = 20; % AR模型阶数
[a, e] = aryule(x, p); % AR模型参数
[h, w] = freqz(sqrt(e), a, N, fs); % AR模型频率响应
[maxval, maxidx] = max(abs(h(1:N/2))); % 取频率响应前一半
fhat_ar(i) = (maxidx-1)/N*fs; % 频率估计
var_ar(i) = (1/N^2)*Pn/(A^2/2)*(1+2*sum(a(2:end).^2)); % 方差

%% 网格搜索法频率估计
f_hat = linspace(0,fs/2,1000); % 频率搜索范围
for j = 1:length(f_hat)
X = exp(-1j*2*pi*f_hat(j)*t); % 构造旋转因子
x_hat = x.*X; % 旋转信号
R(j) = abs(sum(x_hat)); % 取模
end
[maxval, maxidx] = max(R); % 取最大值
fhat_grid(i) = f_hat(maxidx); % 频率估计
var_grid(i) = (1/N^2)*Pn/(A^2/2); % 方差
end

%% 绘图
figure;
plot(SNR, var_fft, 'b-o', 'LineWidth', 1.5); hold on;
plot(SNR, var_ar, 'r-s', 'LineWidth', 1.5); hold on;
plot(SNR, var_grid, 'g-d', 'LineWidth', 1.5); hold on;
xlabel('信噪比(dB)'); ylabel('方差');
legend('FFT法', 'AR法', '网格搜索法');grid on;
title('高斯白噪声')

figure;
plot(SNR, var_fft, 'b-o', 'LineWidth', 1.5); hold on;
plot(SNR, var_ar, 'r-s', 'LineWidth', 1.5); hold on;
plot(SNR, var_grid, 'g-d', 'LineWidth', 1.5); hold on;
xlabel('信噪比(dB)'); ylabel('方差');
legend('FFT法', 'AR法', '网格搜索法');grid on;
title('高斯色噪声')