Calculate a PSD using a time series
I have the following code and the sample data set (attached).
clear; close all; clc
load("sampleData.mat")
whos
dt = sampleData.time(2) – sampleData.time(1); % Time step
n = length(sampleData.eta); % Number of data points
fs = 1 / dt; % Sampling frequency (Hz)
eta_fft = fft(sampleData.eta) / n; % Normalized FFT
% Single-sided Power Spectral Density
PSD_single_sided = (2 * abs(eta_fft(1:round(n/2))).^2) * fs;
f_single_sided = (0:round(n/2)-1) * (fs/n); % Frequency array in Hz
omega_fft = 2 * pi * f_single_sided; % Convert to rad/s
omega = sampleData.omega;
psd_org = sampleData.S;
% Plot the PSD against the spectrum
figure;
plot(omega, psd_org, ‘b’, ‘LineWidth’, 1.5); % Original spectrum
hold on;
plot(omega_fft, PSD_single_sided, ‘r’, ‘LineWidth’, 1.5); % Single-sided PSD from time series
xlabel(‘omega (rad/s)’);
ylabel(‘Spectral Density (m^2 s)’);
legend(‘direct Spectrum’, ‘PSD from Time Series’);
xlim([0 5])
grid on;
I have followed the theretical steps so that the two PSD should match, maybe not perfectly but at least to a good extent. But I’m not getting a decent match, and I think I might be doing the PSD calculations wrong. Signal processing is not my strongest suit, so could anyone tell me what I’m doing wrong?
I also looked into Power Spectral Density Estimates Using FFT, but I’m having a hard time getting the results in m^2s. Technically, the psd_org and the calculated psd should match.I have the following code and the sample data set (attached).
clear; close all; clc
load("sampleData.mat")
whos
dt = sampleData.time(2) – sampleData.time(1); % Time step
n = length(sampleData.eta); % Number of data points
fs = 1 / dt; % Sampling frequency (Hz)
eta_fft = fft(sampleData.eta) / n; % Normalized FFT
% Single-sided Power Spectral Density
PSD_single_sided = (2 * abs(eta_fft(1:round(n/2))).^2) * fs;
f_single_sided = (0:round(n/2)-1) * (fs/n); % Frequency array in Hz
omega_fft = 2 * pi * f_single_sided; % Convert to rad/s
omega = sampleData.omega;
psd_org = sampleData.S;
% Plot the PSD against the spectrum
figure;
plot(omega, psd_org, ‘b’, ‘LineWidth’, 1.5); % Original spectrum
hold on;
plot(omega_fft, PSD_single_sided, ‘r’, ‘LineWidth’, 1.5); % Single-sided PSD from time series
xlabel(‘omega (rad/s)’);
ylabel(‘Spectral Density (m^2 s)’);
legend(‘direct Spectrum’, ‘PSD from Time Series’);
xlim([0 5])
grid on;
I have followed the theretical steps so that the two PSD should match, maybe not perfectly but at least to a good extent. But I’m not getting a decent match, and I think I might be doing the PSD calculations wrong. Signal processing is not my strongest suit, so could anyone tell me what I’m doing wrong?
I also looked into Power Spectral Density Estimates Using FFT, but I’m having a hard time getting the results in m^2s. Technically, the psd_org and the calculated psd should match. I have the following code and the sample data set (attached).
clear; close all; clc
load("sampleData.mat")
whos
dt = sampleData.time(2) – sampleData.time(1); % Time step
n = length(sampleData.eta); % Number of data points
fs = 1 / dt; % Sampling frequency (Hz)
eta_fft = fft(sampleData.eta) / n; % Normalized FFT
% Single-sided Power Spectral Density
PSD_single_sided = (2 * abs(eta_fft(1:round(n/2))).^2) * fs;
f_single_sided = (0:round(n/2)-1) * (fs/n); % Frequency array in Hz
omega_fft = 2 * pi * f_single_sided; % Convert to rad/s
omega = sampleData.omega;
psd_org = sampleData.S;
% Plot the PSD against the spectrum
figure;
plot(omega, psd_org, ‘b’, ‘LineWidth’, 1.5); % Original spectrum
hold on;
plot(omega_fft, PSD_single_sided, ‘r’, ‘LineWidth’, 1.5); % Single-sided PSD from time series
xlabel(‘omega (rad/s)’);
ylabel(‘Spectral Density (m^2 s)’);
legend(‘direct Spectrum’, ‘PSD from Time Series’);
xlim([0 5])
grid on;
I have followed the theretical steps so that the two PSD should match, maybe not perfectly but at least to a good extent. But I’m not getting a decent match, and I think I might be doing the PSD calculations wrong. Signal processing is not my strongest suit, so could anyone tell me what I’m doing wrong?
I also looked into Power Spectral Density Estimates Using FFT, but I’m having a hard time getting the results in m^2s. Technically, the psd_org and the calculated psd should match. fft, psd MATLAB Answers — New Questions
