1

I want to find the median frequency from the FFT result. All I have is binned data I got from FFT (An array).

I'm doing this calculation in C, so I don't have access to the statistic tools from MATLAB or Python. I can't figure out the calculation to the median algorithm. How can I find the bin containing the (estimate) median frequency? Could you please guide me to the right direction? Thanks in advance.

spinmaster
  • 25
  • 2
  • Are you asking for a C function to compute the median of an array? Or how, once you have that function, to compute the median frequency? – Jdip Dec 08 '22 at 16:15
  • @Jdip. Not necessarily C function, I'm looking for the math or algorithm for that. I can figure out the median frequency once I got the which bin contains it. I can't figure out the median of the (binned) array. – spinmaster Dec 08 '22 at 16:17
  • Sorry I read the question slightly wrong.

    The median frequency $\text{MF}$ is the frequency that divides the power spectrum $P$ equally on both sides: $\sum_{j=1}^{\text{MF}}P_j = \sum_{j = \text{MF}}^{M}P_j = \frac{1}{2}\sum_{j=1}^{M}P_j$. It's not a completely trivial algorithm. When I have some time I'll try and explain here ;)

     – Jdip Dec 08 '22 at 17:00

1 Answers1

0

I'm assuming you know the sampling frequency of the data from which you computed the $\texttt{FFT}$. In which case here are the steps to compute the Median Frequency of the Power Spectrum of a signal $x[n]$ (you already have step 1. and possibly step 2. done):

  1. Compute the $\text{length}$-$N$ $\texttt{FFT}$ of $x[n]$ $$X_{\texttt{full}} = \texttt{FFT}(x[n])$$

  2. Keep the positive frequencies only to go from $X_{\texttt{full}}$ to $X_k$. Something like X = Xfull(0:floor(N/2))

  3. Compute the Power Spectral Density (or Power Spectrum, you'd just drop the $f_s$ term): $$P_k = \frac{2|X_k|^2}{f_sN}$$

  4. Compute the cumulative sum vector for $P_k$, call it $C_{k}$. There's no general expression for this, but you can do that easily with a loop. Inside the loop will be something like C(k) = C(k-1) + P(k)

  5. Compute the half power of $C_k$. There's 2 ways to get this: $$C^{\text{1/2}} = \frac{1}{2}\texttt{max}(C_k) = \frac{1}{2}C_{k=\texttt{last}}$$

  6. Find the index $k$ for which $C_k \leq C^{\text{1/2}} \leq C_{k+1}$

  7. Finally, map to a frequency in $\texttt{Hz}$: $$\texttt{Median Frequency} = kf_s/N \texttt{ (Hz)}$$

Jdip
  • 5,980
  • 3
  • 7
  • 29
  • Thank you so much for the effort but I'm having trouble understanding something. In the third step, isn't that supposed to be in a loop or is that not the regular power spectrum? Also what is the Xk symbolises in the third equation? That does not look like a known index in that step? – spinmaster Dec 09 '22 at 03:42
  • $X_k$ is the FFT array (only positive frequencies) and $k$ the indices. $P_k$ is also an array based on $X_k$. Whether it's computed in a loop or not depends on your implementation. In C you'd need a loop. Not in Matlab or Python. – Jdip Dec 09 '22 at 08:22
  • I implemented my own median frequency estimation in MATLAB and verified with the medfreq function. It looks easier now to do again in C. Thank you kind sir for everything. Your comments really made me understand the math. – spinmaster Dec 10 '22 at 13:46
  • 1
    @spinmaster Glad it helped! – Jdip Dec 10 '22 at 14:49