用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('高斯色噪声')