How can I take a fixed number of bins after N-point DFT when N is unknown?

Question

I am working with machine learning for time series classification. I am trying to extract features from the amplitude spectrum.

My current concern is that I cannot tell the length of the signal in advance. Here's one example signal with $N=192$ observations, corresponding to approximately 2 seconds of data collection:

I do not have an upper bound on the length of the signal, although more than 1,000 observations usually indicate I have an outlier.

My approach has been to take the input signal $x=(x_1,x_2,\ldots,x_N)$ and then find the one-sided amplitude spectrum from its N-point DFT:

• $\mathcal{F}(x) = (X_1, X_2, \ldots, X_N)$
• $A(X)=(|X_1|, |X_2|, \ldots, |X_{\frac{N}{2}+1}|$)

Given that my sampling rate is approximately $f_s=100\,\text{Hz}$, I figured I would interpolate the one-sided amplitude spectrum to 51 points. (PS: I made a mistake here; please see clarification below.)

This procedure gives me a 192-point DFT and a one-sided amplitude spectrum with 97 coefficients, which I interpolate to 51 (linear interpolation). For comparison, here's the original amplitude spectrum and its interpolation. The horizontal axis is given in Hz.

They don't look very different. However, the features I take from the spectrum are such as short-term energy, variance, standard deviation, number of peaks, and location of the first peak, as well as the magnitude of specific coefficients (am I right in saying that's the energy at those bins?). So different approaches in finding those 51 coefficients can lead to quite different results.

Also, I have the impression that interpolating the amplitude spectrum is just plain wrong. Looking at the figure, it seems obvious that the interpolation made me lose a lot of energy content in the range near 0-5 Hz.

For disclosure, I did check each related question suggested by Stack Overflow as I wrote mine. The most similar question to mine seems to be "How to combine bins of my DFT", however I believe I have several additional doubts---to start with, I'm eager for criticism on whether I am doing the right approach in taking the absolute value of the Fourier coefficients and then trying to reduce the number of bins, or whether I should zero-pad my signal to a number of observations that is a multiple of 50 (and thus hopefully avoid leakage? I don't know, I'm really struggling to understand the subject and I might be working with several misconceptions here; any corrections will be appreaciated)

---

PS: I made a mistake in my assumptions. I incorrectly assumed Matlab's interp1 function interpolates over a set of query points and averages the values in-between. What it is in fact doing is downsampling.

---

PPS: here's the code for coefficient reduction:

qp = (0:50) * (numel(X) - 1) / (50) + 1;
X50 = interp1(X, qp);

Answer 1

Selection of feature vectors is a bit of an art when a sufficient statistic isn't obvious. If you want bins corresponding to the bandwidth of an equivalent fixed DFT.

I suggest you perform the STFT over your time duration varying signal, and then average the magnitude or if you want to avoid square roots the magnitude squared, over successive frames. It will preserve "energy". This does complicate things by adding issues such as windowing and overlapping but sometimes having more parameters to tweak provides flexibility.

One way of looking at feature set selection is:

P. M. Baggenstoss, "Class-specific feature sets in classification," in IEEE Transactions on Signal Processing, vol. 47, no. 12, pp. 3428-3432, Dec 1999. doi: 10.1109/78.806092 keywords: {probability;signal classification;statistical analysis;Class-specific feature sets;dimensionality problem;likelihood ratios;probabilistic classifiers;signal classification;sufficient statistics;Bayesian methods;Signal analysis;State estimation;Statistics;Training data;White noise}, URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=806092&isnumber=17473

Answer 2

You're almost there, you just need to be consistent in computing the discrete data values from your interpolating function. Suppose you have a continuous function $X(f)$. Let's say you want to get the discrete approximation over $N$ bins, and $f$ runs from $0$ to $f_{max}$. Then \begin{equation} X_k = \frac{1}{\Delta f} \int_{(k-1/2) \Delta f}^{(k+1/2) \Delta f} X(f) df \end{equation} You don't really have $X(f)$ of course. But you can approximate the integral on the RHS by using one of Matlab's numerical integration routines over your interpolating function.

How is this different from what you did? Well, near the 1-2 Hz peak, simply taking the downsampled DFT value misses the peak, as you mentioned. Taking the integral over the entire bin width does include this part of the signal, so you will get a lower but wider peak.

Finally, check your total signal energy. Parseval's theorem states that the Fourier transform preserves total signal energy. Consider signal points $x_i$ with $i=0 \ldots N-1$, and write the DFT as $X_i$. Then \begin{equation} \sum_{i=0}^{N-1} |x_i|^2 = \frac{1}{N} \sum_{i=0}^{N-1} |X_i|^2 \end{equation} So, any consistent resampling procedure has to keep the RHS of this equation constant. Simple interpolation and downsampling won't do the job for sure, but the procedure I've suggested should.