离散傅里叶变换DFFT求频谱和功率谱

请教各位专业人士：

用下面的代码生成数据，其中y(:,2)即为需要用于傅里叶变换的时间信号。
请问:如何在MATLAB中写代码，将该时间信号转变为频谱和功率谱？不同的采样时间长度是否会影响结果？（程序中的结束时刻、时间间隔均可任意修改）

参考公式：
由傅里叶变换转为频谱：（其中的y(t)为时间信号，y*(Ω)为频谱）

img


频谱转为功率谱：（其中的Sy(Ω)为功率谱）

img

clear;
clc;

%参数
I=5.37*10^10;         
delta=1.47*10^8;    
R1=3.17;
R3=-1.61;
R5=-2.76;         
g=9.81;
xs=0.875;
r1=delta*g*R1/I;
r3=delta*g*R3/I;
r5=delta*g*R5/I;
fwind=0.01144;   

%二阶常微分方程求解
odefun2 = @(t,y)[y(2);-(r1*y(1)+r3*y(1)^3+r5*y(1)^5)+fwind*(1-y(1)^2/2+y(1)^4/24)]; 
ts = 0; %计算开始时间
te = 1000; %计算结束时间
h = 0.01;  %时间间隔
t = ts:h:te; %时间段
y = zeros(length(t),2); 
y(1,:)=[xs,0]; %初值

for n = 1:length(t)-1
k1 = odefun2(t(n),y(n,:)');
k2 = odefun2(t(n)+h/2,y(n,:)'+h/2*k1);  
k3 = odefun2(t(n)+h/2,y(n,:)'+h/2*k2);
k4 = odefun2(t(n)+h,y(n,:)'+h*k3);
y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';
end

plot(t,y(:,1),t,y(:,2));

如能解答，万分感激！

该回答引用GPTᴼᴾᴱᴺᴬᴵ

在MATLAB中，将时间信号转换为频谱和功率谱可以使用FFT和pwelch函数。下面是一个示例程序，假设有一个时间信号x，结束时刻为t_end，时间间隔为dt：

% 生成时间序列
t = 0:dt:t_end;

% 生成信号
f1 = 10;  % 信号频率1
f2 = 20;  % 信号频率2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);

% 计算频谱
N = length(x);
Fs = 1/dt;  % 采样频率
X = fft(x)/N;
f = (0:N-1)*(Fs/N);

% 绘制频谱图
figure;
plot(f,abs(X));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Frequency Spectrum');

% 计算功率谱
window = hann(N);  % 汉宁窗
[Pxx,f] = pwelch(x,window,[],[],Fs);

% 绘制功率谱图
figure;
plot(f,Pxx);
xlabel('Frequency (Hz)');
ylabel('Power');
title('Power Spectrum');

在上面的代码中，首先生成了一个时间序列t和一个信号x（这里的信号是两个正弦波的叠加）。然后使用FFT函数计算了信号的频谱，并绘制了频谱图。最后使用pwelch函数计算了信号的功率谱，并绘制了功率谱图。其中，pwelch函数使用了汉宁窗，可以改变窗口大小、重叠率等参数以适应不同的应用场景。

该回答引用ChatGPT

要将信号转换为频谱和功率谱，可以使用MATLAB中的fft和pwelch函数。

首先，使用fft函数计算信号的频谱，代码如下：


% 计算频谱
Fs = 1/h; % 采样频率
L = length(y); % 信号长度
Y = fft(y(:,2)); % 计算信号的傅里叶变换
P2 = abs(Y/L); % 计算双侧频谱
P1 = P2(1:L/2+1); % 取单侧频谱
P1(2:end-1) = 2*P1(2:end-1); % 按照MATLAB官方文档的建议修正幅度
f = Fs*(0:(L/2))/L; % 构建频率向量
figure;
plot(f,P1);
xlabel('频率 (Hz)');
ylabel('幅度');

该代码将信号y(:,2)进行fft变换，得到信号的频谱P1和对应的频率向量f。这里使用了MATLAB的fft函数默认的参数，因此得到的频率向量是以Hz为单位的。如果需要将频率转换为角频率，可以将f乘以2π。

接下来，可以使用pwelch函数计算信号的功率谱密度，代码如下：


% 计算功率谱密度
N = length(y); % 段数，与信号长度相等
window = hamming(round(N/4)); % 汉明窗函数
noverlap = round(N/8); % 重叠点数
[Pxx,f] = pwelch(y(:,2),window,noverlap,N,Fs); % 计算功率谱密度
figure;
plot(f,Pxx);
xlabel('频率 (Hz)');
ylabel('功率谱密度');

该代码将信号y(:,2)进行功率谱密度计算，得到信号的功率谱密度Pxx和对应的频率向量f。这里使用了汉明窗函数，并指定了重叠点数，以便得到更平滑的估计值。

不同的采样时间长度可能会影响结果。在信号的时间域和频率域中，采样时间长度都会影响信号的分辨率。较短的采样时间长度可以提供更高的时间分辨率，但频率分辨率较低。相反，较长的采样时间长度可以提供更高的频率分辨率，但时间分辨率较低。因此，在选择采样时间长度时，需要权衡时间分辨率和频率分辨率的需求，以获得最合适的结果。