66

The main challenge to make something general is to create a list of points that, when joined, is as close as possible to the original image. The first example (elephant) focuses on the mathematical idea. In the second example (Bart Simpson), the code is wrapped in a function. The last example gives some ideas on how to apply the method to a photograph ...


61

It looks like you need some basic facts on numerical Fourier transforms. I am going to assume that you are going from the time domain to the frequency domain. The number of points in the time domain equals the number of points in the frequency domain. As the data is sampled in the time domain i.e. is a set of equally spaced points, then in the frequency ...


46

Here's what I did back in pre-MMA6 times. I believe you can easily adapt this method to work with all the fancy commands that have been added since then. Someone proficient in CDF programming could also make the mask building process interactive with disks and polygons: one window shows the FFT, the user places shapes to hide unwanted features, and another ...


45

I can explain this. The definite flavor of Integrate works with assumptions in a few ways. One is to use them in Simplify, Refine, and a few other places that accept Assumptions, to do what they will in the hope of attaining a better result (it also uses them to determine convergence and presence or absence of path singularities). Those places also get the $...


33

The Fourier transform is defined as: $$ \begin{align*} H(f) &= \int_{-\infty}^\infty h(t) \mathrm e^{2\pi \mathrm i f t}\mathrm dt\\ h(t) &= \int_{-\infty}^\infty H(f) \mathrm e^{-2\pi \mathrm i f t}\mathrm df \end{align*} $$ where $ h(t) $ is the signal, and $ H(f) $ is its Fourier transform; if $ t $ is measured in seconds, then $ f $ is measured ...


32

