# Fourier Like Spectral Analysis with Uneven Intervals and Redesigned DFT Matrix

I intended to use a discrete Fourier transform (DFT) on a time series sampled at uneven intervals. What I did was to calculate a DFT matrix where the elements are the values at the uneven locations like this

{\displaystyle {\begin{aligned}X_{k}&= \frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}x_{n}\cdot e^{-{\frac {2\pi i}{time_{N-1} - time_0}}k \cdot time_n}\\&\end{aligned}}}

where $$time_n$$ is the time location in the sequence and $$time_{N-1} - time_0$$ is the entire time span of all samples.

I realize that orthogonality of the basis vectors in the DFT matrix might get somewhat lost. Is anything more lost? Can I do something more to make the transformation more comparable to an ordinary DFT?

Here is an implementation of it you can try with your own data http://bodavid.github.io/Be10/Be10/public_html/index.html Orthogonality between basis vectors are visualized in a colorized matrix.

• I hope you are aware of the procedures for DFT analysis of unevenly sampled data. There are various ways. Some people prefer NFFT. – Fat32 Jan 2 '19 at 23:13
• No, when I searched I only found recommendations to interpolate. – David Jonsson Jan 2 '19 at 23:18
• Do you want a practical code or theoretical end? (i.e., why did you try to formulate it by yourself ?) – Fat32 Jan 2 '19 at 23:26
• Both. I found it intuitive to adjust the DFT matrix to the points of measurements. – David Jonsson Jan 2 '19 at 23:37
• Ok. So if you need a library, go NFFT for a beginning. Or you can write your own from a few formulas (you can find nonuniform sampling and DFT here too). For a theoretical pursuit, notice that a mere re-locationing of the DFT basis is not sufficient for nonuniform to uniform conversion, as the previously hidden (located at nulling time positions) contributions from all other DFT basis must now also be included as they are no more nulled due to the uneven located data samples. – Fat32 Jan 2 '19 at 23:46

I've run up this alley a couple of times. Let's look at what you've actually accomplished.

1) You picked an interval/DFT frame

2) You built basis vectors with sinusoids the length of the frame, including only the time points where you have samples.

As you observed, this is no longer an orthonormal basis. This means you are going to have to solve a matrix inverse to get your actual coordinates. With a regular DFT, the matrix is the identity matrix so solving the inverse is implied. The coordinates you solve for are very much like bin values at this point.

3) Your basis vectors are actually sampled points from continuous sinusoidal functions. Therefore you can construct a continuous function which includes all your original sample points.

What have you actually accomplished? Well, a very good interpolation of your data.

If it so happens that your signal is composed of tones that have a whole number of cycles, your coordinates/bin values will be spot on. However, if the tone is not a whole number of cycles, the leakage of your uneven sampling will not match the leakage of evenly sampled points, just as leakage doesn't match for different sized DFTs.

So, at this point, if you want to do any serious spectral analysis of your signal using DFT techniques, I agree with the recommendations you have gotten, resample your reconstructed signal at evenly spaced points and work with that.

You will get a different interpolation function for every different interval length. You need to be aware of the "wraparound" nature of the DFT. This means the interpolation function will repeat every interval length. Therefore you will want to make your interval wider than just the width of your samples. When you resample you will likely want to narrow the interval to just around your points. The extrapolations are useless.

Of course, there are much simpler interpolation techniques.

I have derived the DFT for data whihc is sampled in a non uniform manner:

# The DFT Matrix for Non Uniform Time Samples Series

## Problem Statement

We have a signal $$x \left( t \right)$$ defined on the interval $$\left[ {T}_{1}, {T}_{2} \right]$$.
Assume we have $$N$$ samples of it given by $$\left\{ x \left( {t}_{i} \right) \right\}_{i = 0}^{N - 1}$$. The samples time $${t}_{i}$$ is arbitrary and not necessarily uniform.

We're after the DFT of the samples $$\left\{ X \left[ k \right] \right\}_{k = 0}^{K - 1}$$ as it was samples in a uniform manner (Implicitly means the samples in Frequency domain will be uniform as well).

## Deriving the Connection

In the DFT Transform the connection between time and frequency is given by:

