23

You can use logarithms to get rid of the division. For $(x, y)$ in the first quadrant: $$z = \log_2(y)-\log_2(x)\\ \text{atan2}(y, x) = \text{atan}(y/x) = \text{atan}(2^z)$$ Figure 1. Plot of $\text{atan}(2^z)$ You would need to approximate $\text{atan}(2^z)$ in range $-30 < z < 30$ to get your required accuracy of 1E-9. You can take advantage of ...


20

The naive implementation of an $N$-point DFT is basically a multiplication by a $N \times N$ matrix. This results in a complexity of $\mathcal{O}(N^2)$. One of the most common Fast Fourier Transform (FFT) algorithm is the radix-2 Cooley-Tukey Decimation-in-Time FFT algorithm. This is a basic divide and conquer approach. First define the "twiddle factor" ...


20

The answer is the same to the question: "Why do we need computers to process data when we have paper and pencil?" DTFT as well as the continuous-time Fourier Transform is a theoretical tool for infinitely long hypothetical signals. the DFT is to observe the spectrum of actual data that is finite in size.


19

http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629 DFT, size 16 FFT, size 16 The difference in complexity is pretty evident from that, isn't it? Here's how I understand FFT. First off, I would always think about Fourier transforms foremostly as transforms of continuous functions, i.e. a bijective mapping $\operatorname{FT} ...


15

Here is a picture to add to Robert's good answer demonstrating the "re-use" of operations, in this case for an 8 point DFT. The "Twiddle Factors" are represented in the diagram using the notation $W_N^{nk}$ which is equal to $e^{j2\pi \frac{nk}{N}}$ Note the path shown and the equation underneath shows the result for the frequency bin X(1), as given by ...


14

This is a well-studied problem, dating back from the mid 90s (DARPA/NIST broadcast transcription challenges). Search for "speech/music segmentation" or "audio segmentation" and you'll find thousands of research papers. There are two broad approaches to solve this problem: Supervised classification Train a speech/music classifier, using a standard machine ...


13

if you want a cheap and dirty optimized power-series expansion (the coefficients for Taylor series converge slowly) for sqrt() and a bunch of other trancendentals, i have some code from long ago. i used to sell this code, but no one has paid me for it for nearly a decade. so i think i'll release it for public consumption. this particular file was for an ...


12

TL, DR: world pervasive algorithms (FFT-related)! The continuous Fourier transform, the Discrete-time Fourier transform (DTFT) and the Discrete Fourier transform (DFT) share conceptually similar traits (regarding energy, convolution, shift, scale, etc.) The DFT unveiled a very scalable and fast algorithm to put those concepts into practice: the FFT. It is ...


11

Thanks to the plot in Olli Niemitalo's answer I got convinced that the formula given in the book has a sign error. The non-linearity used for fuzz or distortion is always some type of smoothed clipping function, which compresses the input signal. So small input amplitudes experience little change whereas high input amplitudes are (more or less) softly ...


11