img = Import["ExampleData/lena.tif"]; Image[img, ImageSize -> 300] data = ImageData[img];(*get data*) {nRow, nCol, nChannel} = Dimensions[data]; d = data[[All, All, 2]]; d = d*(-1)^Table[i + j, {i, nRow}, {j, nCol}]; fw = Fourier[d, FourierParameters -> {1, 1}]; (*adjust for better viewing as needed*) fudgeFactor = 100; abs = fudgeFactor*Log[1 + Abs@...


29

This is how I would do this. Define frequencies and sampling rate precisely. Then use Periodogram because it takes SampleRate as an option and rescales frequency axis automatically. Read up Docs on Periodogram - see examples there. data = Table[{t, Sin[2 Pi 697 t] + Sin[2 Pi 1209 t]}, {t, 0., 0.1, 1/8000.}]; ListLinePlot[data, AspectRatio -> 1/4, ...


26

I think there are at least three elements to consider here: FourierTransform and Fourier, by default, output results in different forms Plotting Sin[x] UnitStep[x] is not the same as Sin[x] and behaves differently when used in conjunction with Fourier and FourierTransform Plot does not handle DiracDelta elegantly The signal processing form of the Fourier ...


21

I appreciate very much, that you wrote up such a nicely formatted question, although this is your first post. Therefore, let's put the comments into an answer. Your first issue was that you used ( ) where you should have used [ ]. That's maybe not obvious for starters and I have seen this mistake very often. There are different types of braces and it's ...


20

It always takes me a while to remember the best way to do a numerical Fourier transform in Mathematica (and I can't begin to figure out how to do that one analytically). So I like to first do a simple pulse so I can figure it out. I know the Fourier transform of a Gaussian pulse is a Gaussian, so pulse[t_] := Exp[-t^2] Cos[50 t] Now I set the timestep ...


19

Based on the MATLAB documentation, it would appear that this is accomplished by simple zero-filling. As such, you can obtain the same result in Mathematica using Fourier[PadRight[list, n, 0.], FourierParameters -> {1, -1}] where list is your signal and n is the desired length. For a multidimensional FFT, replace n with {n1, n2, ...}, where n1, n2, &...


18

update Just to clean things up a bit, we can use the discussion here to make a couple functions that help extract the frequency data from this dataset. I define two functions findPeriod and reconstruct: Clear[findPeriod]; findPeriod[data_, threshold_] := Module[{fs, s1, s = {}, i, a0f, af, pf, pos, fr, frpos, fdata, fdatac, n, per}, n = Length[data];...


18

Physicist chiming in - Hi!. I believe there has been some confusion here. It seems to me that OP is meaning to plot an Airy disk which was studied by G.B. Airy but is not given by the Airy function. It is given by the Fourier transform of the indicator function of the unit circle, which actually happens to be a Bessel function (see e.g. wikipedia). If I ...


17

I think perhaps you need codes like this: Func[x_] := Sin[x]; tmin = 0; tmax = 10; \[CapitalDelta]t = (tmax - tmin)/100; tgrid = Table[t, {t, tmin, tmax, \[CapitalDelta]t}]; wgrid = RotateRight[(2 \[Pi])/(tmax - tmin)* Range[-((Length@tgrid - 1)/2), (Length@tgrid - 1)/2], ( Length@tgrid - 1)/2]; ListLogLogPlot[{wgrid, (tmax - tmin)/Sqrt[2 \[Pi]*...


17

CUDALink allows you to write custom kernels but unfortunately it doesn't allow to use its CUDAMemory in other functions. I found two methods to deal with it. Documented method You can use LibraryLink and call CUDA functions as in a regular CUDA program. For example, you can write something like fourier.cu. Then you can compile it and load "fourier_single" ...


16

As mentioned by @Rahul, you have not sampled your sine wave often enough and have introduced artifacts due to aliasing. The frequency of Sin[500 x]=Sin[2 Pi f x] is $f=500/(2\pi)$, which is about 80 Hz. At least two samples per cycle are required to avoid aliasing, hence the default $x$ interval of 1 in {x,0,100} must be reduced to less than about $1/(2*80)=...


16

Finally I found the most promising algorithm proposed in this really good reference Manuel Guizar-Sicairos and Julio C. Gutiérrez-Vega, "Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields," J. Opt. Soc. Am. A 21, 53-58 (2004). The authors call the algorithm pth-order quasi-discrete Hankel Transform (pQDHT) ...


16

Looking at your plotted data you can see about 40 cycles of the dominant frequency, this tells you that the peak will appear somewhere around the 40th element of the DFT. That's in the region where your plot of the DFT is clipped, so it's no wonder you can't see the peak. Looking at the relevant part of the DFT you can see the peak quite clearly: ...


16

I will take the data as a time history. First I assume that the x-values are equally spaced and work out the time increment and frequency increment and then plot the data. tinc = data[[2, 1]] - data[[1, 1]]; finc = 1/(tinc Length[data]); ListPlot[data] This looks like almost two cycles of a sine wave with a frequency of about 0.1 and an amplitude of 17. ...


16

I like @Vitaliy's answer, but here's another approach using Fourierinstead of Periodogram. time = 2; tinc = 0.001; sampls = Table[Sin[n*(2 Pi) 4], {n, 0, 2, tinc}]; nyquist = 1/(2 tinc) len = Length@sampls; ListLinePlot[Sqrt[4/len] Abs@Fourier[sampls], PlotRange -> {{0, 10}, All}, DataRange -> {0, (len - 1)/time}] Briefly, I construct a sample set ...


16

First, I have a few improvement suggestions for your Fourier code: The bright vertical and horizontal lines you see in your Fourier image are the sharp gradients at the borders of the image (because the Fourier transform assumes a periodic image). So you should get rid of the black border at the bottom: img = Import["http://i.stack.imgur.com/bIUkE.png"]; ...


16

Note that the expression returned by Sum is correct and equals $x(1-x)$ for $0 \leq x \leq 1$. I assume your question is how to simplify the expression into $x(1-x)$? I was able to hack a solution, and unfortunately I don't think it scales very well to other expressions. But here goes: First, evaluate the sum: sum = 1/6 - Sum[Cos[2 x Pi n]/(Pi n)^2, {n, 1, ...


16

First, import the audio and extract usable data from it: audio = Sound[ SampledSoundList[ Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"], 22050]] audioDuration = Duration[audio]; audioSampleRate = AudioSampleRate[audio]; data = AudioData[audio][[1]]; Second, use PeakDetect to see which points are peaks (= 1) and which points are ...


14

Here is an explicit way to calculate the frequency corresponding to each element of the output of the Fourier command. The frequencies will depend on two values: the sampling interval and the number of points in the data analysis. ssf = RotateRight[Range[-n/2, n/2 - 1]/(n sampInt), n/2]; where n is the number of points analyzed and sampInt is the time ...


14

Since you are using this as an excuse to learn to code in Mathematica, I'll try to help with that in mind. Sometimes I find it nice to design "from the top down", document from the top down, split in many small functions with no state, and go testing them "from the bottom up". As you gain confidence, you will perhaps use coarser functions, and not test every ...


14

For the reasons mentioned above, I wrote the following "shell" for these functions: (* Weak test for singularity *) findSingular[coe_, k_, assump___] := Piecewise[{{{}, # === k}}, #] &[ k /. Union@Solve[{assump, Denominator@Together@coe == 0, k ∈ Integers}, k]] generalTerm[expr_, {t_, a_, b_}, n_: C][coefunc_, termfunc_, assump___] := With[{...


13

You could simply remove the vertical waves (e.g. by subtracting the median of each column) and histogram modification. Using @bill s's cropped image: img = Image[ ImageData[ ColorConvert[Import["http://i.stack.imgur.com/EvjuW.png"], "Grayscale"]][[;; , ;; , 1]]]; (* remove alpha channel *) columnMedian = Median[ImageData[img]]; medianRemoved = # ...


13

I also would expect Mathematica to simplify all Fourier transformed derivatives equally, but it may be understandable that the simplifications are harder to see when the derivative is not taken with respect to the innermost Fourier transform variable. To work around this problem, you could change the order of integrations for the Fourier transform to ...


13

Taking the Fourier transform is easy and fun! Let's strip away some of the complexities. First, remove the color from the image, since this just complicates things (you can always take the transform of each color channel separately). img = ColorConvert[Import["ExampleData/lena.tif"], "Grayscale"]; Here are the magnitude and phase of the Fourier transform: ...


Only top voted, non community-wiki answers of a minimum length are eligible