matlab实验-经典谱估计（自相关法）

内容：计算并画出下述信号的功率谱

x(t)=y(t)+A_1*sin(2*π*f_1*t)+A_2*sin(2*π*f_2*t)+A_3*sin(2*π*f_3*t)

其中，归一化频率f_1=0.1,f_2=0.28,f_3=0.3,y(t)是白噪声u(t)经过FIR滤波器的输出，滤波器参数为B=[1,-0.1,0.09,0.648].

要求：采用自相关法，画出x(t)的功率谱。

说明： 1.产生白噪声，功率为0.1(randn),通过FIR滤波器产生y(t)(filter)

            2.产生x(t)

            3.画信号的功率谱（freqz）

            4.用自相关法画x(t)的功率谱(fft)

            5.用Welch方法画x(t)的功率谱

MATLAB代码实现：

% MATLAB 滤波器设计和信号处理

% 参数设置
A1 = 1.0; A2 = 1.0; A3 = 1.0; % 振幅
f1 = 0.1; f2 = 0.28; f3 = 0.3; % 归一化频率
B = [1, -0.1, 0.09, 0.648];    % 滤波器系数
N = 1024;                      % 信号点数
fs = 1;                        % 采样频率
t = (0:N-1) / fs;              % 时间向量

% 生成白噪声
rng(0); % 随机数生成器的种子，确保可重复性
u = randn(N, 1) * sqrt(0.1);

% 通过FIR滤波器滤波
y = filter(B, 1, u);

% 生成信号x(t)
x = y + A1 * sin(2 * pi * f1 * t') + A2 * sin(2 * pi * f2 * t') + A3 * sin(2 * pi * f3 * t');

% 使用freqz计算FIR滤波器的频率响应
[h, w] = freqz(B, 1, 8000, fs);
figure;
subplot(1, 1, 1);
plot(w, 20 * log10(abs(h)));
title('FIR滤波器的频率响应');
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');


% 使用FFT计算x(t)的功率谱密度（自相关法）
Rxx = xcorr(x, 'biased');
psd_fft = abs(fft(Rxx, N));
f_fft = linspace(0, fs/2, N/2);
figure;
subplot(1, 1, 1);
plot(f_fft, psd_fft(1:N/2));
title('x(t)的功率谱密度 - FFT（自相关法）');
xlabel('频率 (Hz)');
ylabel('功率');

% 使用Welch方法计算x(t)的功率谱密度
[pxx, f] = pwelch(x, hamming(256), 128, 256, fs);
figure;
subplot(1, 1, 1);
plot(f, pxx);
title('x(t)的功率谱密度 - Welch方法');
xlabel('频率 (Hz)');
ylabel('功率');

运行结果：