How to calculate appropriate DFT without typical sampling frequency?

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

[Paul Hinze]

This is how it looks at the moment

Answer 1

[Sudarsanan A K]

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:

• "pwelch()": https://mathworks.com/help/signal/ref/pwelch.html
• "periodogram()": https://mathworks.com/help/signal/ref/periodogram.html

I hope this helps you to understand the PSD estimation.