Thread: A question about handling Fourier transform results

Not gonna lie, that post looked so intense I was like, "Dude...."

Like, seriously, that's intense. What'd you use to make all those plots? Like MatLab or something?

Gnuplot.

Time-domain versus frequency-domain, and sampling issues, are quite important, if the correctness of your results matter. It is quite closely related to the concepts Dr. Alvy Ray Smith pointed out in his 1995 article A Pixel is Not a Little Square.

In this thread the issue is that the discrete Fourier transform can be used for more than one purpose. Usually, it is used to exactly define the frequency characteristics of a periodic signal. However, it can also be used to sample the frequency profile of a nonperiodic signal -- a wave packet --; the Nyquist-Shannon sampling theorem tells us that that sampling is sufficient (to describe the original signal).

My post proves that. I used three different wave packets, with wildly different frequency-domain profiles, to show that it does not work for just some signals, it works for all kinds of signals.

I felt that Brewbucks response sounded so convincing at the face of it, that I'd have to show absolute proof to convince anyone that he was utterly wrong.

It was important to show he was wrong, because DFT is the correct method to analyse the frequency characteristics of a discrete time-domain nonperiodic signals. If someone believes his post, they'll wander off elsewhere, looking for a solution in the wrong places.

---

To generate the pulse images (pulse-curve.svg and pulse-dft), I used

Code:

./dft $[7*5000] 7 pulse > pulse.out
gnuplot -e 'set term svg enhanced size 800,480
            set output "pulse-curve.svg"
            set title "Pulse"
            unset xtics
            unset ytics
            set xrange [0:1.0/7.0]
            set yrange[-1.1:1.1]
            plot "pulse.out" u 1:2 notitle w lines'

gnuplot -e 'set term svg enhanced size 800,480
            set output "pulse-dft.svg"
            set title "Detail of pulse DFT"
            unset xtics
            unset ytics
            set xrange [0:0.003]
            set yrange [-1300:1300]
            plot "pulse.out" u 1:3 t "|F7|" w lines lc 1, \
                 "pulse.out" u 1:4 t "Re F7" w lines lc 2, \
                 "pulse.out" u 1:5 t "Im F7" w lines lc 3, \
                 "pulse.out" u 1:6 t "|F|" w points pt 6 lc 1 \
                 "pulse.out" u 1:7 t "Re F" w points pt 2 lc 2 \
                 "pulse.out" u 1:8 t "Im F" w points pt 1 lc 3'

For the sawpulse-curve.svg and sawpulse-dft.svg I used

Code:

./dft $[3*5000] 3 sawpulse > sawpulse.out
gnuplot -e 'set term svg enhanced size 800,480
            set output "sawpulse-curve.svg"
            set title "Sawpulse"
            unset xtics
            unset ytics
            set xrange [0:1.0/3.0]
            set yrange[-1.1:1.1]
            plot "sawpulse.out" u 1:2 notitle w lines'
gnuplot -e 'set term svg enhanced size 800,480
            set output "sawpulse.svg"
            set title "Detail of sawpulse DFT"
            unset xtics
            unset ytics
            set xrange [0.00237:0.0107]
            set yrange [-800:900 ]
            plot "sawpulse.out" u 1:3 t "|F3|" w lines lc 1, \
                 "sawpulse.out" u 1:4 t "Re F3" w lines lc 2, \
                 "sawpulse.out" u 1:5 t "Im F3" w lines lc 3, \
                 "sawpulse.out" u 1:6 t "|F|" w points pt 6 lc 1, \
                 "sawpulse.out" u 1:7 t "Re F" w points pt 2 lc 2, \
                 "sawpulse.out" u 1:8 t "Im F" w points pt 1 lc 3'

I picked 5000 samples for the length of the original wave above, as that was the maximum sample window length for the OP in this thread. So, just to make sure the graphs were representative of the issues at hand. Multiple birds with one stone, and so on.

