%{
计算功率谱的四种方法，前三种方法采用归一化频率，第四种方法采用原始频率
%}

matlab">clear;clc;close all;
% 生成一个示例时间序列（例如正弦波加上一些随机噪声）
t = 0:0.01:100-0.1;  % 时间向量
% x = sin(t) + 0.5 * randn(size(t));  % 正弦波加噪声
x = sin(20*pi*t);%% 使用 pwelch 计算功率谱
% 输入信号 x，窗口大小为 256，重叠为 128，采样频率为 1 Hz
[pxx, f] = pwelch(x, 256, 128, [], 1);% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 Welch 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');%% 使用 periodogram 计算功率谱
[pxx, f] = periodogram(x, [], [], 1);% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 Periodogram 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');%% 直接使用fft
% 计算信号的 FFT
X = fft(x);% 计算功率谱
pxx = abs(X).^2 / length(x);  % 计算功率谱
f = (0:length(x)-1) / length(x);  % 频率轴% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 FFT 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');%% 使用自相关函数的方法
% 计算自相关函数
[R, lags] = xcorr(x, 'biased');
% 计算功率谱
N = length(x);  % 信号长度
Fs = 100;  % 采样频率% 获取自相关函数的一半（因为自相关函数是对称的）
R_half = R(N:end);% 计算自相关函数的傅里叶变换，得到功率谱
S = abs(fft(R_half));
f = (0:length(S)-1) * (Fs / length(S));  % 频率轴% 可选：将功率谱转换为dB
S_dB = 10 * log10(S);
figure;
plot(f, S_dB);  % 显示dB功率谱
xlabel('Frequency (Hz)');
ylabel('Power (dB)');
title('Power Spectrum using Autocorrelation');