Right algorithm for fourier transform on physical heights

Question

I have data from a LIDAR unit that I would like to get the spectral density of. Unfortunately, the only thing I remember from my Fourier analysis class are the methods that I know will not work.

The data comes from a 1D LIDAR scan of a (mostly flat) surface, which returns radial distance at evenly spaced $d\theta$'s. I can convert this data to x-y data which looks distinct depending on the type of surface (e.g. grass looks a lot choppier than concrete). I think Fourier analysis would be a good way to distinguish between the different types. Unfortunately the data has some problems that makes it difficult to analyze:

• It is unevenly spaced in x.
• It has some missing values.
• It is not strictly a function in x (if there is some overhang, an earlier part of the sweep might hit under the overhang to get a distant x value, and a later part could strike the overhang, to get a closer value.)

It has occurred to me to do a Fourier analysis on the original $r-\theta$ data, but since I am concerned with the deviations around perfectly flat ground, I am not sure how I could use it.

Answer 1

Your question seems to have a couple of nested issues: First off, the computation of the Power-Spectral-Density is as straight forward as the computation of the signals' Discrete Fourier Transform, (DFT), followed by its absolute magnitude squared. The computation of the $O(N^2)$ DFT is accomplished very efficiently using the $O(NlogN)$ FFT algorithm, that comes in canned form from many different libraries.

Thus, if $x[n]$ is your original domain signal, (sampled in space, time, or whatever other quantity), the PSD is given by:

$$ PSD[k] = |\sum_{n=o}^{N-1} x[n] \ e^{\frac{-j \ 2 \pi n k}{N}}|^2 = |X[k]|^2 $$

(Where $X[k]$ is the Discrete Fourier Transform (DFT) of $x[n]$). (Also note, the above formula is for same-length DFTs).

The second aspect of your quandary seems to rest on:

1. What domain to best represent your orignal LIDAR measurements in.
2. How to deal with non-evenly spaced samples.
3. How to deal with missing samples.

Regarding (1): I do not see why you cannot retain the data in its original $r - d \theta$ form. You are after all after radial distance per spatial sample, and this is exactly what LIDARs do for you. You want to then use the DFT as a transformation that will ostensibly magnify and separate different surfaces for you, by virtue of the repetitiveness in the spatial domain, corresponding to particular deltas in the Fourier domain. (Grass with high spatial repetitivity, VS concrete with low repetitivity).

Regarding (2) & (3): Both those problems can easily be solved using linear regression, and/or interpolation, whichever one suits your fancy. That is, you would collect all your $r- d \theta$ data, replete with missing samples and non-even ones. What you seek however is a uniformly spaced sampling on a $d \theta$ grid. You can do a simple polynomial fit via Least-Squares-Estimation, (LSE). That is, your objective will be to decipher the co-efficients of a best-fit polynomial (of degree constrained by the LIDAR physics and spatial bandwidth of the target surface), that best fits your data.

For example, let $d \theta$ be the LIDAR independent variable. You have elected to use a degree-2 polynomial ($D = 2$) to best fit your data in the Least Squared (LS) sense. You have collected a bunch of points from your LIDAR, at some given $d \theta$'s and corresponding radial distances $r$. You would like to find the co-efficients $p_i$, that give you the line of best fit. That is, you would like to find the $p_i$'s in:

$$ r = p_0 + p_1 \ (d \theta) + p_2 \ (d \theta)^2 $$

Thus let $\bf{r}$ be a $1$ x $N$ vector of recorded radial distances to the scanned surface. Let $\boldsymbol{d \Theta}$ be a $(D+1)$ x $N$ data matrix, such that the first row is all $1$s, the second row is composed of each $d \theta$ for each value of $n$, the third row is composed of each $(d \theta)^2$ for each value of $n$, etc. $N$ of course, is the total number of data points from the LIDAR scan of the surface. Finally, let $\bf{p}$ be the $(D+1)$ x $1$ vector of co-efficients to be solved for. Then, the LSE solution is:

$$ \bf{p} = (\boldsymbol{d \Theta} \ \boldsymbol{d \Theta}^T)^{-1} \boldsymbol{d \Theta} \ \bf{r}^T $$

With the coefficient vector $\bf{p}$ in hand, you can now go back and solve for $r$'s at any $d \theta$ grid. That is, simply construct a uniformly spaced $d \theta$ grid, and solve for corresponding $r$'s. Then, you will be ready to transform this regressed result into the Fourier domain, and look for characteristic spatial periodicity corresponding to different materials.