$$x \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N} n } \tag{1}$$

In $$\eqref{EqnIdft}$$ we use $$n$$ for modeling the sample index in time. We usually build samples in time as $$x \left[ n \right] = x \left( n {T}_{s} \right)$$ where $${T}_{s}$$ is a uniform sampling interval.
Hence we could write:

$$x \left( n {T}_{s} \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} n {T}_{s}} \tag{2}$$

In $$\eqref{EqnIdft2}$$ we added explicit scaling of time. This is a known property of Fourier transform family which scales the domain in order to normalize the transform.

Now, there is nothing which blocks us from using arbitrary time:

\begin{align*} \tag{3} x \left( t \right) & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} t} && \text{} \\ & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k {F}_{s}}{N} t} && \text{Since {F}_{s} = \frac{1}{{T}_{s}} } \end{align*}

As can be seen $$\eqref{EqnIdft3}$$ makes sense as it goes through each element according to its frequency and sums to give the output at time $$t$$. We can go step farther and generalize it for cases we don't have uniform sampling frequency.
The average sampling frequency is given by $$\bar{F}_{s} = \frac{N}{ {T}_{2} - {T}_{1} }$$. Let's define $$T = {T}_{2} - {T}_{1}$$ and we'll get:

$$x \left( t \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} }$$

Which is many ways resembles the DTFT Transform equation which does the same in the other direction, transforming uniform discrete samples in time domain to arbitrary frequency (Within a frequency interval) in Frequency Domain:

\begin{align*} \tag{4} X \left( f \right) & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi f {T}_{s} n } && \text{} \\ & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi \frac{f}{ {F}_{s} } n } && \text{Since {F}_{s} = \frac{1}{{T}_{s}} } \end{align*}

We see the same scaling, $$\frac{f}{ {F}_{s} }$$ which scales the continuous $$f$$ relative to the interval of frequencies $${F}_{s}$$ which is equivalent to $$\frac{t}{ T }$$ which scales $$t$$ relative to the time interval of the continuous signal.

## The Transform Matrix

So, given the set of time indices $${\left\{ {t}_{i} \right\}}_{i = 0}^{N - 1}$$ the transformation matrix, from frequency domain to time domain, is given by:

$$D \in \mathbb{R}^{N \times K}, \; {D}_{i, k} = {e}^{ j 2 \pi k \frac{ {t}_{i} }{T} }$$

## The Model

In vector form the model is:

$$x = D y$$

Where $$y \in \mathbb{C}^{K}$$ is the vector of the frequency coefficients in uniform grid, $$x$$ is the samples in time (Non Uniform, Or at least no assumption of uniformity) and $$D$$ as defined above.
Since in our model we're after $$y$$ the answer is given by:

$$y = {D}^{\dagger} x$$

Where $${D}^{\dagger}$$ is the Pseudo Inverse Matrix of $$D$$.

## Implementation & Results

The code is as following:

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

samplingFrequency = 101; %<! [Hz]
samplingInterval = 1 / samplingFrequency; %<! [Sec]
startTime = 1; %<! [Sec]
endTime = 4; %<! [Sec]
timeInterval = endTime - startTime; %<! [Sec]

numSamples = round(samplingFrequency * timeInterval);
numSamplesTT = round(1.2 * numSamples);

signalFreq = 2; %!< [Hz]

% The uniform time grid
vT      = linspace(startTime, endTime, numSamples + 1);
vT(end) = [];
vT      = vT(:);

% The non uniform time grid - Reconstruction
vTT = endTime * rand(numSamplesTT, 1);
vTT = sort(vTT, 'ascend');

% The non uniform time grid - DFT
vTD = linspace(startTime, endTime, (10 * numSamples) + 1);
vTD(end) = [];
vTD = vTD(sort(randperm(length(vTD), numSamples)));
vTD = vTD(:);

% The uniform frequency grid
vF      = (samplingFrequency / 2) * linspace(-1, 1, numSamples + 1);
vF(end) = [];
vF      = vF(:);

vK = [-floor(numSamples / 2):floor((numSamples - 1) / 2)];
vK = vK(:);

%% Generate Data