Finally, for the noise, I used

Code:

./dft $[1024*1024] 1024 noise > noise.out
gnuplot -e 'set term svg enhanced size 800,480
            set output "noise-curve.svg"
            set title "Noise"
            unset xtics
            unset ytics
            set xrange [0:1.0/1024.0]
            set yrange[-1.1:1.1]
            plot "noise.out" u 1:2 notitle w lines'
gnuplot -e 'set term svg enhanced size 800,480
            set output "noise-dft.svg"
            set title "Detail of noise DFT"
            unset xtics
            unset ytics
            set xrange [0:32.0/1024.0]
            set yrange [-40:50]
            plot "noise.out" u 1:3 t "|F1024|" w lines lc 1, \
                 "noise.out" u 1:4 t "Re F1024" w lines lc 2, \
                 "noise.out" u 1:5 t "Im F1024" w lines lc 3, \
                 "noise.out" u 1:6 t "|F|" w points pt 6 lc 1, \
                 "noise.out" u 1:7 t "Re F" w points pt 2 lc 2, \
                 "noise.out" u 1:8 t "Im F" w points pt 1 lc 3'

which is at the interesting limit -- squaring the number of samples. (I could have used 5000*5000 here, but without wisdom, that DFT takes over a minute -- didn't check, just aborted the calculation at one minute point -- whereas the 1024*1024 one is just a few seconds, on this laptop. Even the 1024*1024 one has way too much data; 128*128 would have been just as illustrative.)

As you can see, only the file names and x and y ranges vary. I removed the x and y axis labels and tics, because the data is represented normalized, and anyway the axis values are not important; the fact that the samples match the higher-precision, zero-padded DFTs exactly, is what matters.

Using ImageMagick and NetPBM,

Code:

convert -background white input.svg output.pnm
pnmquant 32 output.pnm | pnmtopng -compress 9 output.png
rm -f output.pnm

produces nice, antialiased 32-color PNG images out of the SVG images.

Of course, I didn't do the above by hand. I used make using the following Makefile:

Code:

CC	:= gcc
CFLAGS	:= -Wall -O2 -fomit-frame-pointer
LDFLAGS	:= -lm -lfftw3
AWK	:= awk
GNUPLOT	:= gnuplot
CONVERT	:= convert
PNMQUANT := pnmquant 2>/dev/null
PNMTOPNG := pnmtopng -compress 9 2>/dev/null

.PHONY: all clean run

all: clean run

clean:
	rm -f dft
	rm -f pulse.out pulse-input.png pulse-dft.png
	rm -f sawpulse.out sawpulse-input.png sawpulse-dft.png
	rm -f noise.out noise-input.png noise-dft.png

dft: dft.c
	$(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@

run-pulse: dft
	@echo
	./dft 35000 7 pulse > pulse.out
	@$(AWK) 'NR > 1 && NF >= 6 && $$3 != $$6' pulse.out
	@echo

	@rm -f pulse.svg pulse.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "pulse.svg" ; set title "Pulse" ; unset xtics ; unset ytics; set xrange [0:1.0/7.0] ; set yrange[-1.1:1.1] ; plot "pulse.out" u 1:2 notitle w lines'
	@$(CONVERT) -background white pulse.svg pulse.pnm
	@$(PNMQUANT) 32 pulse.pnm | $(PNMTOPNG) -compress 9 > pulse-input.png
	@rm -f pulse.svg pulse.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "pulse.svg" ; set title "Detail of pulse DFT" ; unset xtics ; unset ytics; set xrange [0:0.003] ; set yrange [-1300:1300]; plot "pulse.out" u 1:3 t "|F7|" w lines lc 1, "" u 1:4 t "Re F7" w lines lc 2,  "" u 1:5 t "Im F7" w lines lc 3, "" u 1:6 t "|F|" w points pt 6 lc 1, "" u 1:7 t "Re F" w points pt 2 lc 2, "" u 1:8 t "Im F" w points pt 1 lc 3'
	@$(CONVERT) -background white pulse.svg pulse.pnm
	@$(PNMQUANT) 32 pulse.pnm | $(PNMTOPNG) -compress 9 > pulse-dft.png
	@rm -f pulse.svg pulse.pnm

run-sawpulse: dft
	@echo
	./dft 15000 3 sawpulse > sawpulse.out
	@$(AWK) 'NR > 1 && NF >= 6 && $$3 != $$6' sawpulse.out
	@echo

	@rm -f sawpulse.svg sawpulse.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "sawpulse.svg" ; set title "Sawpulse" ; unset xtics ; unset ytics; set xrange [0:1.0/3.0] ; set yrange[-1.1:1.1] ; plot "sawpulse.out" u 1:2 notitle w lines'
	@$(CONVERT) -background white sawpulse.svg sawpulse.pnm
	@$(PNMQUANT) 32 sawpulse.pnm | $(PNMTOPNG) -compress 9 > sawpulse-input.png
	@rm -f sawpulse.svg sawpulse.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "sawpulse.svg" ; set title "Detail of sawpulse DFT" ; unset xtics ; unset ytics; set xrange [0.00237:0.0107] ; set yrange [-800:900 ] ; plot "sawpulse.out" u 1:3 t "|F3|" w lines lc 1, "" u 1:4 t "Re F3" w lines lc 2,  "" u 1:5 t "Im F3" w lines lc 3, "" u 1:6 t "|F|" w points pt 6 lc 1, "" u 1:7 t "Re F" w points pt 2 lc 2, "" u 1:8 t "Im F" w points pt 1 lc 3'
	@$(CONVERT) -background white sawpulse.svg sawpulse.pnm
	@$(PNMQUANT) 32 sawpulse.pnm | $(PNMTOPNG) -compress 9 > sawpulse-dft.png
	@rm -f sawpulse.svg sawpulse.pnm

run-noise: dft
	@echo
	./dft 1048576 1024 noise > noise.out
	@echo

	@rm -f noise.svg noise.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "noise.svg" ; set title "Noise" ; unset xtics ; unset ytics; set xrange [0:1.0/1024.0] ; set yrange[-1.1:1.1] ; plot "noise.out" u 1:2 notitle w lines'
	@$(CONVERT) -background white noise.svg noise.pnm
	@$(PNMQUANT) 32 noise.pnm | $(PNMTOPNG) -compress 9 > noise-input.png
	@rm -f noise.svg noise.pnm
	$(GNUPLOT) -e 'set term svg enhanced size 800,480 ; set output "noise.svg" ; set title "Detail of noise DFT" ; unset xtics ; unset ytics; set xrange [0:32.0/1024.0] ; set yrange [-40:50] ; plot "noise.out" u 1:3 t "|F1024|" w lines lc 1, "" u 1:4 t "Re F1024" w lines lc 2,  "" u 1:5 t "Im F1024" w lines lc 3, "" u 1:6 t "|F|" w points pt 6 lc 1, "" u 1:7 t "Re F" w points pt 2 lc 2, "" u 1:8 t "Im F" w points pt 1 lc 3'
	@$(CONVERT) -background white noise.svg noise.pnm
	@$(PNMQUANT) 32 noise.pnm | $(PNMTOPNG) -compress 9 > noise-dft.png
	@rm -f noise.svg noise.pnm

run: run-pulse run-sawpulse run-noise

The dft.c generator program is written so that it's very easy to add new waveforms -- just the function itself, a new if statement to select that function at run-time, and a new print line to name that function if the user runs the program without parameters.

I didn't write all that just to prove one point; I hoped the program would help others explore the properties of DFTs, comparing the DFT of a periodic signal, to the DFT of that same signal padded with zeroes, "simulating" a nonperiodic signal.

---

One subproblem (of MD simulation sonification) I'm working on right now involves one example method of generating noise whose spectrum describes the statistical distribution of some discrete variables -- in particular, the kinetic, potential, or total energies of particles in a molecular dynamic simulation.

It is not a difficult problem at all, but to be useful, one has to take into account the logarithmic sensitivity of human hearing, the psychoacoustic model (hearing sensitivity is highly dependent on the frequency, so you'll want to start with that profile as your "white noise"), and acoustic glitches: I'm not generating just one profile, but trying to generate continous sound based on a sequence ("movie") of statistical snapshots taken as a simulation progressed. Plus, the snapshots are not guaranteed to be taken at constant intervals. And I'd like to synch it to a visual animation, too.

The purpose of this is that in thermodynamic equilibrium, that is the Maxwell-Boltzmann distribution, and the human ear is much more sensitive to qualitative changes in it, than comparing two curves by eye. (I have experimented on just using the difference of the distributions, but the problem there is the lack of anything to compare to: yes, you can tell the differences, but it's very hard to tell whether those differences matter, or are interesting.)

Simply put, the sound generated this way is a pretty good indicator whether the system is in thermodynamic equilibrium, or if there are some interesting processes ongoing, and perhaps even what kinds those processes might be (based on their velocity characteristics compared to the overall temperature of the system).

Originally Posted by Mario F.

Also, since you are trying this out based on sound modulation, you must be aware we are particularly keen to sound illusion (see shepard's tones), and other false representations when sounds aren't distinct by more than one or two variables.

Yup. Artefacts on modulation are so bad (very audible!) I'm actually looking at doing this via filtering ("noise shaping"), probably using quite a few separate filters successively.

(I could tell you a tale how I tried to use a windowed inverse Hartley transform for generating such sound, and it kinda works (handwaving here), but the artefacts from the windowing function are extremely audible, and very distracting. Spent quite a bit of time on it, too, with no good results.)

You're absolutely correct that this specific example of sonification is of very limited use. My thought is that it is useful when first observing a simulation run, say when you first open some molecular simulation in a visualization program like ovito. Shaped noise might be useful in describing the overall kinetic state of the system in a qualitative sense; an intuitive assist for locating the interesting/problematic points in a simulation. It is also just one example method of several I'm developing.

Originally Posted by Mario F.

I would still choose any visual representation to detect changes over time.

I'm trying to approach sonification and visualization as complementary methods, not alternatives. Information density is a real problem, so I'm hoping that using more than one sense, I can convey a better "picture" of what is happening.

Anyway, to return to the original topic in this thread:

The simulations I use to check my ideas include both elastic (material stressed but returns to original state) and plastic deformation (material "breaks", does not return to original state). The boundary between the two is very interesting, and observable at the atomic level. You get defects that start to cluster together and won't go away when the strain that caused them goes away. (It would be very interesting to be able to sonify these, but the difficulty is recognizing them based on the numerical data -- at room temperatures, nothing is static in an atomic simulation.) Compared to the rest of the material, these defects move fast. Any sudden changes -- thus higher frequencies in the frequency domain -- are likely to be related to these defects. Because there are both elastic and plastic changes, comparing the higher frequency components of the measurements across events, is likely to tell you how the changes differ across events -- i.e., how the material evolves, and suffers from fatigue.

As a computational physicist, one must accept that we use imperfect computational methods. To simulate large systems, we have to use approximate potentials (classical potential models) instead of say electron density calculations. To a numerical problem, there is no perfect answer. (Actually, the first thing you have to do with any numerical result is analyse whether it makes any practical sense or not.)
Approximations and mathematical tools like Nyquist-Shannon sampling theorem are indispensable.

This is why I felt so strongly about Brewbucks dismissal of DFT wrt. nonperiodic signals. I felt that it was like someone claiming you couldn't weld with explosives; that explosives only blow stuff apart. I thought it was important to refute that, in a manner which is plausible to any programmer, without any background in this stuff.