3
$\begingroup$

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:

Sample signal with N=192

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.

Amplitude spectrum comparison

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);
$\endgroup$
  • $\begingroup$ Your plots show that you are downsampling or smoothing somehow, not just interpolating. Interpolation always produces more data points. Please clarify this part of your procedure. $\endgroup$ – Dave Kielpinski Sep 15 '17 at 1:29
  • $\begingroup$ Thanks for that question. I am indeed downsampling, but unintentionally. Effectively, when I can drop every other bin, that's what Matlab is doing. Otherwise it averages certain bins and drops others. Looks even worse than what I thought I was doing wrong. How should I proceed? $\endgroup$ – ishouldknowtheanswerbutidont Sep 15 '17 at 14:52
  • $\begingroup$ I can probably answer your question, but it depends on the details of your downsampling procedure. I understand you are using Matlab interp1 but is there more to it? Can you post this part of your code? $\endgroup$ – Dave Kielpinski Sep 15 '17 at 18:18
  • $\begingroup$ It's quite simple. The Fourier coefficients go from X[1] through X[N/2 + 1], so I just scale the sequence of points 0..50 accordingly and call interp1 to query X on the scaled points, then add 1 because Matlab index from 1. But again, I am not required to use interp1 or anything similar. I just want to estimate the frequency content at specific bins, regardless of the length of the series. $\endgroup$ – ishouldknowtheanswerbutidont Sep 18 '17 at 12:06
  • $\begingroup$ About interpolation of the frequency data (not saying it is the right path to take), see: Question: Interpolation of magnitude of discrete Fourier transform (DFT) $\endgroup$ – Olli Niemitalo Sep 18 '17 at 12:26
2
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.