vX  = cos(2 * pi * signalFreq * vT);
vFx = fftshift(fft(vX));

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = subplot(1, 2, 1);
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Reference Signal']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
'FontSize', fontSizeTitle);

hAxes           = subplot(1, 2, 2);
hStemObj = stem(vF, abs(vFx));
set(hStemObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['DFT of the Reference Signal']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
'FontSize', fontSizeTitle);

if(generateFigures == ON)
saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - Reconstruction

mD = exp(1j * 2 * pi * (vTT / timeInterval) * vK.') / numSamples;

% Reconstruction according to the model
vY = real(mD * vFx);

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTT, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - DFT of the Non Uniformly Sampled Data

vY  = cos(2 * pi * signalFreq * vTD);

mD = exp(1j * 2 * pi * (vTD / timeInterval) * vK.') / numSamples;
vFy = pinv(mD) * vY;

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTD, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hStemObj    = stem(vF, abs([vFx, vFy]));
set(hStemObj, 'LineWidth', lineWidthNormal);
% hLineSeries     = plot(vTT, vY);
% set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['DFT of the Uniform Signal & Non Uniform Signal']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end


Results are:

## Summary

In this post we derived how to estimate the Uniform DFT of a Non Uniform Time Series by solving linear system of equations.

The full code is available on my StackExchange Signal Processing Q32137 GitHub Repository.

## Remark: Why Do We Need to Apply fftshift() on the DFT of the Signal?

Indeed in the Reconstruction part we use fftshift(). The shallow answer is easy, we also build the vector vK as symmetric around zero.
But there is a deeper reason for that. In the DFT when we use uniform sampling in Frequency Domain and Time Domain Magic happens without us seeing it explicitly.

When we defined the term $$\frac{k}{ N {T}_{s} } n {T}_{s}$$ we replaces $$n {T}_{s}$$ with $$t$$ hence we prevent the term $${T}_{s}$$ to cancel itself. Now setting $${F}_{s} = \frac{1}{{T}_{s}}$$ means that we multiply by $$k$$ and we get frequencies which are out of the Nyquist Frequency.
In most cases when we that happens the Modulo property of the exponent comes in and we get the correct negative value of the frequency in the range $$\left[ -\pi, \pi \right]$$. Yet when $$t$$ is arbitrary we can think that $${F}_{s}$$ is changing per sample which means when we go farther than $$\pi$$ the modulo doesn't bring us to the correct answer.

First, as intuition, always think the DFT is defined on the $$\left[ -\pi, \pi \right]$$ interval and it is continuous. So as long as you work on this range things works as intended. This intuition can come from the Fourier Series and Discrete Fourier Series (DFS).

Let's try explaining it using a concrete example. Let's examine the exponent term from the derivation:

$$2 \pi \frac{k}{N {T}_{s}} n {T}_{S} = 2 \pi \frac{k}{{F}_{S}} \frac{{F}_{s}}{N} n = 2 \pi \frac{k b}{{F}_{s}} n$$

Where $$b$$ is the Bin Resolution in the Frequency domain. Now given the signal is:

$$x \left( t \right) = \cos \left( 2 \pi f t \right) \Rightarrow x \left( n {T}_{s} \right) = \cos \left( 2 \pi f {T}_{s} n \right) \Rightarrow x \left[ n \right] = \cos \left( 2 \pi \frac{f}{ {F}_{s} } n \right)$$

For $${F}_{s} = 100$$ [Hz] and $$N = 100$$ (Which means $$b = 1$$) we will have delta at $$k = 2$$ and $$k = 98$$. For $$k = 98$$:

$$2 \pi \frac{98}{{F}_{s}} n$$

This is clearly above the Nyquist frequency ($$\frac{{F}_{s}}{2}$$) and only for $${F}_{s} = 100$$ its modulo is $$-2$$ which is correct. But in the model above, since we have arbitrary $$t$$ one could think we have changing $${F}_{s}$$ which means we don't get the correct value.

This means the actual equation should be:

$$x \left( t \right) = \frac{1}{N} \sum_{k = \left \lfloor - \frac{K}{2} \right \rfloor }^{ \left \lfloor \frac{K - 1}{2} \right \rfloor } X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} }$$