How to calculate appropriate DFT without typical sampling frequency?

2 次查看（过去 30 天）
Paul Hinze
Paul Hinze 2022-7-8
回答： Sudarsanan A K 2023-10-13
Dear community,
I am not good in signal processing. I need to calculate a single sided power spectral density of a signal, which is called 'g' in my code, but i am not able to get it completely correct (Sampling frequency issues i gues).
%% Investigate window function
clear
close
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Params %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tleng = 50; % simulated time length [ms]
dt = 0.01; % time step [ms]
% time vector
t = (0:round(Tleng/dt)-1)*dt; % time vector
spike = zeros(size(t));
spike(1) = 1; % Make amplitude
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Generate alpha function %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau = 3.2; Amp = 1; window_shape = {'alpha'};
g = Synapse (spike, dt, tau, Amp, window_shape);
subplot(2,1,1)
plot ( t , g )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% DFT (Internet example which is probalbly not correct) %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs=100000; %sampling frequency
NFFT = 1024;
X = fft(g,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(NFFT*NFFT);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
subplot(2,1,2)
plot(f/1000,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (kHz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Function for Synapse %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function g = Synapse (spikes, dt, tau, Asyn, window_shape)
% This synapse is build to stimulate the IF-Model.
% *** INPUT PARAMETERS ***
% spikes: : spike vector (to show the number of inputs at each step)
% dt : size of time step [ms]
% c : input amplitude [coulomb]
% tau : time constant [ms]
% Asyn : synaptic input amplitude
% window_shape : determination of window
%
% *** OUTPUT PARAMETERs ***
% g : simulated synaptic input waveform [nS]
% Select window shape
window_shape = cell2mat(window_shape);
switch window_shape
case 'alpha'
% creating output vector
g = zeros(1,length(spikes));
% initial values
x = 0;
y = 0;
% factors used for calculation
format long % Genrate more digits
eulers_number = exp(1); % Generate eulers number
format short % Reduce number of digits for better speed
tau = tau / eulers_number; % Standardise with factor. This factor is based on eulers number
aa = Asyn * exp(1.0) / tau; % amplitude factor for normalizing the input
ee = exp(-dt/tau); % decrease factor used at each time step
% step-by-step calculation
for t = 1:length(spikes)
y = ee*y + aa*spikes(t);
x = ee*x + dt*y;
g(t) = x; % store data
end
end
end

请先登录，再回答此问题。

回答（1 个）

Sudarsanan A K
Sudarsanan A K 2023-10-13
Ran in:
Hello Paul,
I understand that you are trying to calculate a single-sided power spectral density (PSD) of a signal and finding difficulty in understanding the relevance of sampling frequency in that. From the code that you provided; I also note that the PSD is computed using the FFT method.
The sampling frequency is inversely dependent on the time step or sampling interval. A smaller time step leads to a higher sampling frequency, providing finer frequency resolution and better capturing of high-frequency components. A larger time step results in a lower sampling frequency, leading to a coarser frequency resolution and potential loss of high-frequency information. Therefore, the choice of the time step affects the accuracy and resolution of the PSD estimation.
Please find the code attached below, which demonstrates the estimation of PSD using the FFT method.
% Define the sampling time
Ts = 0.001; % Sampling time of 1 ms
% Calculate the sampling frequency
fs = 1 / Ts; % Sampling frequency
% Generate a simple sinusoid signal
t = 0:Ts:1; % Time vector (1 second duration)
f = 10; % Frequency of the sinusoid (Hz)
x = sin(2*pi*f*t); % Sinusoidal signal
% Calculate the FFT-based PSD
N = length(x); % Length of the signal
X = fft(x); % Compute the FFT
X = X(1:(N+1)/2); % Select one-sided spectrum
Pxx_fft = (1/(fs*N)) * abs(X).^2; % Calculate the PSD
Pxx_fft(2:end-1) = 2*Pxx_fft(2:end-1); % Double the amplitude of non-DC and non-Nyquist frequencies to conserve the total power
frequencies_fft = 0:fs/length(x):fs/2; % Frequency axis for single-sided spectrum
% Plot the time domain signal and the PSDs
figure;
% Plot the time domain signal
subplot(2,1,1);
plot(t, x);
grid on;
xlabel('Time (s)');
ylabel('Amplitude');
title('Time Domain Signal');
% Plot the PSD
subplot(2,1,2);
plot(frequencies_fft, pow2db(Pxx_fft));
grid on;
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
title('Power Spectral Density');
Feel free to use your own signal in place of the sinusoid for the PSD estimation.
Additionally, you can refer to the following link for better understanding the PSD estimation using FFT method:
https://mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html
Further, there are MATLAB functions that ease your job for estimating the PSD, which can be found at:
I hope this helps you to understand the PSD estimation.

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by