The basic technique to place a mono source in stereo is called constant power panning. If you want to place a mono source at angle $\theta$ you can just use $A_\mathrm{amp}$ and $B_\mathrm{amp}$ as amplitudes for your channels: $A_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} + \sin{\theta})$ $B_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} - \sin{\...


9

the general polynomial form is: $$\begin{align} f(u) &= \sum\limits_{n=0}^{N} \ a_n \ u^n \\ \\ &= a_{\small{0}} + \Bigg(a_{\small{1}} + \bigg(a_{\small{2}} + \Big(a_{\small{3}} + \,... \big(a_{\small{N-2}} + (a_{\small{N-1}} + a_{\small{N}} \,u \,)u \, \big)u \ ...\Big)u \, \bigg)u \, \Bigg)u\\ \end{align}$$ the latter form is using Horner's ...


9

PROLOGUE My answer to this question is in two parts since it is so long and there is a natural cleavage. This answer can be seen as the main body and the other answer as appendices. Consider it a rough draft for a future blog article. Answer 1 * Prologue (you are here) * Latest Results * Source code listing * Mathematical justification for ...


9

You mention in a comment that your target platform is a custom IC. That makes the optimization very different from trying to optimize for an already existing CPU. On a custom IC (and to a lesser extent, on an FPGA), we can take full advantage of the simplicity of bitwise operations. In addition, to reduce the area it is not only important to reduce the ...


8

essentially, in computing the naive DFT directly from the summation: $$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{j 2 \pi \frac{nk}{N}} $$ there are $N$ table lookups for the twiddle factor $ e^{j 2 \pi \frac{nk}{N}} $, $N$ complex multiplications, and $N-1$ additions. and that's just for one value of $X[k]$ and one instance of $k$. then the naive DFT ...


8

Question 1 If you apply it over the entire length of the array, the length of the FFT would be the length of the array. But, the FFT is more efficient if the length is a power of two, so it is common to pad 0's onto the end of the signal until its length is a power of 2. Overly simple example... x = [3.4, 2.56, 1.3] x has a length of 3, the next power of ...


8

1. Logarithms and exponents to avoid multiplication To completely avoid multiplication, you could use $\log$ and $\exp$ tables and calculate: $$I^2 + Q^2 = \exp\!\big(2\log(I)\big) + \exp\!\big(2\log(Q)\big).\tag{1}$$ Because $\log(0) = -\infty,$ you'd need to check for and calculate such special cases separately. If for some reason accessing the $\exp$ ...


8

Given two complex numbers $z_1=a_1+jb_1$ and $z_2=a_2+jb_2$ you want to check the validity of $$a_1^2+b_1^2>a_2^2+b_2^2\tag{1}$$ This is equivalent to $$(a_1+a_2)(a_1-a_2)+(b_1+b_2)(b_1-b_2)>0\tag{2}$$ I've also seen this inequality in Cedron Dawg's answer but I'm not sure how it is used in his code. Note that the validity of $(2)$ can be checked ...


7

The beat and onset detection algorithms used at the Echo Nest are probably variants/improvements of the techniques developed by Tristan Jehan in his Ph.D. This is not the only approach, and I would recommend you to try first: Getting an onset detection function using spectral flux or complex amplitude. Using this algorithm to detect beats (you can improve ...


7

If you haven't seen it, the "Quake square root" is simply mystifying. It uses some bit-level magic to give you a very good first approximation, and then uses a round or two of Newton's approximation to revise. It might help you if you're working with limited resources. https://en.wikipedia.org/wiki/Fast_inverse_square_root http://betterexplained.com/...


6

You might try Vincent Falco's Collection of Useful C++ Classes for Digital Signal Processing. The StackOverflow Question A C++ library for IIR filter offers a few more suggestions. Finally: the hard part for Butterworth/Chebychev filters is really the design, not the implementation. You might consider doing the design in Matlab (or Gnu Octave), or an ...


6

Unfortunately you did not mention what your channel is so I assume that merely AWGN is present. Furthermore I assume symbol-by-symbol detection is desired since a decoder for the Hammingcode is following later. So a single symbol $x$ is received as $y = x + n$ where $x \in [1,j,-j,-1]$ and $n$ is distributed according to a zero-mean complex normal ...


6

because the code markup for SE seems to work like shit, i'll try to answer this more directly, specifically for the $\sqrt{x}$ function. yes, a power series can quickly and efficiently approximate the square root function, and only over a limited domain. the wider the domain, the more terms you will need in your power series to keep the error sufficiently ...


6

You could also approximate the square root function by using Newton's Method. Newton's Method is a way of approximating where the roots of a function are. It is also an iterative method where the result from the previous iteration is used in the next iteration until convergence. The equation for Newton's method to guess where the root is of a function $f(...


6

Note: I originally posted this answer for the Stack Overflow copy of this question, before realizing that it had also been asked here. It somewhat duplicates pichenettes' answer, but I felt it still worth (re)posting here, since it includes some extra details. (Whether those details are useful or not, I'll leave for you and the OP to judge.) If you know ...


6

You want a method that removes noise while preserving edges. This cannot be achieved well by linear filtering, as you noticed yourself. I know of two approaches that might work well for your problem. The first is median filtering, where samples inside a window are replaced by their median. The following plot shows the result of median filtering with a window ...


6

I just wanted to point out that if you're planning to use these formulas in your code, you can get the exact same results with fewer calculations by using an angle $\theta$ between 0 and 90 degrees and simply calculating $A_{amp} = \sin(\theta)$ and $B_{amp} = \cos(\theta)$. You may have run across these before (they seem to be more commonly referenced in ...


6

There's a whole area of signal processing dedicated to optimal filtering. In pretty much every case I've seen the filtering problem is formulated with a convex cost function. Here's a freely available book on the subject - Sophocles J. Orfanidis - Optimum Signal Processing.


6

You need to build a time varying delay, where you can modulate the delay amount over time. The peak delay modulation is a function of your maximum desired frequency shift and the modulation frequency. This is not trivial since it will require fractional sample delays with some kind of interpolation algorithm. You can't round to the nearest integer delay ...


6

Here I expected $y(n)$ is to be computed by convolving $x(n)$ with $h(n)$, but in the equation given by Wikipedia it is shown as a matrix multiplication $y(n) = h^H(n).x(n)$. Are these two operations(convolution and matrix multiplication) same here?. The system is an FIR system, so the vector multiplication here is equivalent to convolution --- for ...


6

Not sure if this has a name, but it is a nonlinear low pass filter that uses different smoothing constants depending on the input signal deviation from the filtered output. Small deviations are typically assumed to indicate consistency with the smooth estimation and result in little adaption to the input, while large deviations indicate a relevant state ...


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