plotting power spectrum density

I have 2 .dat file with me one of amplitude and one of time and I want to do the powerspectrum analysis of the data. Can anyone please help me in creating a time vector for doing FFT

Answer 1

[Hassaan]

% Load the data from the .dat files
amplitude = load('amplitude.dat');
time = load('time.dat');

% Check if time data is provided, otherwise create time vector
if isempty(time)
    Fs = 1000; % Sampling frequency in Hz
    Ts = 1 / Fs; % Sampling interval in seconds
    time = (0:length(amplitude)-1) * Ts; % Time vector
end

% Compute FFT
L = length(amplitude); % Length of the signal
Y = fft(amplitude); % Compute the FFT
P2 = abs(Y / L); % Two-sided spectrum
P1 = P2(1:L/2+1); % Single-sided spectrum
P1(2:end-1) = 2 * P1(2:end-1);

f = Fs * (0:(L/2)) / L; % Frequency vector

% Calculate Power Spectrum Density (PSD)
PSD = (1 / (Fs * L)) * (abs(Y(1:L/2+1)).^2);
PSD(2:end-1) = 2 * PSD(2:end-1);

% Plot the PSD
figure;
plot(f, 10 * log10(PSD));
title('Power Spectrum Density');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
grid on;

Notes

1. Sampling Frequency (Fs): Ensure that you know the correct sampling frequency of your data. If it is not 1000 Hz, replace Fs = 1000 with the correct value.
2. Units: The PSD is plotted in decibels (dB/Hz). The conversion is done using 10 * log10(PSD).

-----------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[Kashish]

I already have the time

[Kashish]

I already have these files and I have used this code for extraction of data

fileinfo = dir('amp-A-016.dat');

num_samples = fileinfo.bytes/2;

fid = fopen('amp-A-016.dat', 'r');

v = fread(fid, num_samples, 'int16');

fclose(fid);

v = v * 0.0000374;

%%%Time Data input

fileinfo = dir('time.dat');

num_samples = fileinfo.bytes/4; % int32 = 4 bytes

fid = fopen('time.dat', 'r');

t = fread(fid, num_samples, 'int32');

fclose(fid);

% t = t / frequency_parameters.amplifier_sample_rate; % sample rate from header file

t=t/20000;

[Kashish]

what would be the time vector if I already have the time data

please help me

[Hassaan]

% Load the amplitude data
fileinfo = dir('amp-A-016.dat');
num_samples = fileinfo.bytes / 2;
fid = fopen('amp-A-016.dat', 'r');
amplitude = fread(fid, num_samples, 'int16');
fclose(fid);
amplitude = amplitude * 0.0000374;

% Load the time data
fileinfo = dir('time.dat');
num_samples = fileinfo.bytes / 4; % int32 = 4 bytes
fid = fopen('time.dat', 'r');
time = fread(fid, num_samples, 'int32');
fclose(fid);
time = time / 20000; % Adjust according to your sample rate

% Compute the sampling frequency (Fs) from time data
Fs = 1 / mean(diff(time)); % Sampling frequency in Hz

% Compute FFT
L = length(amplitude); % Length of the signal
Y = fft(amplitude); % Compute the FFT
P2 = abs(Y / L); % Two-sided spectrum
P1 = P2(1:L/2+1); % Single-sided spectrum
P1(2:end-1) = 2 * P1(2:end-1);
f = Fs * (0:(L/2)) / L; % Frequency vector

% Calculate Power Spectrum Density (PSD)
PSD = (1 / (Fs * L)) * (abs(Y(1:L/2+1)).^2);
PSD(2:end-1) = 2 * PSD(2:end-1);

% Plot the PSD
figure;
plot(f, 10 * log10(PSD));
title('Power Spectrum Density');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
grid on;

[Kashish]

thank you the code is working what if I want to do multitaper powerspectrum density method