# Tag Info

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 ...

24

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, ...

20

The Fourier transform is defined as: $$H(f)=\int_{-\infty}^\infty h(t) e^{2\pi i f t}dt\\ h(t)=\int_{-\infty}^\infty H(f) e^{-2\pi i f t}df$$ where $h(t)$ is the signal, and $H(f)$ is it's Fourier transform, if $t$ is meassured in second, then $f$ is measured in Hz. The discrete Fourier transform is defined as: $$... 18 The problem here is that Mathematica doesn't recognize {x, y, z} as some kind of a vector object that should be treated as grouped together; instead, it substitutes in three independent variables, and probably starts integrating them one by one. The result is a very complicated integral. If you do the coordinate transformation yourself, you can reproduce ... 18 Here's a possible starting point for a solution. It splits the sample list into chunks and measures the Norm of the sample Differences in each chunk, and then does the FFT on that data. bpmplot[snd_, bpmmax_: 300] := Module[{samples, minfreq, signal, fft}, samples = snd[[1, 1, 1]]; minfreq = snd[[1, 2]]/Length[samples]; signal = (Norm[Differences[#]]) ... 18 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 ... 17 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, ... 17 There is the function NFourierTransform[] (as well as NInverseFourierTransform[]) implemented in the package FourierSeries`. The function, as with the related kernel functions, takes a FourierParameters option so you can adjust computations to your preferred normalization as needed. For your specific normalization, you apparently want the setting ... 17 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 ... 16 So, you have a function F(x,y) = f_x(x)g_y(y) + g_x(x)f_y(y), and you want to recover f_x,g_y,g_x,f_y. If you've tabulated the values of F(x,y) in a matrix \mathbf F with entries f_{ij} = F(x_i,y_j), then this amounts to decomposing the matrix as$$\mathbf F \approx \mathbf f_x\mathbf g_y^T + \mathbf g_x\mathbf f_y^T,$$where \mathbf f_x,\mathbf ... 15 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 ... 14 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 + ... 14 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. ... 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 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"]; ... 13 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 ... 13 It's not a bug, it's a feature Exact integration returns 1/Sqrt[2 Pi] Integrate[(1 + Sqrt[x])^2 Exp[I k x], {x, -Infinity, Infinity}, Assumptions -> {k \[Element] Reals}] Integrate::idiv: "Integral of E^(I\k\x)\ (1+[Sqrt]x)^2 does not converge on {-Infinity,Infinity}." However we can multiply by Exp[-b Abs[x]] and then put b -> 0 ... 13 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 = ... 13 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: ... 13 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 ... 12 Fourier uses FFT when possible 12 Fourier[list] computes the discrete Fourier transform of list. I assume it uses the FFT when it can. 12 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 = # ... 12 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) ... 12 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 ... 12 I believe the frequency mismatch arises because the endpoints of your 200 point series are offset. The first point has amplitude 0.1, the last 1.5584. As others mention, the Fourier transform assumes periodicity. So the signal you are transforming has a sine component, a cosine component, and a step function offset of the first and last points. The Fourier ... 11 In Mathematica there is a designated function for this, UnitBox, PiecewiseExpand[UnitBox[x]] which gives expected result without assumptions: FourierTransform[UnitBox[x/a], x, k, FourierParameters -> {1, -1}] Abs[a] Sinc[(a k)/2] There is actually a set of designated functions: FourierTransform[UnitTriangle[x/a], x, t, FourierParameters ... 11 If you take a look at the documentation, Mathematica's symbolic Fourier transform function, FourierTransform, computes$$\hat f(k) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{ikx}\mathrm{d}x You can discretize some piece of this integral by limiting $x$ and $k$ to values $x_1 + (r-1)\Delta x$ and $(s-1)\Delta k$ respectively, where \$\Delta ...

11

In the definition of s you're summing from k==0. Since the summand has a term 1/k this gives a divide-by-zero error when calculating the partial sums. The sum should in fact start from k==1 (the zeroth coefficient is taken care of by the constant term in front of the sum). The first few approximations then look like s[n_, x_] := 8/4 + 3/(9 \[Pi]) Sum[(6 ...

11

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 ...

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