Understanding FFTs for simple Sin / Cos this is my first question in this forum, and although I read several threads on this side and googled a lot I could not find the answer to my question (maybe it is too basic)?

For anyone reading this thread with a similar question: The code is now adapted in a way that the answer is partly given in this question. For the explanations look in the answers below.

I found this one, but it only covered my answer partly:
Understanding the FFT phase spectrum with a simple example
If you have good literature that covers this topic, I'm grateful if you can share it.

So, let's assume I have a very basic signal in the time domain, for example:
5*sin(3*2*pi*x-2) with amplitude=5, frequency=3 and shift of 2.

import matplotlib.pyplot as plt
import numpy as np
npts  = 100
tmax  = 10
dt = tmax / npts[![enter image description here][1]][1]
fs = 1 / dt
t     = np.linspace(0, tmax-dt, npts)
y     = 5 * np.sin(3 * 2*np.pi * t - 2)
N = y.size

freq = np.fft.fftshift(np.fft.fftfreq(N)*fs)
yfft = np.fft.fft(y)
yfft = np.where(abs(yfft) < 1.0e-10, 0, yfft)
magn = np.fft.fftshift(np.abs(yfft*dt)/(N*dt))
phase = np.fft.fftshift(np.angle(yfft))
psd = np.fft.fftshift(np.abs(yfft*dt)**2/(N*dt))
yfft = np.fft.fftshift(np.fft.fft(y))

How do I get the frequency 3 and the shift of 2 to show?
Is that even possible?

fig, ax = plt.subplots(3, 2)
fig.tight_layout()
ax[0,0].plot(t, y)
ax[0,0].set_title('Time-Domain')
ax[0,1].plot(freq, yfft)
ax[0,1].set_title('Frequency-Domain')
ax[1,0].plot(freq, magn)
ax[1,0].set_title('Magnitude')
ax[1,1].plot(freq, phase)
ax[1,1].set_title('Phase')
ax[2,0].plot(freq, psd)
ax[2,0].set_title('PSD')

Best regards
Jonas

Creating the Frequency Vector
The arrangement of the output of fft() depends on whether you use an odd or even number of points for your fft. I think this post nicely summarises how the frequencies are arranged. Have a look at it.

Since you are using an even number of points, the Nyquist frequency, $F_N = F_s/2$, is present in the output of your fft, and is the maximum value in your frequency vector. Your sampling frequency is $N/T=100/10=10\;\textrm{Hz}$, and so your maximum detectable frequency in the vector will be $F_N=5\;\textrm{Hz}$.
Secondly, no matter whether you use even or odd points, the DC value is always the first entry in the fft output.
The frequency resolution $df = F_s/N$ is the separation between fft output bins in frequency space. For you, this will be $0.1\;\textrm{Hz}$.
Finally, it is also most common to display the frequency spectrum in a symmetric, centered way, with the DC frequency in the middle - the function fftshift does this for you.

Accounting for all this, your frequency vector corresponding to your fft output (after you have applied fftshift() to centre it) will be
$$
f = [-F_N : df : (F_N-df)]
$$
There is some discussion about fftshift in this thread if you are interested. Notice that after fftshift() the component corresponding to the Nyquist, $F_N$, is placed at the beginning of the vector, and so is considered a negative frequency from the point of view of fftshift() - there is no Nyquist at the end of the vector in the positive frequencies (then in contrast, for symmetry, the DC component is viewed as a positive frequency).
Plotting the Spectrum with Correct Scaling
Personally, I prefer to deal with these things in terms of the power of a wave (rate of change of energy) instead of the amplitude of a component in a fft bin. This just helps me to keep track of the units, and to keep things physical. The following code calculates the power spectral density, $S_{xx}$, for your signal using a simple  periodogram estimate (I assume your signal is in volts - may or may not be true - but could be anything: temperature, stock prices, etc). Sorry, but my code is in Matlab, not Python, but hopefully you can follow the steps without problem.
N = 100;   % Number of samples
T = 10;    % Record window duration
dt = T/N;  % Sampling period
Fs = 1/dt; % Sampling frequency

t = 0:dt:(T-dt); % Time vector for sampling

% generate samples at the specified times
x = 5*sin(2*pi*3*t - 2); % units: [V]

df = Fs/N; % frequency resoltuion (bin width in frequency space)

% generate frequency vector for 2-sided spectrum (NOTE, this arrangement
% only works for even number of points - otherwise, use f = -(Fs/2-df/2):df:(Fs/2-df/2))
f = -(Fs/2):df:(Fs/2-df);

% Calculate Fourier transform (approximating CFT), and shift DC term to centre
X = fftshift(fft(x))*dt; % units: [V sec]

X((abs(X)<1e-10)) = 0; % kill values below threshold, so phase is well-behaved

% Calculate power spectral density using periodogram estimation
Sxx = (X.*conj(X))/(N*dt); % units: [V^2 / Hz]

figure; stem(f,Sxx)      % Plot power spectral density
figure; stem(f,angle(X)) % Plot phase

The variable $X$ in the code is a discrete approximation of the continuous Fourier transform integral, which is why it is multiplied by the sampling period $dt$ (the approximation is a Riemann sum). Then the periodogram can be calculated from this - again, in a discrete version of its continuous counterpart, which is explained very nicely I think in this answer. Note, Matlab also has a built-in function  periodogram(), which does the calculation as well.
You can see that there are peaks correctly at $\pm 3\;\textrm{Hz}$, matching exactly your input signal (remember the sin wave is actually a sum of complex exponentials).
-- Note that in signal processing, the "energy" is defined as [signal squared multipled by time], and so therefore power is given as [energy per time], which is then [signal squared]. This is why you see "power" in my answer given in units of $[V^2]$. To find the physical power in [Watts] you would simply need to scale by the resistance of whatever the load is you are driving. This is usually set to 1 in signal analysis --
Verifying Parseval's Theorem
You can also verify the scaling of the fft output by checking that the power and energy in both the time domain and frequency domain are equal, as they should be according to Parseval's Theorem. Remember that if you integrate the power spectral density over all frequencies, you should obtain the total power that was present in the input signal. For example, using
Energy_timedomain = sum(x.*conj(x))*dt  = 125  [V^2 sec]
Power_timedomain = Energy_timedomain/T  = 12.5 [V^2]
Energy_freqdomain = sum(X.*conj(X))*df  = 125  [V^2 sec]
Power_freqdomain = sum(Sxx)*df          = 12.5 [V^2]

Note that if you take a longer time window duration $T$ for your input signal then the "energy" will increase (because there is more of it), but the average power will stay the same (because that is the rate of change of energy).
Recovering the Amplitude
The value of each peak in the power spectral density plot is $62.5\; [\textrm{V}^2 / \textrm{Hz}]$. To convert this to power, we need to multiply by the bin width (essentially integrating the power spectral density over that bin). In your case, the bin width is $df=F_s/N = 10\;\textrm{Hz}/100 = 0.1\;\textrm{Hz}$. Therefore, the total power in the input signal at that frequency is
$$ P = 2 \times 62.5\;[V^2/\textrm{Hz}] \times 0.1\;[\textrm{Hz}] = 12.5\;[V^2] $$
where the factor of $2$ is because there are 2 peaks - half of the power in the negative part of the spectrum and half in the positive part. Note that if your input signal is real, then the output spectrum will always be symmetric, and you will often see the negative frequencies (which don't contain any extra information) discarded. If you choose to throw away the negative frequencies then you should multiply the power in the unique positive frequencies that are left by 2 to compensate (but don't multiply the DC or the Nyquist by 2, because there are only one of each of those in the fft output - they don't have negative counterparts).
We know that for sinusoids the power and amplitude $A$ are related by $P = A^2/2$ (again, just setting any load resistance to $1$). Therefore, rearranging this, your amplitude of the 3 Hz wave is
$$A = \sqrt{2\times 12.5} = 5\;[V]$$
Recovering Phase
As Hilmar said, the phase of $-2$ can be recovered by taking the phase value at $3\;\textrm{Hz}$, which is 2.712, and adding $\pi/2$ then subtracting $2\pi$. The $\pi/2$ is because your input signal is a sine wave, which essentially includes a phase shift already because
$$
\sin(2\pi f_0 t + \phi) = \cos(2\pi f_0 t + (\phi - \pi/2)) = \cos(2\pi f_0 t + \phi_{\textrm{fft}})
$$
so you need to add $\pi/2$ to the phase returned by the fft() to get your input phase
$$\phi = \phi_{\textrm{fft}} + \pi/2 = 2.712 + \pi/2 = 4.283.$$
The final $2\pi$ subtraction is because after doing this you are left with something greater than $\pi$, and so you just need to unwrap the phase to bring it back into the range $[-\pi ... +\pi]$.
A FINAL NOTE OF CAUTION
It all works nicely in this case because all of your signal energy falls exactly in a single fft bin, due to the combination of the sampling rate, number of points, time window, and your input frequency wave (see here). Things are more difficult if you have any windowing or zero-padding applied to your input signal before you pass it to fft(). There are some enlightening answers here and here, when you get that far.


Loss of precision due to low-pass filtering in physical units Suppose I have a $10{\times}10$ pixel image, where each pixel represents 1 mm in physical units, and I apply a Gaussian low-pass filter with $\sigma = 5$. Can I make any statements such as "the maximum level of precision in the original image was 1 mm, while it is $x$ mm in the filtered image."

In general, no.  In fact, if the low-pass filter characteristics were well-chosen in the sense of getting the best post-filtering signal to noise ratio, then the amount that you might degrade any such measurements would be minimal.

In part this conclusion comes from my "well-chosen" qualifier.  Roughly speaking (and I can only speak roughly because your problem is so vaguely stated), if you've got a true feature in the object that you're imaging that resolves to an image feature as sharp as your 1-pixel line then unless you have absolutely tremendous noise your $\sigma = 5$ Gaussian filter is poorly chosen.

I would assume that for looking at tissues, you're looking for features that are significantly larger than 0.05mm, and in x-rays you're getting significant white noise -- hence, a filter that smooths out the noise without smoothing out the features that you can actually reasonably expect to be there is only going to make things -- including measurement precision -- better, or at least no worse.

(But if I had a shop accident and had some AWG 36 gauge copper wire embedded in me, then I would presume that your $\sigma = 0.25\mathrm{mm}$ filter may blur the edges too much -- but that's a different problem, and thank goodness that copper wire shows up very strongly on an x-ray).

In general, what you can do is to set out the problem with more specifics: i.e., what's a typical feature you're looking at, how do its edges naturally merge into the surroundings, what's its contrast in the image, and how much noise is there in the image.  With all of these in hand, then you can start asking how much a given filter will impact tasks such as identifying the location of an edge (or, since the subject is mammography, presumably identifying the size of a lesion).


Do FFT-based filtering methods add intrinsic latency to a real time algorithm? In a current audio processing application, I am working entirely in the time domain using IIR filters. This is so I can use an audio buffer of just 3 or 4 samples and can guarantee that there will be some form of audio output within that frame.

If I were to switch to FIR based filters (to reduce filter ringing time as outlined in this question), I would probably need to use an FFT-based overlap-add technique to maintain processing efficiency. As I see it, I would have to fill an FFT frame buffer with samples before I would be able to do the transforms. This would add an intrinsic latency to the algorithm proportional to the number of samples used to calculate the FFT. Am I missing a trick?

You are correct. FFT based processing adds inherent latency to your system. However there are ways to tweak this.

Let's assume you have an FIR filter of length "N". This can be implement FFT-based using the standard overlap add or overlap save method, where the FFT length would be 2*N. Overall system latency will also be roughly 2*N: you need to accumulate a frame of N samples and then while you are accumulating the next input frame, do the math on the current frame. By the time the second frame has been accumulated (time offset 2*N), the first frame is ready to go out. If you have a fast CPU you can speed this up somewhat by using different alignment for input and output frames but that's typically more bother than it's worth.

You can also break down the filter into K smaller blocks of length M, i.e. N = K*M. The FFT only needs to be done over M input samples and the delay and accumulation over the multiple filter sections is done in the frequency domain. Bill Gardner described a few flavors of that here: http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF
It's often referred to as a "Block Convolver".

This allows to basically trade-off latency against efficiency and it provides a continuum between the direct FIR and the full-sized overlap add method. A nice side-effect of the Block Convolver is that the total filter length doesn't have to be a power of 2. For example you can implement a filter of 768 tabs as 6 blocks of 128 without any significant loss of efficiency.


What is intuitively a syndrome in error correction? In error correction for example for Reed-Solomon codes there are computed syndromes for received data in decoding process. What is intuitively a syndrome (I don't mean the mathematical expression for it)?

The word syndrome was introduced into error-control coding theory
and practice (by Hagelbarger) very early on. It is a term borrowed from
medicine: a syndrome is a group of symptoms, e.g. stomachache, nausea,
fever, lack of appetite, which together are diagnostic (or highly suggestive)
of the presence of a specific disease or condition -- in this instance, appendicitis.  Note that the presence of a single symptom is usually not
considered diagnostic of a specific disease, it is the presence of the
syndrome or a collection of symptoms.
This notion in medicine is almost exactly mirrored in
the error-control coding world. If the syndrome is identically
zero (there are no symptoms whatsoever), nothing needs to be
done, and the received word is accepted as correct (no disease
is diagnosed and no treatment is prescribed).
A nonzero syndrome indicates that the received word has errors
in it (the presence of symptoms means that something is not right
with the patient), and in most cases, the syndrome is also diagnostic
of the locations of the errors and what needs to be done to correct them
with the (Reed-Solomon) decoder substituting for the physician's
thought processes (from syndrome, find the error locations and values)
and actions such as prescribing medicines or performing surgery
(correct the errors, fill in the erasures, etc.).  In some
cases, the fact that the syndrome is nonzero indicates the presence
of errors but the decoder is not able  to determine the error
locations and values (and so is unable to correct the errors).
This is called decoder failure.  In even fewer cases, the syndrome
indicates a specific error pattern as the most likely cause but in fact the
actual error pattern is something else. Correcting what the decoder
thinks is the error pattern introduces more errors into the
received word, thus making it "wronger" than before, but
this is not detectable at the receiver. This is not called decoder
malpractice but decoder error or undetected error, and is,
fortunately, a very rare event, much rarer than decoder failure.


Homography valid for patch within the image but not the entire image? The homography I have is non-degenerate (det!=0), and generated from a valid planar pose. When I use it to warp the four corners of image, it returns something like the following:
1 0
3 2

instead of something like
0 1
3 2

where 0 represents top-left corner, 1 top right, 2 bottom left, 3 bottom right. It doesn't follow the clock-wise order anymore, and it's twisted.
The weird thing is, if I apply it to a local patch within the image, i.e. where the plane is, the returned result is valid.
How can this happen? Shouldn't it always return a valid quadrangle?

Are you sure that the points are all in front of the camera (z > 0), if you apply a a projection like $(x, y) = (x/z, y/z)$ then points in a 3d quadrilateral where vertices 0 and 1 have negative z, you could end with a 2d figure like you described


How make FFT when buffer size is not 2^L I still try to understand FFT :)
I know I can do something like that:

float bufferSize = 500.0f;
int correction = log2(bufferSize);
bufferSize= pow(2, correction);

But then I get buffer size 256. So it doesn't use all samples in the buffer. But I've heard there are some methods to make FFT with buffer size for example 500. And use all possible samples.
Could anyone explain me how to do that?
Great thanks in advance.

So I need to develop my question :)

I’ve just read Rabiner, Lawrence R., and Bernard Gold. "Theory and application of digital signal processing."  Of course not all but the fragment titled “A Unified Approach to the FFT”.

And there is about dividing N-point sequence to make matrix. I even implemented in C++ exact example from that book with N=60 points of sequence. But I am not sure if I understand the idea (maybe it’s because English is not my first language).

Firstly in the book they say to make matrix 60=5 X 12, and further they say to pay attention that we can divide each 12 points to make matrix 3 X 4. But although it is said, the further calculations are only for 5 X 12. Why not 5 X 3 X 4?
Is it only for make explanations easier? Or what?

And more important question is what exactly that matrix division give me?
Stanley Pawlukiewicz says about my example (when N=500) to make prime factors 2 2 5 5 5. So as I understand it’s 5 dimensional matrix, 2 X 2 X 5 X 5 X 5? Am I right?

But to make DFTs on that matrix I still need to make 500x500=250000 calculations, so it’s like regular DFT, not any FFT. I could imagine I can make 125 X 4 points FFTs (4 is 2 to the power of 2). But if it’s the idea, why we are not talking about matrix 125 X 4, but instead we are talking about matrix 2 X 2 X 5 X 5 X 5?

In which point here is FFT?

The prime factors of 500 are:

2     2     5     5     5

fftw can do this fast.

If you want to roll your own in:

Rabiner, Lawrence R., and Bernard Gold. "Theory and application of digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc., 1975. 777 p. (1975).

there is a section that shows how to do a 1 dimensional dft in two dimensions.

see my answer

Combine FFT's of shorter length than sample data to get spectrum of all data

pasting the Matlab code from there to here, and taking M=20 N=25 for your problem.

clear all
M=3;
N=32;
x=linspace(0,10,M*N);
X=reshape(x,N,M).'; % read in as rows
Twiddle=zeros(size(X));
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)(k-1)/(NM));
end
end
X=fft(X) % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')
figure(3)
plot(x,'linewidth',2)
title('time series')

FFTW is probably faster.


Which parameters defines a polar plot? I need to convert a polar plot into a short array of 3-5 values. This array should describe the polar plot as best as possible.

For example, I know that I can use the beamwidth which is also known as the half power. This is the angle between the left and right -3 dB point:

I also know that I can use the the ratio between the surface of the polar plot versus the whole 360° plane.

Thus, so long I only have two values in the array:

Beamwidth
Ratio polar plot vs plane

But are there other values that can describe the polar plot ?
And also, how can I have a value which describes second biggest lobe ?

Some examples would be the null to null bandwidth(s) as well as the sidelobe levels of the first/second sidelobes.

Example vector of 5 numbers:
[mainlobe gain, first null location, 1st sidelobe gain, second null, 2nd sidelobe gain]

From this you could produce an estimate for the 3 dB points of each lobe using some "reasonable" formula, i.e. halfway between the peak and the null.


Minimum phase conversion with Cepstrum method, how to scale the result? I am trying to convert a zero phase spectrum (magnidue response curve with zero phases) to a minimum phase spectrum, because I need a totally causal impulse response for FFT spectral filtering, and having an IR devoid of anticausal components allows me to implement a zero latency partitioned scheme.
By now I experimented with two similar methods based on the real cepstrum but my question is inherent to the second method.
In the first method you compute the cepstrum by taking the ifft of the log of the magnitude curve, negate the second half of the cepstrum (actually I rather use a cosine function to have a more graceful transition and less ringing artifacts), take the fft of the result and then instead of computing the exp of the resulting spectrum I simply keep the phase combining it with the original magnitude. This method works well enough despite the resulting IR doesn't always fade to zero.
In the second method you zero the second half of the cepstrum instead of negating it (here too I rather prefer a rised cosine for a more graceful transition), double the result, take the fft and then exponentiate the result.
Here is the problem. The resulting spectrum is uniformly scaled by a variable amount which seems depending on the shape of the original magnitude spectrum, so I cannot compensate by scaling by a fixed value because I don't know the relationship between the resulting scaling factor and the original magnitude spectrum, and I can't definitely have a global level which is dependent on it! I searched but I found very little information around about the min phase cepstrum method. How can I scale the final result properly? I realize this may be a very specialistic topic. Thanks

One possibility is to take the FFT magnitude of both, run a linear regression between the two magnitude vectors, and adjust the gain until the slope of the regression fit is 1.


Is it possible to recover a waveform from spectrograms of magnitude *and* phase? It is my rough understanding that a typical spectrogram discards phase information by squaring only the real-valued part (magnitude) of the STFT (I know it's not quite this simple, but feel free to correct me if I'm wildly wrong here).

While many software libraries (e.g., scipy) will, by default, produce a spectrogram of magnitude, they often also make it possible to instead compose a separate spectrogram for phase. My question is quite simple: if I had both of these matrices, could I reconstruct the original signal?

I suspect the answer to this question is 'not quite'. If so, what part of the process is problematic and are there any changes that could be made to make this goal achievable? I am currently using STFT.real and STFT.imag (where STFT is a 2D matrix for the short time Fourier transform of some signal) and plotting of them separately. While I can reconstruct the original signal using this approach (which I need to be able to do), these plots are not as easy to read as actual spectrograms.

(For instance, I know that squaring the magnitude is clearly a problem for reconstruction as it destroys negative numbers).

Real/imaginary or modulus/phase are two representations of a complex number that carry the same level of information. Then, a STFT is a redundant mapping from a space of functions over a 1D variable (time) onto a space of functions over a 2D variable (time and frequency). Under mild conditions on the window, there are an infinity of inverses, due to the redundancy. In the continuous setting, the analysis ($a$) and synthesis ($s$) windows only need to overall a little (more pedantic: $0 < \int a(t).s(t)dt <+\infty$). In the discrete form, this can be a little more involved, yet there are several optimized inverses for "almost" every sound redundant analysis filter-bank).

So, if for each point of the $(t,f)$ plane, you have the magnitude in the phase, you can recover the waveform (apart possibly at the extremeties). A side effect of the huge redundancy, inversion even remains (approximately) possible under mild conditions if you only have the modulus or the phase.


Fourier series of cycloid What is the Fourier series representation of a cycloid?

The parametric representation of the curve is as follows.
$$
t=\dfrac{\theta-\sin\theta}{\pi}\\
x=\dfrac{1-\cos\theta}{\pi}
$$
The period is $2$, so the coefficients of the complex exponential Fourier series should be
$$
c_n=\dfrac{1}{2}\int_0^2\!x(t)e^{-jn\pi t}\,\mathrm{d}t
$$
and that's where I got stuck. I performed substitution
$$
t=\dfrac{\theta-\sin\theta}{\pi}\implies\mathrm{d}t=\dfrac{1-\cos\theta}{\pi}\,\mathrm{d}\theta\\
c_n=\dfrac{1}{2}\int_0^{2\pi}\!\left(\dfrac{1-\cos\theta}{\pi}\right)^2e^{-jn(\theta-\sin\theta)}\,\mathrm{d}\theta
$$
but don't know how to integrate this monster. Is there any other way?

The Fourier series of the cycloid can be expressed in terms of the Bessel functions of the first kind:
$$J_n(x)=\frac{1}{\pi}\int_0^{\pi}\cos(nt-x\sin t)dt,\qquad n\in\mathbb{Z}\tag{1}$$
Using the cycloid parameterization
$$y(t)=1-\cos t,\qquad x(t)=t-\sin t\tag{2}$$
which results in a period of $2\pi$ and a maximum value of $2$, the Fourier series of $y(t)$ as a function of $x$, referred to as $f(x)$, is given by
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{
\begin{align}f(x)&=\frac32+\sum_{n=1}^{\infty}\frac{J_{n+1}(n)-J_{n-1}(n)}{n}\cos(nx)\\&=\frac32-2\sum_{n=1}^{\infty}\frac{J'_{n}(n)}{n}\cos(nx)\end{align}}\tag{3}$$
where $J'_n(x)$ is the derivative of $J_n(x)$ w.r.t. $x$.
The figure below shows the cycloid and its Fourier series approximation according to $(3)$ using the first $20$ coefficients in the sum:

Proof:
Let's consider the real-valued Fourier series. Since $f(x)$ is even, all sine coefficients vanish and we get
$$f(x)=a_0+\sum_{n=1}^{\infty}a_n\cos(nx)\tag{4}$$
with
$$a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)dx\tag{5}$$
and
$$a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos(nx)dx,\qquad n>0\tag{6}$$
With $dx=(1-\cos t)dt$ we obtain
$$a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi}y(t)(1-\cos t)dt=\frac{1}{2\pi}\int_{-\pi}^{\pi}(1-\cos t)^2dt=\frac32\tag{7}$$
and
$$\begin{align}a_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}(1-\cos t)^2\cos(nt-n\sin t) dt\\&=\frac{2}{\pi}\int_{0}^{\pi}(1-\cos t)^2\cos(nt-n\sin t) dt,\qquad n>0\tag{8}\end{align}$$
where the last equality follows from the fact that the integrand is even.
Expanding
$$(1-\cos t)^2=1-2\cos t+\cos^2 t=\frac32-2\cos t+\frac12\cos 2t\tag{9}$$
and using $\cos(\alpha)\cos(\beta)=\frac12[\cos(\alpha-\beta)+\cos(\alpha+\beta)]$, the integrand in $(8)$ can be rewritten as
$$\begin{align}(1-\cos t)^2\cos(nt-n\sin t)=\frac32\cos(nt-n\sin t) -&\\\big[\cos((n-1)t-n\sin t) + \cos((n+1)t-n\sin t)\big]+\\\frac14 \big[\cos((n-2)t-n\sin t) + \cos((n+2)t-n\sin t)\big]\tag{10}\end{align}$$
Plugging $(10)$ into $(8)$ and using the definition of the Bessel function $(1)$, we can write
$$a_n=3J_n(n)-2\big[J_{n-1}(n)+J_{n+1}(n)\big]+\frac12\big[J_{n-2}(n)+J_{n+2}(n)\big],\qquad n>0\tag{11}$$
This can be further simplified using the recurrence relation (10.6.1) for $J_n(x)$:
$$J_{n-1}(n)+J_{n+1}(n)=\frac{2n}{x}J_n(x)\tag{12}$$
which can be used to eliminate $J_{n-2}(n)$ and $J_{n+2}(n)$ in $(11)$, and which finally results in
$$a_n=\frac{J_{n+1}(n)-J_{n-1}(n)}{n},\qquad n>0\tag{13}$$
Since (10.6.1)
$$J_{n-1}(x)-J_{n+1}(x)=2J'_n(x)\tag{14}$$
where $J'_n(x)$ denotes the derivative of $J_n(x)$ with respect to $x$, the coefficients $a_n$ can also be expressed as
$$a_n=-\frac{2J'_n(n)}{n},\qquad n>0\tag{15}$$
Eqs $(7)$, $(13)$, and $(15)$ establish the result $(3)$.


Converting a pitch bend (MIDI) value to a "normal" pitch value I am trying an offline LPC analysis and synthesis using the rt_lpc (real-time LPC) implementation as given here. There are some functions within the program that can be used for an offline LPC analysis and synthesis program. The rt_lpc code is mostly meant for composers etc., which means that it uses a MIDI or glottal pulse input. There is some code in the program that converts a pitch value (obtained from auto-correlation) to the bend value. The relation is given as follows:

    pitch = (int)( Stk::sampleRate() / midi2pitch[ananya.data[1]] ) /
                            pow( 1.0653f, bend/64.0f*11.0f );
                    power *= ananya.data[2] / 64.0f;
bend = ge.data[1] / 128.0f + ge.data[2] - 64;

Some points:
a. ananya is an object of type MidiMsg that seems to be populated on the fly.
b. ge is also an object of type MidiMsg that also seems to be populated on the fly.
c. These two objects are created when the program is run in real-time (mine is an offline version that runs selective parts of it in a main.cpp of my own.

My question(s):
1. What is pitch bend?
2. How do I convert a pitch value into a bend value? What are the mapping relations?

I have googled for solutions but found no clear answer.

At the core of MIDI is a representation of music as discrete note events, each of those having a static pitch. This is perfect for representing music as played on keyboard instruments. You can convert any frequency corresponding to a note on the tempered scale into a MIDI note number, using:

$69 + 12 \times \log_2 \frac{frequency}{440}$

Under the assumption that the MIDI receiver is calibrated for A4 = 440 Hz.

This representation is alright for piano music, but the problem is how to represent pitches which are not mapped to the tempered scale (non-western music, non musical sounds), and how to represent pitch variations over the duration of a note (glissando, vibrato).

This is done in MIDI by using "pitch bend messages" which instruct the synthesizer to shift the pitch of the currently played note by a small interval. Most synthesizers are calibrated by default for +/- 2 semitones over the course of the pitch bend message range (0 .. 16383). 8192 corresponds to no pitch bending - the emitted pitch is exactly that of the note value. The mapping between the pitch bend value and the frequency shifting ratio is given by:

$\frac{f_{emitted\_note}}{f_{note\_message}} = 2 ^ {\frac{pitchbend - 8192}{4096 \times 12}}$

You can thus get the frequency of a note played by a syntheszier from the following formula:

$440 \times 2^{\frac{note - 69}{12.0} + \frac{pitchbend - 8192}{4096 \times 12}}$

Where note is the 7-bits MIDI note number of the last received Note On message ; and pitchbend is the 14-bits value of the last received Pitch bend message. A synthesizer starts with its pitch bend register set to 8192, and this value is also reset during the reception of a "Reset all controllers" message.

Let us take the following example. You want to express a flute trill with the following frequency trajectory: 500 Hz, 510 Hz, 500 Hz, as MIDI messages.

The base note number is:

$round(69 + 12 \times \log_2(500 / 440)) = 71$.

So you send a "note on" message with note# equal to 71. This is equivalent to a pitch of:

$440 \times 2 ^ {(71 - 69) / 12} = 493.88$

Which is the nearest pitch on the tempered scale. You need to send a pitch bend message to raise the pitch by a factor of:

$\frac{500}{493.88} = 1.0124$

And get your 500 Hz. The corresponding pitch bend value is:

$round(8192 + 4096 \times 12 \times log_2 1.0124) = 9065$

To get your 510 Hz, the pitch bend value would be:

$round(8192 + 4096 \times 12 \times log_2 \frac{510}{493.88}) = 10469$

So your sequence of MIDI messages for 500, 510, 500 Hz would be:

NOTE 71
PITCH BEND 9065
...
PITCH BEND 10469
...
PITCH BEND 9065

You can think of the MIDI note number as the "integral" part of the pitch ; and the pitch bend as a redundant "fractional" part of the pitch.


Redistributing Color in a RGB Image According to a Gaussian Distribution I haven't done this stuff in a while.

If I have an image $I$ I can equalize the histogram of the image using some opencv procedure, it's already defined. Equalizing an histogram means essentially to construct a function such that applied to each pixel $p \in I$ the new value it's like it had been sampled from a uniform distribution, assuming gray scale image for such explanation.

I don't remember if what I want to do has a specific name, but I kind of remember it was possible, but suppose I have a gaussian distribution (instead of a uniform one), say $\mathcal{N}(0,\Sigma)$, and I want to re-distribute the colors according to such distribution.

Is there a way to do this in opencv?

After you equalize the histogram you can think of your data as a stream of variables $ {X}_{i} $ where $ X \sim U \left[ 0, 1 \right] $.

Now all you need is to transform samples of Uniform Random Variable into Gaussian Variable.
You should do that by applying the Inverse CDF of Gaussian Distribution.
Basically applying the Inverse Transform Sampling method.
For Normal Distribution easy choice would be the Box Muller Transform.

After you apply this, you can just multiply each RGB pixel in the Cholesky Factorization of $ \Sigma $ to get exactly what you need.


Image processing, recognizing a small feature in a larger image I am trying to write an image processing program to recognize bubbles in oil.
It has been suggested I try computing the convolution of the image and an image of a typical bubble.
i.e. : ifft(fft(a).fft(b)) will have high peaks near bubbles.

where:
ifft : inverse fourier transform
fft : fourier transform
a : the image
b : kernel image (bubbles)

Is this a good way to solve this problem? Also, for the dot product to work, I need the images to be the same dimension. This is not likely to be the case.

What has been suggested to you is called template matching. If you don't mind re-using code (and working with C++), you can find an implementation in the OpenCV library.

A simpler alternative (that I would start with) would be to work with generalized Hough transform, but that will require some pre-processing first by detecting the edges in the image.


Why don't DSLR cameras default to png instead of jpeg? I'm working with a new DSLR camera for a Computer Vision application and I noticed most of them default to saving in the lossy JPEG format. Initially, this caused some problems in my application and I started working with the RAW format, and saved them as png. I wondered what was the reason behind defaulting to JPEG as the compression encoding. Wouldn't it be better to support png as a default format due to its lossless quality and compression factors? (I assume DSLR users like to have quality and compression both). I'm curious as to whether this has a scientific reason I have not thought of.

Its force of habit and good enough. People were using JPEG's in their digital cameras before, so they stuck with it, and its plenty fine for viewing (and it produces pretty small file sizes as well). You're free to turn to some other compression method your camera supports.

Note that the RAW files often contain a lot more information than what a PNG would have anyway (like a bunch of stuff about the sensor), which can be exploited in programs like photoshop for editing, so from that perspective, you may as well have the RAW rather than the PNG (especially since the RAW is rather space efficient to begin with). You can also put in things specific to your camera there if you want.

Note that many cameras have the option to shoot RAW+JPEG so you simultaneously get both, which is good for viewing (use the JPEG) and editing (use the RAW).


Why should one use windowing functions for FFT? So I just revised my pitch calculation algorithm using a harmonic product spectrum algorithm. I was just curious about why this explanation of Harmonic Product Spectrum states that you need to implement a Hanning Window to the data set. What would be the effect of implementing other Window functions on a data set (and then FFTing it)? Which Windowing function is actually the best for frequency detection? Here are the relevant methods I have used in my code:

/**
 * Calculates the Frequency based off of the byte array,
 * @param bytes The audioData you want to analyze
 * @return The calculated frequency in Hertz.
 */
private int getFrequency(byte[] bytes){
    double[] audioData = this.bytesToDoubleArray(bytes);
    audioData = applyHanningWindow(audioData);
    Complex[] complex = new Complex[audioData.length];
    for(int i = 0; i<complex.length; i++){
        complex[i] = new Complex(audioData[i], 0);
    }
    Complex[] fftTransformed = FFT.fft(complex);
    //return calculateFrequency(fftTransformed);
    System.out.println("Max size:" + (fftTransformed.length*getFFTBinSize(fftTransformed.length)/4));
    return calculateFundamentalFrequency(fftTransformed,4);
}

private double[] applyHanningWindow(double[] data){
    return applyHanningWindow(data, 0, data.length);
}

private double[] applyHanningWindow(double[] signal_in, int pos, int size)
{
    for (int i = pos; i < pos + size; i++)
    {
        int j = i - pos; // j = index into Hann window function
        signal_in[i] = (double)(signal_in[i] * 0.5 * (1.0 - Math.cos(2.0 * Math.PI * j / size)));
    }
    return signal_in;
}

/**
 * Harmonic Product Spectrum
 * @param fftData
 * @param n
 * @return
 */
private int calculateFundamentalFrequency(Complex[] fftData, int n){
    Complex[][] data = new Complex[n][fftData.length/n];
    for(int i = 0; i<n; i++){
        for(int j = 0; j<data[0].length; j++){
            data[i][j] = fftData[j*(i+1)];
        }
    }
    Complex[] result = new Complex[fftData.length/n];//Combines the arrays
    for(int i = 0; i<result.length; i++){
        Complex tmp = new Complex(1,0);
        for(int j = 0; j<n; j++){
            tmp = tmp.times(data[j][i]);
        }
        result[i] = tmp;
    }
    //Calculates Maximum Magnitude of the array
    double max = Double.MIN_VALUE;
    int index = -1;
    for(int i = 0; i<result.length; i++){
        Complex c = result[i];
        double tmp = c.getMagnitude();
        if(tmp>max){
            max = tmp;;
            index = i;
        }
    }
    return index*getFFTBinSize(fftData.length);
}

The FFT can only be performed over a limited chunk of data. The basic math is based on the assumption that the time domain signal is periodic, i.e. your chunk of data is repeated in time. That typically results in a major discontinuity at the edges of the chunk. Let's look at a quick example: FFT size = 1000 points, Sample Rate = 1000 Hz, Frequency resolution = 1Hz. If you have a 10 Hz sine wave you don't have a discontinuity since exactly 10 periods fit into your FFT window and the values (and derivatives) at the edges are the same. The FFT of this siganl will be zero except for a single value in bin #10. This also works for a 11 Hz sine wave just as well. However, for 10.3 HZ sine wave you end up with a lot of discontinuity and the FFT will have energy in all bins with maximum around 10 or 11 and then "skirts" that roll trail off to the sides. So a small change in frequency results in a massive change in the FFT picture.

Windowing is used to avoid this: Windows make sure that the data at the edges are zero, so there is no discontinuity. However multiplication in the time domain is convolution in the frequency domain and that results in widening of spectral lines and also in side lobes. Teh choice of window controls the trade offs between main lobe width and side lobe spacing and height. Your application specific requirements determine what window to use and there are dozens of choices. Hanning is just one of them. It's basically "the window of choice if you don't have any better ideas". Personally, I prefer Kaiser windows as they have a continuous parameter that can control the window behavior over a wide range.

In general FFT is not a great method for pitch detection. For most audio signals the maximum in the spectrum is NOT the fundamental (typically harmonics have higher energy), in order to get decent resolution you need a long pieces of data but that makes the algorithm very slow and sluggish to respond to changes. Much better options are phase looked loops, delay looked loops, autocorrelations, max/min tracker, zero  crossing tracker, etc.


Applying Convolution (In 2D) Theorem Swaps Quadrants For educational purposes I implemented the DFT and inverse DFT for images using OpenCV. Applying the DFT to an image and taking the inverse DFT of the spectrum yields the original image, so this works. I do not rearrange the spectral image from $ [0,N) $ to $ [-N/2, N/2) $.

Now, using the convolution theorem I want to filter an image. I use a 3x3 averaging kernel and pad it with zeroes to the same size as the input image. I then take the DFT of it, and do a piecewise multiplication of the image spectrum and the kernel spectrum.

The inverse DFT of the resulting spectrum yields an interesting result: The image is filtered (smoothed) correctly, however the quadrants are swapped. Interestingly, if I filter it twice (or any even number), the quadrants are properly arranged. This works for both

$$ \mathcal{F}^{-1}\{ (\mathcal{F}\{I\}\mathcal{F}\{H\})\mathcal{F}\{H\} \} = $$
$$ \mathcal{F}^{-1}\{ \mathcal{F}\{I\}(\mathcal{F}\{H\}\mathcal{F}\{H\}) \} $$

with the image $ I $ and kernel $ H $. This makes sense to me since convolution is associative.

What is the cause of the swapped quadrants?

It's probably due to how you padded the filter before DFT'ing it.  I bet you had it in a corner instead of in the middle.  Since convolution via DFT is circular convolution, it's not so much of a quadrant swap as it is a wrap around.  I believe that if you put the non-zero taps in the middle that it will fix the problem.


Distance along beam calculation I measure vibration along a beam, 74m long. I have 2 sensors either side of the beam and measure vibration for 1 second ($t= 1s$) at 9000 Hz sampling frequency. I measured speed to 12000 m.s.

I'm picking up features along the beam, representing discontinuities along the beam. I would like to convert the time vector to distance, so that everything which happens in $t$ is now distance along the beam ($x$ axis) and $y$ axis represents the discontinuity feature.

This is in MATLAB.

My attempt:

Consider the following diagram showing the beam with length L, sensors A and B on each side of the beam and an arbitrary position x representing the location of the discontinuity. First consider the case of a signal source emanating from x equally toward A and B, with zero loss in the beam. The overall distance of the beam itself is such that the delay from A to B does not cause the phase of a signal at it's sampling rate to exceed $2\pi$ (otherwise alias solutions will exist). Smaller phases at this frequency would be allowable, but this will allow for the best precision.

Thus similar to the transfer function of any cable the transfer function from x to A over the normalized radian frequency of $0$ to $2\pi$ has a magnitude of $1$ and a linear phase that extends from $0$ to $-\frac{2\pi}{x/L}$. Similarly the transfer function from x to B would have a magnitude of $1$ and a phase that extends from $0$ to $-\frac{2\pi}{(L-x)/L}$

These are given as:

$$H_{xA}(\omega) = e^{-j\omega L/x}$$

$$H_{xB}(\omega) = e^{-j\omega L/(L-x)}$$

And the ratio of these two is given as:

$$H(\omega) = H_{xA}(\omega)/H_{xB}(\omega) = e^{-j\omega L (L-x)/x}$$

With $\omega$ as the normalized radian frequency in units of radians/sample ($\omega$ extends from $0$ to $2\pi$ with $2\pi$ representing the sampling rate assuming the beam is sized such that the phase delay of the sampling rate was $2\pi$ (or this can be modified accordingly for shorter delays).

Consider a single tone given as $e^{j \omega_n}$ that is originating at $x$. The signal received at A compared to the signal received at B is given by $H(\omega_n)$ to be $e^{-j\omega_n L (L-x)/x}$.

This would appear at A as $e^{-j\omega_n L/x}$ and at B as $e^{-j\omega_n L/(L-x)}$.

Specifically if we followed the guidelines above on length versus sampling rate, the phase of this comparison signal is unique to $\omega_n$ and $x$ for a given $L$ as $\omega_n L/(L-x)$.

Given that and to the extent we can model the discontinuity as a source that emanates equally in all directions, I would recommend the following approach:

Measure A as $a(t)$ and compute the DFT to be $A(\omega)$

Measure B as $b(t)$ and compute the DFT to be $B(\omega)$

Compute the comparison signal $R(\omega) = A(\omega)/ B(\omega)$

Compare the phase of $R\omega$ to the phase of $H\omega$ and solve for x for each bin:

$$arg(R(\omega)) = \omega L/(L-x)$$

$$x(\omega) = L\bigg(  1 - \frac{\omega} {arg(R(\omega))} \bigg)$$

Where $\omega = 2\pi k/N$ for DFT bins $k = 0$ to $N-1$.

This approach can be extended for lengths with phase shifts less than $2\pi$ at $f_s$ as well as actual transmission loss factors (since lossless was assumed here).

Note that two or more discontinuities that react to the same frequency will not be resolvable: their relative phase contributions will add and so will therefore appear to be coming from another location (such as images in a reflection).

We can calculate this for your specific example:

$$f_s = 9000 \text{ Hz}$$
$$L = 74 \text{ m}$$
$$v = 12000 \text{ m/s}$$

The wavelength in meters of the sampling frequency is $v/f_s$ = $12000/9000 = 1.33 \text{ m}$

This means a the maximum length of the beam is 1.33 meters for a unique solution over the full sampling rate given this sampling rate given. However all frequencies up to point where their delay would be $2\pi$ would be resolvable and this frequency would be:

$$f_{max} = \frac{v}{L} = \frac{12000}{74} \approx 162.2 \text{ Hz}$$

Given this, there is no reason to sample higher than twice this frequency but it doesn't hurt to do so. We would expect to see repeating images of everything from DC to 162.2 Hz at the higher frequencies.


How to apply an anti-aliasing filter before downsampling I have an eeg signal with two channels (f3m2 and f4m1) which is divided into epochs. I want to augment the data by taking every 5th sample. I can re-use the discarded samples to create 5 versions of the data (take every 5th sample starting at sample 0, then take every 5th sample starting at sample 1, etc.). I have written my own code to downsample the data starting at sample 0, 1 etc. but I know that in order to avoid aliasing I need to first apply a low pass filter.
How do I decide what filter to use (FIR, IRR, etc.)? How do I decide which cutoff frequency to use? How would I implement this?
My code is in python and I have been considering using the scipy.signal library (but am open to any python libraries).

Anti-aliasing filtering is applied just as any other LTI filtering: If your input data is $x[n]$, and the impulse response is $h[n]$, then your output will be
$$y[n] = x[n] \star h[n] $$
where $\star$ is the convolution operation, a.k.a. the anti-aliasing filtering in this context.
Your impulse response $h[n]$, ideally, corresponds to a lowpass brickwall filter, with a cutoff frequency of
$$ \omega_c = \frac{\pi}{M} $$ radians per sample, where $M=5$ is the downsampling ratio in your particular application. You may design $h[n]$ by any suitable method, but for maintaining high accuracy, you better use tight specs on the desiged filter, this can be achieved by a sufficiently long FIR linear-phase filter, or similar IIR filter. The latter will be more efficient to implement but less accurate due to its nonlinear phase.
Then you can decimate $y[n]$ by $M$ as $$v[n]= y[Mn]$$
to  get the downsampled sequence.


Correlating accelerometer data I am getting a signal from a device. This signal contains a non-random noise signal in it. I want to remove that noise from the signal. This noise is correlated to the motion the device is having. So now I've attached a accelerometer to the device to measure the motion produced by the device. I'm getting the data of X, Y, Z-axis from 3-axis accelerometer.

Now here is the process I thought to remove the noise.First get a motion signal from the 3 axis data which correlate to the motion noise present in the main signal. Then pass the main signal and the correlated noise signal through a LMS filter to get a filtered data. I've already designed a LMS filter but I'm stuck at how to get the motion signal(from the data of the 3-axis from accelerometer) which will correlate to the motion noise present in the main signal.

Can anyone help me with this plz??

Motion can be described using displacement with time, velocity or acceleration. To get these quantities from the accelerometer data the following link may help.
http://www.chrobotics.com/library/accel-position-velocity


Represent DFT coefficients with respect to Continuous time-Fourier series coefficients Does anyone know how to represent the Discrete Fourier transform (DFT) coefficient, $X[k]$, with respect to the Continuous time-Fourier series (CT-FT) coefficient, $X_k$? I come to the conclusion as $X[k] = NX_k$. $\\\\$
As we know Fourier series is for a periodic signal. So DFT coefficient is just a sampled version of the Fourier series coefficient. Let $P$ be the period of CT signal $x(t)$, and $N$ be the sampled period of $x[n]$. We sampled $x(t)$ by $T=P/N$ as sample period. From $$x[n] = x(nP/N)\\=\sum_{k=-inf}^{inf}X_k e^{2\pi*k/P*(nP/N)}\\ =\sum_{k=-inf}^{inf}X_k e^{2\pi*kn/N}$$
Since for DFT, $$x[n] =1/N\sum_{k=-0}^{N-1}X[k] e^{2\pi*kn/N}$$, we can say that $X[k]=NX_k$ for $k=0,1,...,N-1$? $$\\\\$$Is my interpretation correct?

No, it's not that simple. In general, sampling will introduce aliasing. The correct derivation is as follows. Let $T$ be the period of the continuous-time signal $x(t)$, and let $T_s$ be the sampling period satisfying $T=NT_s$ with integer $N$. It follows that the sampled sequence $x_d[n]=x(nT_s)$ is also periodic with period $N$.
We have
$$x(t)=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{T} t}\tag{1}$$
For the sampled sequence we obtain
$$x_d[n]=x(nT_s)=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{T} nT_s}=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{N} n}\tag{2}$$
where we've used $T_s=T/N$. Note that the term $e^{j\frac{2\pi k}{N} n}$ is $N$-periodic in $k$, so we can rewrite $(2)$ as
$$x_d[n]=\sum_{k=0}^{N-1}\sum_{l=-\infty}^{\infty}c_{k+lN}e^{j\frac{2\pi k}{N} n}\tag{3}$$
Comparing Eq. $(3)$ to the IDFT
$$x_d[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi kn}{N}}\tag{4}$$
yields
$$X[k]=N\sum_{l=-\infty}^{\infty}c_{k+lN}\tag{5}$$
Form Eq. $(5)$ it is obvious that the DFT coefficients $X[k]$ are a (scaled and) aliased version of the Fourier coefficients $c_k$ of the continuous-time function.


Intuition about independent signals Given is this Wiener filter:

From this we take \begin{equation}
x[k]-a x[k-1]=v[k]
\end{equation}

$v(k)$ is assumed to be a white gaussian noise.

In the textbook it is then stated that

  The input $v[k]$ at time $k$ and the output $x[k − 1]$ at time $k − 1$ are statistically independent. Thus, $E\{v[k]x[k − 1]\} = 0$.

My intuition tells me thats wrong. The equation above can also be written as
\begin{equation}
x[k-1]=\frac{v[k]-x[k]}{-a}
\end{equation}

in which I see that $x[k-1]$ depends on $v[k]$ and is therefor not independent of it.

Where is the mistake in my thoughts? I guess I look at it the wrong way. I'd like to build a better intuition for problems like this. They always mention it so easily in these textbooks, and I need hours to verify everything.

Think it this way; assuming $x[-1] = 0$, then recusively compute the output $x[k]$ for $k \ge 0$ such as

$$\begin{align}
x[0] &= v[0] \\
x[1] &= a x[0] + v[1] \\
x[2] &= a x[1] + v[2] = a^2 x[0] + a v[1] + v[2] = a^2 v[0] + a v[1] + v[2] \\
x[2] &= \left( a^2 v[0] + a v[1]  \right) + v[2] \\
... &= ... \\
x[k] &= \left( a^k v[0] + a^{k-1}v[1] + ... + a v[k-1] \right) + v[k] \\
x[k] &= a ~ x[k-1] + v[k] \\
\end{align}
$$

As you can see at the bottom two lines, the output at time $k$ uses all input samples up to time $k-1$ and the new input sample $v[k]$. Since $v[k]$ is declared as an independent random process, then all of its random variables $v[0],v[1],v[2],...,v[k]$ are independent from each other. Then the current input $v[k]$ will be independent from the past output $x[k-1]$, as that output does not contain any contributions from the input sample $v[k]$;

So what's the illusion in this equation
\begin{equation}
x[k-1]=\frac{x[k]-v[k]}{a}
\end{equation}

which seems to suggest a dependency between $x[k-1]$ and $v[k]$ ?

If you consider the causal solution of the difference equation as outlined recursively above, then it's seen that the current output sample $x[k]$ contains a contribution from current input sample $v[k]$ and plus all past input contributions from $v[0],v[1],...,v[k-1]$ but then the difference $x[k]-v[k]$ will cancel out the presence of the current input $v[k]$ at the right hand side of the equation. Hence the left hand side $x[k-1]$ will be actually depending on the past inputs alone and not on the current input (which apprently explicitly shows up there but is actually implicitly cancelled by the difference of $x[k]$ and $v[k]$).


How to use MATLAB's $\tt xcorr$ to find index of a section inside a whole recording? Trying to find the index offset of a short section of an audio I have taken out of my recording with cross-correlation, but this method sometimes succeeds finding the index of the short section sometimes fails, seems that it depends which sections are being compared and the length of the section made.

My question is what makes the xcorr lose its accuracy ? And if is it dependable ?

clear all;
close all;

[A,fs]=audioread('hendrix.wav'); % some random audio

offset=36644-(784/2); % here I define the offset for my selection.

B=A(offset:offset+784); %here I make selection from A 784 samples long...

[c,lags] = xcorr(A,B);  % c is the correlation, should have a clear peak

c=c((round(length(c)/2):round(length(c)))); %positive coefficients

[ff,ind]=max(c); % finding the index of the offset

M=A(ind:ind+(length(B))-1); % the selection from A based on the new index

plot(B);
hold;
plot(M);

Tried with taking absolute value of the maximum ? I have updated the code, see if that works for you...

[ff,ind]=max(abs(c)); % finding the index of the offset

M=A(ind:ind+(length(B))-1); % the selection from A based on the new index

plot(B);
hold;
plot(M);


How to generate a signal which frequency are obtained? I want to generate a signal which frequency are obtained.
For example,I want the frequency domain similiar to this one:

So I use the following code:

fs =44100;
t = 0:1/fs:4;
f = 800 + 500*sin(2*pi*t);
x = sin(2*pi*f);
plot(t,x)

But using Adobe Audition to visualize its spectrogram.

It looks different from my hope.
So how can I generate a signal if already have freqency?

Your code is wrong because when you want to get a signal whose frequency changes with time, then the argument of the sine you generate is given by the integral of the function you want to generate in frequency (see https://en.wikipedia.org/wiki/Chirp). That's why, in your case, your f function should be:

$f(t)=\int_0^t800+500\sin{2\pi t}=\phi_0 + 800t -\frac{500}{2\pi} \cos(2\pi t)$

If you set the zero phase $\phi_0=0$, then you would obtain the signal you are looking for. Apart from that, take into account that if you just want that signal in the positive part of the spectrum of your double sided spectrum, then you should do $e^{if(t)}$, because if you just take the sine of $f(t)$ then the signal would be replicated in your negative side of the spectrum.


Echo cancellation supporting long delays and without frequency domain processing AEC algorithms mostly rely on LMS adaptive filtering, i.e. you update FIR filter coefficients then perform the filtering. Theoretically, the FIR must be as long as the maximum echo length you want to cancel. For instance to cancel delays up to 500ms on a 48kHz signal, you'll need a 24000 point FIR.
When your memory and processing power limitations make it so you can neither afford to perform 24000 MACs per processed sample, nor to use FFT-based fast convolution algorithms, is there a way around for cancelling such potentially high length echoes in a more affordable way, given delay is unknown and potentially variable?

I was wondering maybe about some other algorithm running in parallel that could assess the approximate delay, then use an adaptive length delay line + a shorter adaptive FIR filter (up to a few hundreds taps is OK)

Does this make sense? Any other neat approach to suggest?

I believe you logic is reasonable. The fact that there are companies offering algorithms for long echo tail (even above 1 second) with reasonable performance indicates that an algorithm similar to your suggested logic can work.


Why a Large Diaphragm Microphone is limited to lesser frequency range Here Above are the frequency responses of different size microphone(condenser)

what actually is the relation between the microphone capsule size and frequency response ?

  The conclusion is that a large diaphragm microphone will have a more limited frequency range than a small diaphragm. This is illustrated in the graphs

i found some sites explaining the reason but i did not get the reason well, can any one illustrate for me plz

There are typically two reasons: First, a larger diaphragm is also heavier and hence it's inertia is bigger, which slows it down at higher frequencies (but also provides for significantly higher max SPL handling).

The main reason, however, is the relationship between wavelength and diaphragm size. The microphone basically "integrates" the sound pressure over the spatial area of the diaphragm. For a plane wave, the sound pressure is the same everywhere, but in a diffuse field that's not the case anymore. The wave length at 10 kHz is about 34 mm. So there are points on the diaphragm that are half a wavelength apart and it's entirely possible that one point sees positive pressure and the other negative pressure. These can cancel each other out (to some extent).


How are signal bandwidth and MSPS related? In various software defined radios, there are three important parameters to set when receiving: frequency, bandwidth and MSPS (million samples per second?). What does it mean to receive with a triple of parameters $(freq, b, s)$ (freq-frequency, b-bandwidth, s-MSPS)? Radio listens to frequencies in $[f-b/2,f+b/2]$ range, divides this range into $s$ chunks and transmits the sampled data to the PC?

How to calculate $s$ for a given $b$, to avoid aliasing?

Am I right that MSPS has nothing to do with the Nyquist theorem?

The need for modulation:

Your voice signals generally lies in the frequency range 1kHz to 4kHz and music lies in the
range 20Hz to 20Khz. Say your locality has 4-5 AM broadcast stations. Now all the AM broadcast stations cannot transmit their content as is because if all the stations transmit at the same time there will be signal interference. Also, for efficient radiation of electromagnetic energy the radiating antenna should be of the order of the wavelength of signal radiated. If the signal is 20Hz-20kHz the length of antenna has to be of the size of kilometers. So, they place the information carried by 20Hz-20kHz at higher frequencies. Every station places its their content in their own band of frequencies. Like one AM station will operate on 100Mhz with their band of operation being 100Mhz - 20Hz and 100Mhz+20Khz and the second one at 105Mhz with their band of operation being 105Mhz - 20Hz and 105Mhz+20Khz and so and so forth. The process of placing the frequency content from one spectrum on to another is called modulations. Now this information is radiated in to free space by an antenna at the station.

Receiver:

The receiver has a tuning knob, which lets you select a band of frequencies. Say, you have tuned your receiver to receive station 3 which is operating at 110 Mhz. You will select all the frequencies  in the range 110Mhz - 20Hz and 110Mhz+20Khz. The signal being present in the spectrum 110Mhz - 20Hz and 110Mhz+20Khz is brought back to 20 Hz to 20kHz using a process which reverses the operation in encoder side. This is called demodulation. This signal is ready to hear.

For a simple modulation-demodulation read this

Finally, if you want to post-process the signal before listening like reducing the noise, increasing bass/treble, you have to filter the signal.

Software-Defined Radio:

In a software defined radio, there will be a RF receiver which does the job of the receiver mentioned above. The message signal that was recovered is then converted to digital samples using an analog-to-digital converter. The analog to digital converter has 2 steps, 1: Sampling at rate above or equal to Nyquist and 2. Quantization. These digital samples are processed using a software on a pc/embedded system. You can apply filters on these samples. These digital/discrete samples is then converted to analog using a D/A Converter. Main step of a D/A converter is interpolation.

Theoretically, in an ideal system, it is sufficient to sample the signal at Nyquist rate. Before play back, an ideal low pass filter will recover the message/baseband signal with out any information loss. But such a design is not feasible.

At a low sampling rate, the interpolation filter (low-pass filter) will not be able to interpolate intermediate samples effectively. However if  the sampling rate is increased you can have better reconstructed message signal. Hence at higher msps you must receive a better quality signal than at low msps. However, the quality you perceive will become negligible after certain point as your ear might not be able to negotiate the difference.

I believe lowering the msps below Nyquist can introduce cross talk. That doesn't mean that you will receive neighbor radio stations because they got clipped at the RF receiver stage. The cross talk you receive is because of frequency warping.

In the above description i am speculating  that sampling happens on the demodulated signal not on the modulated signal. Because demodulation can be easily achieved via an envelope detector for AM or a phase locked loop for FM. Demodulation in digital domain can be costly as sampling a high frequency signal at nyquist rate and operating at such high bandwidths can be resource consuming.

Also no FM station operates at Ghz. Giga is for satellite. 100 Mhz for FM. and Khz for AM.


Confusion regarding frequency of Discrete Time Periodic Signal While Analyzing the DFT Plot for a signal with N samples, say we find a peak in the magnitude of the DFT plot at index $k$, this implies that our signal has a high amount of similarity with an exponential phasor rotating at $\omega = (2\pi/N)k$ or equivalently a period of $N/k$ samples.
But since we are in discrete time, the fundamental period must be an integer i.e. let the fundamental period be $T$ samples
$ \implies T/m = N/k$ (where T and m are relatively prime)
Hence shouldn't the period be $Nm/k$ and hence the frequency would be $(2\pi/N)(k/m)$  ? (where $m$ is the smallest integer that makes $T$ a positive integer.)
Thank you

Your assumption that the fundamental period has to be an integer since we are in discrete time is wrong. It can be any real number, nobody says, there has to be an integer number of periods in a signal of finite length. Also, if there is no $k$ representing this exact period, the peak in the spectrum will smear across neighboring bins.


Overlap/add time-domain audio frames: How does normalization/scaling work with overlap greater than 50%? I am attempting to overlap-add time-domain audio frames in a real-time context (with a continuous stream of data coming-in). Each frames consist of 1024 samples and have been previously "Hann" windowed in view of performing FFT/IFFT operations.

When the hopsize is set to exactly 50% (here, 512 samples), everything works fine as the frames "concatenate" in a way that it creates a flat-top. Consequently, there isn't any amplitude build-up occurring:

However, as soon as I try to set the overlap to a greater value, let's say 75% (=768 samples), the frames expectedly stack-up and their respective amplitudes accumulate in a troublesome way:

Unfortunately, I couldn't find any clear documentation on how to properly normalize overlapping frames based off the hopsize value, the framesize and the window type.

What would be the proper way of scaling my audio signal ?

Thanks in advance

Let's assume continuous time (rather than discrete time).

If you do not process the windowed data at all, you would like the output (the sum of the windowed frames) to be equal to the original signal. Allowing scaling of the output by a constant scaling term, this is only possible if the sum of all of the time-shifted window functions is constant over time.

A typical window function (including Hann) is a sum of harmonic cosine terms. Let's assume you only use such windows. At the end of the window, the cosine terms have derivative discontinuities that can only be "cured" by stitching to the end of the window the start of another window which has identical values of all derivatives. This, and that you probably want the windows to be distributed evenly, requires that the hopsize is $\frac{1}{N}$ window lengths, with $N$ integer. The truncated cosine terms in each of the first $N$ windows are extended to complete cosines by concatenation of every $N\text{th}$ window.

The sum $f(T)$ of the windows distributed this way is:

$$f(T) = \sum_{n = 0}^{N - 1}\sum_{h = 0}^{H - 1} a_h\cos\left(2\pi h\left(T-\frac{n}{N}\right)\right),$$

where $h$ is the harmonic number, $H$ is the number of harmonics in the window function, $a_h$ are the window function specific coefficients, $T$ is time expressed so that a time difference of $1$ means one window length. $h$ and $n$ are the running harmonic and hop numbers.

Let's try window densities $N = 1\dots3$ with a window that has two harmonics ($H=2$). First $N=1:$

$$f(T) = a_0 + a_1\cos(2\pi T)$$

then $N=2:$

$$f(T) = a_0 + a_1\cos(2\pi T) + a_0 + a_1\cos\left(2\pi \left(T - \frac{1}{2}\right)\right) = 2a_0$$

and $N=3:$

$$f(T) = a_0 + a_1\cos(2\pi T) + a_0 + a_1\cos\left(2\pi \left(T - \frac{1}{3}\right)\right) + a_0 + a_1\cos\left(2\pi \left(T - \frac{2}{3}\right)\right) = 3a_0$$

This could be for example the Hann window with $a_0 = \frac{1}{2}$ and $a_1 = \frac{1}{2}:$

Figure 1. Hann windows sum to constant with a hopsize of 1/2 or 1/3 window lengths, but not with a hopsize of 1 window length.

Although it is not as evident as with the 1/2 hopsize, the $1\text{st}$ harmonic cosines cancel each other completely with the 1/3 hopsize. The cancellation can be illustrated by representing the phases and amplitudes of the harmonic terms as complex numbers:

Figure 2. Complex number representation of the amplitudes and phases of the 1st harmonic terms in shifted concatenated Hann windows, with a hopsize of $1/3$ window lengths. The rotational symmetry of the pattern means that the first harmonics cancel, because the sum of the three complex numbers must be zero.

For window functions with $H$ harmonic cosine terms, if we tabulate for different window densities $N$ the factors appearing in front of each $h\text{th}$ harmonic cosine term in $f(T),$ a pattern emerges:

$$\begin{array}{r|lllllllllll}&\rlap{h}\\
N&0&1&2&3&4&5&6&7&8&9\\
\hline
1&1&1&1&1&1&1&1&1&1&1\\
2&2&0&2&0&2&0&2&0&2&0\\
3&3&0&0&3&0&0&3&0&0&3\\
4&4&0&0&0&4&0&0&0&4&0\\
5&5&0&0&0&0&5&0&0&0&0\\
6&6&0&0&0&0&0&6&0&0&0\\
7&7&0&0&0&0&0&0&7&0&0\\
8&8&0&0&0&0&0&0&0&8&0\\
9&9&0&0&0&0&0&0&0&0&9\end{array}$$

The factors in the 2-term window example appear in the two first columns and three first rows.

For $H \le N,$ all harmonics except for the $0\text{th}$ harmonic cancel completely. This means for example that the 4-cosine-term ($H = 4$) Blackman–Harris window could be used with a 1/4 hopsize or greater ($N \ge 4$), requiring only an output scaling by $\frac{1}{4a_0} = \frac{1}{4\times0.35876}.$


Can I compute my IIR faster given limited possible input values? I have an IIR routine in C++ I'm trying to speed up.

Currently I'm using a DF2 implementation. Since my input signal will always be either 1.0 or -1.0 is it possible for me to optimize the algorithm, removing some calculations or precomputing values for given blocks of input?

Thanks!

Depends on the hardware running your C++ code.  If you have a slow multiply unit and fast branching, then you might be able to replace a few slower multiply-adds with faster adds or subtracts.  But on most modern GHz CPUs, the add/sub branch penalty will likely be slower than a floating-point mul-add by +-1.0.


What is the need of IIR filter? When FIR have generally good/desirable characteristics? FIR filters are generally stable and have linear phase response. These two characteristics are often desired in dsp

So why do we need and use IIR filters,despite the fact that IIR filters often have issues in these two characteristics of stability and linear phase response ?

The biggest advantage of IIR filter is that you can achieve the filter specifications with much lesser order compared to FIR.
In a recent question (Lowpass Hann filter in Python) a transition band $\omega_{cutoff-passband}=0.05\pi$ to $\omega_{stopband}=0.136\pi$ was achieved with just $N_{IIR}=4$, but using FIR filter with Kaiser window, the order $N_{FIR}=55$. This results in much less complexity in implementation (multipliers, adders) reducing area required for implementation on a chip. In the above exmaple, for a given output, IIR requires 4 multipliers and 5 addition operations. For FIR, 55 multiplications and 54 addition operations are needed. Of course, there are provisions for optimizing this, but you get the idea..
Also, for the given example, in the pass band the phase is nearly linear for the IIR filter, so it results in much less distorted output for the pass band frequencies (stop band frequencies will get distorted but they are attenuated away by the filter). In IIR, you have the additional flexibility of placing poles as per your requirement to fine tune the filter attenuation/ripple. So it may seems FIR has advantage in terms of pure linear phase but IIR isn't far behind.
Code used to compare the order
clc
clear all

fc1=1200;
fc2=3000;
fs=44100;
rp=1;
att=40;
[n,wp]=cheb1ord(fc1/(fs/2), fc2/(fs/2), rp, att);
[b,a]=cheby1(n,rp,wp,'low');
freqz(b,a)
[n1,Wn,beta,ftype] = kaiserord([fc1 fc2],[1 0],[0.9 0.01],fs);
hh = fir1(n1,Wn,ftype,kaiser(n1+1,beta),'noscale');
freqz(hh,1,1024,fs)

Zoomed view near transition band

As you can see Chebyshev Type-1 provides a much higher attenuation of frequencies beyond pass band cutoff ($\text{1200Hz}$).


Is digital up and down sampling linear/causal/time-invariant? So I am trying to determine whether the process of up sampling (zero padding) and downsampling are linear, causal and/or time-invariant.

Based on some resources I found online, I am under the impression that it is linear, but not time-invariant. However I am unsure about causality (since it depends neither on future nor past values?).

Upsampling and downsampling are both causal, since the output at time $m$ depends only on samples you've seen up time $m$.

It is also linear (apply upsampling/downsampling to a x[n]+y[n] and see that you get a upsample/downsample x[n] + upsample/downsample y[n]).

Downsampling not time invariant though - what happens if you downsample by 2 and consider inputting x[n] and x[n] delayed by 1 time unit? Now, what about upsampling?


normalized sum of squared differences I want to search for a patch in a region of image, and I will use the normalized sum of squared differences.

I know that the sum of the squared differences is:

$SSD = \sum {\left(F - I\right)^2}$

Where:
$F$: is the reference patch
and
$I$ is the patch from the image to match with $F$

My question is:  what does it means normalized? is it the same as normalized correlation?

In an attempt to solve the question on why to normalize, and implicitly how to normalize:

$F$ is your reference patch, $I$ is the patch under inspection.

So make $I$ only consist of the maximum possible value in your image format. Then, your $SSD$ becomes very large, although $I$ isn't "similar" to $F$ at all.

Now, for comparison, set $I=F$. In theory, this should give a large (if making any sense at all, the largest) $SSD$ possible, right?

So, obviously, your $SSD$ formula as is isn't useful at all for comparing things, because an $I$ with a high average amplitude will just come out as a "winner", always. So, you need to get the energy of the $I$ out of your formula first. That's easy: Just divide $I$ by the sum of the squares of all its values.


Looking for cycles of periods longer than the input signal length There was a similar question asked here, however I would like to focus on some specifics of this problem. Let me present the Python code sample to illustrate the situation when only a part of the period of the input signal is analyzed using FFT:

N = 1200
H = 4
fig = plt.figure()
sig = np.zeros(N)
for h in range(H):
        sig += [(h+1)*np.sin(2*np.pi*(h+1)*x/N) for x in (range(N))]
ax1, ax2 = (None, None)
for h in range(H):
        if (h):
                fig.add_subplot(H, 2, 2*h+1, sharex=ax1)
        else:
                ax1 = fig.add_subplot(H, 2, 2*h+1)
        plt.plot(sig[:int(N/(h+1))])
        if (h):
                fig.add_subplot(H, 2, 2*h+2, sharex=ax2)
        else:
                ax2 = fig.add_subplot(H, 2, 2*h+2)
        fft = pd.DataFrame()
        fft['coef'] = np.fft.fft(sig[:int(N/(h+1))])
        fft['freq'] = np.fft.fftfreq(int(N/(h+1)))
        fft['per']  = 1/fft['freq']
        fft['mag']  = fft['coef'].apply(lambda z: np.sqrt(z.real**2 + z.imag**2))
        rside = fft['freq'] > 0
        plt.bar(fft[rside]['per'], fft[rside]['mag'])
fig.show()

The parameter H specifies the number of composite sines, their magnitude and the number of divisions of the original signal. For each division we get 1 row of subplots: the input signal in the left column and the period spectrum in the right column. The X axes in each column are shared, so the plots are aligned. Only positive frequencies are considered. The example with H = 4 and the signal length at N = 1200 should produce the set of plots like on the picture below.

The FFT results in the top row are spot on, but on other plots from the right column we can see that the spikes representing the correct dominant cycles are missing, but some additional smaller spikes appear at periods much smaller than the smallest one used to generate the test signal.

In a general case when given a sample of discrete observations or a time series data I would like to have a way to test how would some bigger cycles fit in the spectrum, ie. if we consider the Mercury's distance from the barycenter in that other thread someone suggested that orbital cycles of other planets should be considered, all of which are well known and some are longer than 8 years (the size of the input data). I was searching online, but due to lacks in theoretical (and practical) knowledge on the matter I probably wasn't using the right keywords, so I couldn't find any helpful information. These are the questions I would like to ask:

Can the small cycles that FFT 'found' in the period slices provide clues about the bigger ones that were actually in effect, ie. could there be a correlation in frequency or phase between the small fake cycles and the real big ones?
Is there a way to expand the period domain or in other words assume there are bigger cycles at play and include them in the spectrum? Can Fourier Transform be somehow manipulated to do so?
Lastly is Fourier Transform the right tool or at least a starting point to deal with such issue?

If you plot frequency instead of period, you might be able to see that some of the increasing artifacts are due to the shrinking rectangular window on your signal.

An FFT decomposes the entire length of its input into sinusoids.  You have created a flat spot in your data input, and it takes a linear combination of a lot of sinusoids to represent those extended flat spots.  So the FFT of your latter shorter signal segments becomes more and more about the flat spot than your signal.

Added: There are other methods for estimating the parameters of a known small number of sinusoids within a known distribution of noise, such as MUSIC and ESPIRIT, that may be more suitable for estimation of sine waves not equal to or possibly shorter than the lowest DFT basis vectors.


Books on medical imaging As the title suggests , I'm looking for good books on medical imaging . It would be a great help if you can provide a link to the books website .

Thanks.

I quite like the Medical Imaging Signals and Systems by Jerry L. Prince and Jonathan Links. If you're looking for a DSP approach to medical imaging, this is a good start. It's also quite good on the physics of medical imaging. the format of the book follows one medical imaging modality at a time. It first teaches you the relevant physics of how the imaging itself works for that modality, and then goes into the signal processing side of things in the following chapter.


A good mathematical explanation of Gibbs phenomenon I was explaining to someone how Fourier series work in context of constructing signals that are not everywhere differentiable, e.g. square waves, sawtooth waves, etc. When I mentioned the Gibbs phenomenon however, I realized that I never really learned of why it happens. In fact, as the story goes, not everyone even realized that it's an actual mathematical property of infinite series of periodic signals and not a computational fluke, and it turns out that most proofs are fairly laborious and elaborate.

After reading several of them, I started realizing why such phenomenon might occur, but I have a background in real and complex analysis, topology and so on. The question is can I fully explain and rigorously prove Gibbs phenomenon mathematically to someone with only the basic undergraduate calculus courses in their arsenal (or any other general prerequisites for an undergraduate signal processing course)? If so, then how?

The book "Dr. Euler's Fabulous Formula: Cures Many Mathematical Ills", by P. Nahin, Princeton University Press, leads up to and contains an explanation of the Gibbs phenomena which might be suitable for someone with a good undergraduate university level math background.


Why does the Hilbert-Huang transform not suffer from the uncertainty principle, and what are the implications of this? I have been studying the Hilbert-Huang transform, which I understand is a frequency analysis technique that uses the adaptive IMFs basis components to decompose a signal using EMD.

I understand that frequency analysis techniques such as the STFT (short-time Fourier transform) and Wavelet transform involve the convolution of a signal with a predefined basis in order to extract frequency information.

Several articles from Huang and others state "convolution processes involve integration, which makes the results suffering the limitation imposed by the uncertainty principle".

http://web-static-aws.seas.harvard.edu/climate/pdf/Zhaohua.pdf
https://pyhht.readthedocs.io/en/latest/tutorials/hilbert_view_nonlinearity.html

I also understand that the frequency of each IMF is extracted by finding the time rate of change of phase of the analytic signal for each IMF (which would involve differentiation of the phase with respect to time).

Why does this differentiation mean that the uncertainty principle is avoided?

I have a vague idea that it's due to operators associated with time and frequency that don't commute, but I have struggled to find any sources that approach this from a signal processing/mathematical perceptive rather than in a quantum mechanical sense.

Furthermore, if this kind of signal processing circumvents the uncertainty principle, does that mean time and frequency events can be resolved with arbitrary precision?

(In my head, this, in turn, would have implications in quantum mechanics, but I imagine this would be outside the topic of this site.)

Any insights and/or articles are gratefully received.

I'm not that familiar with the Hilbert-Huang transform, so I won't comment on that.

You seem to be under the impression that the DFT "suffers from the uncertainty principle".  This is not true.  Don't feel bad about thinking that as it is a widespread misunderstanding.  This is a quote from an email to me from a well known expert in the field after I tried to get them to look at my exact frequency formula for a single (non-integer) pure real tone in a DFT:

  An "exact" computation of frequency based on FTs is not possible for
  some very simple reasons:  The Fourier Transform is itself an
  estimator, and the time-frequency uncertainty principle means
  frequency can't be determined exactly without an infinitely long
  observation window.  So essentially all frequency computations of time
  series or signals are estimates, it's just that some are more accurate
  or more efficient than others.

  So claiming that you have an exact solution just isn't going to get
  the attention of everyone interested in the field.

The first paragraph is entirely hogwash.  The last sentence was prophetic, hence I started my blog so no one could dispute the math.  You will find my original exact solution here:

Exact Frequency Formula for a Pure Real Tone in a DFT

I have discovered, and written about, several more.  You will find them among my articles.  There are also a few in the pipeline, including a correction to Macleod's three bin formula which makes it exact also.

The FT does apply to the uncertainty principle in physics, but in a different manner.  There they are talking about the transform of a Gaussian curve which happens to be an eigenfunction of the FT (but not the DFT, close though) and how it is the "narrowest function" in the transform.

I'm not an expert in Physics either, but I dabble.  I consider this to be my greatest mathematical discovery ever (exclipses the simple math of the frequency formulas by a mile):

Off Topic: Refraction in a Varying Medium

This is the vector form of Snell's law in a continuously varying index of refraction medium.  It is a theoretical particle model so its applicability to light is still a subject of study for me.  I believe it has the potential to illuminate subatomic behavior (i.e. the quantum world) but I haven't gotten there yet.


Quantization processs In a quantization process what happens when the signal range is more than the quantizer range?
Suppose i have a sinusoidal signal which ranges from -5 to +5 but the quantizer which has four levels ranges from -3 to 3.now for the samples which are greater than 3 and less than -3 what will happen?

In a typical quantization like that seen in an ADC the output would saturate any value > 3 to 3 (also known as clipping). Note that if you're considering a software quantization, it's also possible that some quantizers (although IMO poorly designed ones) could cause out of range values to wrap around and switch signs.


Programming the IDWT for image processing I want to program the 2D inverse discrete wavelet transform (only 1 level) in the case of image processing. In the matlab website there's this diagram:

now, I want to program the IDWT with haar wavelets. What are the kernels of the low pass filter and the high pass filters?

here is one version of the Haar analysis/synthesis filter bank pair. Oddly, the synthesis filters, up to a scale factor, are mirrors of  the analysis filters.

Other options are:

$\pm 1/\sqrt{2}$ instead or $1/2$ and $1$ (orthogonal version);
reverse: $[1 \;  \pm 1]$ at analysis and $[\pm 1/2 \; 1/2]$ at the synthesis.


Why is my DFT/FFT always 90° out of phase? I'm doing an FFT using Python and Numpy on one machine, and C# on another. I'm using some dummy data that mimics how I'll eventually be gathering data from sensors in the C#/UWP application. The two methods yield consistently identical results on the dummy data, which is great. However, the phase at the sine wave's frequency is always -90°.

Here's the code. Ignore how it essentially does the same thing 7 times, eventually those arrays will be filled with data from 7 different sensors. Sampling frequency is 4096 Hz, signals are gathered for 1 second and are on the interval [-32768,32768]. Note that the x-axes on the Magnitude and Phase graphs are zoomed in to the center frequency plus or minus several Hertz.

import numpy as np
import matplotlib.pyplot as plt
import math

centerf = 1350
span = 2*np.pi*centerf*np.linspace(0,4095,4096)/4096
testData = 32768*np.sin(span)
testData = np.array([math.trunc(x) for x in testData])

plt.close("all")
sensors = ["Channel 1", "Channel 2", "Channel 3", "Channel 4", "Channel 5", "Channel 6", "Channel 7"]

data = []
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)

theTime = np.linspace(0,len(data[0]),len(data[0]))/len(data[0])

thisFFT = [np.fft.fft(x) for x in data]

f, a = plt.subplots(7, 1,figsize=(10,8),sharex=True)
f2,a2 = plt.subplots(7,1,figsize=(10,8),sharex=True)
f3,a3 = plt.subplots(7,1,figsize=(10,8),sharex=True)

mag = [np.abs(x)/2048 for x in thisFFT]
phase2 = [np.arctan2(x.imag,x.real)*180/np.pi for x in thisFFT]

mag = [x[:2048] for x in mag] # removed mirrored upper half
phase2 = [x[:2048] for x in phase2]

for x in mag:
    x[0] = 0 # remove DC component

[aa.plot(theTime,d,c='r',lw=0.5) for (aa,d) in zip(a,data)]
[aa.plot(d,c='g',lw=0.75) for (aa,d) in zip(a2,mag)]
[aa.plot(d,c='b',lw=0.75) for (aa,d) in zip(a3,phase2)]
for ch,ax in zip(sensors,a.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a2.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a3.flat):
    ax.set(ylabel=ch)
a.flat[0].set(title="Amplitude over 1s")
a2.flat[0].set(title="FFT Magnitude")
a3.flat[0].set(title="FFT Phase")
a.flat[5].set(xlabel='Time')

[aa.axvline(x=centerf,lw=0.25) for aa in a2]
[aa.axvline(x=centerf,lw=0.25) for aa in a3]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a2]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a3]

f.tight_layout()
f2.tight_layout()
f3.tight_layout()

print("Magnitude at centerf: {c:1.3f}.".format(c=mag[0][centerf]))
print("Phase at centerf: {c:1.3f}°.".format(c=phase2[0][centerf]))

The output from the print statements is

Magnitude at centerf: 32767.397.
Phase at centerf: -90.000°.

The magnitude is fine, the graph always shows a nice peak right where it should be.

My question: since the input data is a sine wave that always starts at sin(0)=0 and the FFT algorithm uses a rectangular window, shouldn't the phase at the fundamental frequency be 0°? Is there something about phase that I'm not understanding?

The Fourier Transform (and the DFT is the discrete form of this with a finite length sequence) is simply a correlation of the input signal with each frequency bin, of the continuous form:

$$X(\omega) = \int_{t=0}^Tx(t)e^{-j\omega t}dt$$

(If we consider a finite time duration, however I will stay with the continuous time version of the Fourier Transform as it will simplify the explanation).

Here we see that for any freqency $\omega$ we are doing a correlation of $x(t)$ with $e^{j\omega t}$  (correlation is a complex conjugate multiplication followed by integration, or summation when discrete).

As a correlation we see that the result will be maximum when x(t) is fully correlated to $e^{j\omega t}$ over the given time interval, and similarly 0 when it is fully uncorrelated (orthogonal).

Now make note of Euler's identity given as:

$$e^{-j \omega t} = cos(\omega t) - j sin(\omega t)$$

And replacing that in our first formula we can expand the Fourier Transform given into two correlations:

$$X(\omega) = \int_{t=0}^T x(t)cos(\omega t)dt - j\int_{t=0}^T x(t)sin(\omega t)dt$$

From this we see that when $x(t) =sin(\omega t)$, (and duration T covers one full cycle or an integer number of cycles), the second term will be fully correlated, while the first term will be zero, and vice versa when $x(t) =cos(\omega t)$.  Any correlation to the second term results in an imaginary result IF x(t) is real.  Any phase shift on x(t), as in $x(t) = sin(\omega t + \phi)$ where $\phi$ is anything except $\pm N\pi/2$ for N any integer, will result in a complex result with $X(\omega)$ having real and imaginary components.

Note the discrete form of this, the DFT, is as given below where the same result can be explained, but may not be as easy to see for some that are less familiar with it.

$$X[k] = \Sigma_{n=0}^{N-1}x[n]e^{-jn\omega_o \frac{k}{N}}$$


Generate noise with arbitrary covariance I am doing a time series change detection experiment.  One of the experiments includes generating some artificial data for model evaluation.

I came across a problem in generating noise data with some arbitrary covariance at every 100 time steps.

While randn from Matlab allows me to generate the value, I am never able to figure how to change the covariance of the data.  Is there any example that I can refer to?

If you want a Gaussian vector $X$ with $E(X)=0$ and $E(X\,X^T)=C$, then find a decomposition $C=AA^T$, for instance a Cholesky decomposition, compute a vector $U$ with standard normal uncorrelated random numbers as components and set $X=AU$.


Understanding FFT for a vibrational analysis experiment I am trying to understand how FFT works in order to measure the vibrations of a 3D printer. I will be doing this with an Arduino Nano and a MPU6050 accelerometer.
In order to program the Arduino correctly I need to understand some things about FFT first. So various sources state that FFT works with 2^n samples, for example "The FFT-algorithm works with a finite number of samples. This number needs to be 2^n where n is an integer (resulting in 32, 64, 128, etc)". What does this actually mean? From my understanding, if I will be measuring accelerations at a constant sample rate of 600Hz for 10 seconds, I will end up with 6000 samples. Does this mean that I will have to remove some of the measured data in order to end up with 4096 samples, which meets the 2^n samples condition? Or fill up my samples with zeroes to end up with 8196 samples, which also meets the stated condition?
Thank you!

The following is written assuming you do the signal analysis on your PC, and only use the Arduino to acquire the samples and send them to your PC.
Your Arduino Nano's microcontroller doesn't have the memory necessary to store a second of data.

various sources state that FFT works with 2^n samples

FFTs in length of powers of 2 are especially computationally efficient to compute. Since you probably don't care about how efficient your computation is, you don't care. The FFT is just a fast way to calculate the DFT of a vector, and worst case, any length is possible when just calculating the DFT plain from the textbook formula for the DFT.

This number needs to be 2^n

That's plain wrong. FFTs exist for many lengths, usually for any length that can be decomposed into prime factors < 17. (As far as I know. Maybe we can do better today.)

From my understanding, if I will be measuring accelerations at a constant sample rate of 600Hz for 10 seconds, I will end up with 6000 samples.

That sounds right.

Does this mean that I will have to remove some of the measured data in order to end up with 4096 samples, which meets the 2^n samples condition? Or fill up my samples with zeroes to end up with 8196 samples, which also meets the stated condition?

you can do either, but you don't have to; 6000 = 2·3·10³ = 2⁴·3·5³ and there's multiple FFT implementations for such lengths. Again, you don't even have to use an FFT to calculate the DFT.
So, I'd really rather adjust the FFT (DFT) length to your measurement problem than the other way around.


Frequency response of $\mathrm{sinc}[n]$ In this image the frequency response of a discrete time filter given as $h[n]$. Can someone explain how the magnitude of the frequency response is found ?

First of all, there is a mistake in that image: there should be an $n$ in the denominator instead of a $t$.

Check that your signal can be expressed as

$$h[n]=2\frac{\sin\left(\frac{\pi}{2}n\right)}{\pi n}=\frac{\sin\left(\frac{\pi}{2}n\right)}{\frac{\pi}{2}n} \tag{1}$$.

Using Euler's formula we can express that sine as the sum of two exponentials:

$$h[n]=\frac{1}{\pi/2}\frac12\frac{e^{jn\pi/2}-e^{-jn\pi/2}}{jn}=\frac{1}{\pi}\frac{e^{jn\pi/2}-e^{-jn\pi/2}}{jn}$$

Let me multiply and divide that expression by $2$, which doesn't change anything but will be useful later:

$$h[n]=\frac{1}{2\pi}2\frac{e^{jn\pi/2}-e^{-jn\pi/2}}{jn} \tag{2}$$

Notice that there exists an integral that returns that exact function of $n$ and, on top of it, that integral represents the inverse Fourier transform of some function:

$$\frac{1}{2\pi}\int\limits_{-\pi/2}^{\pi/2}2e^{j\omega n} \ \mathrm{d}\omega=\mathcal{F}^{-1}[H(e^{j\omega})]=h[n] \tag{3}$$

where $H(e^{j\omega})$ is periodic in $2\pi$ and, for each period:

$$H(e^{j\omega}) =
\left\{
 \begin{array}{ll}
  2  & \mbox{if } |\omega| \leq \pi/2 \\
  0 &   \ \mathrm{otherwise}
 \end{array}
\right.$$

So the Fourier transform of $h[n]$ is $H(e^{j\omega})$, because the inverse Fourier transform of the latter is the former. And that's the proof.

However, most of the times one does not do all these calculations and just looks up the transform of $\mathrm(1)$. (In the link's table, look up for your signal and replace $W=2$. Then use the linearity of the Fourier transform to get to the result you are looking for.)


Vehicle detection in matlab I am new to matlab and computer vision. I used the code given in the link to detect cars. But the problem is the shadow is also being detected and tracked. Can anyone tell me how to modify the code so that only car are detected. Link:http://in.mathworks.com/help/vision/examples/detecting-cars-using-gaussian-mixture-models.html
Thanks in advance

What I would try is to :

Get the HSL (Hue Saturation Lightness) decomposition of the original picture
Only keep the L component
Multiply the image you have from the last steps with the image you get after the morphological operation (the black and white one, after deleting the noise)
Do some thresholding on that
Use what you obtained from that instead of the image after the morphological operation for the following part of the algorithm.

I don't know well Matlab so I can't provide function names, but I guess you can find it easily with a bit of googling. And please keep me updated if you try this, I would like to know if it works well or not.


Radar Imaging Method I've been working with an FMCW Fast Chirp radar transceiver.  Mainly, I've been looking at the Doppler range map but really it is not suitable as there is only one point per range/ speed and I want to build a better picture of the surroundings.

I came across the following article.  The below picture struck me as excellent:

Does anyone know how I could generate a picture as shown on the right using radar?

You still only get one range and one velocity line per radar observation.

You can, however,  "scan" a region (that's why airport towers and ships have these rotating antennas) and get another axis of information.

The same can be achieved using digital beamforming, which, however, requires multiple signal chains, and specific antenna arrays. For a two-dimensional picture, you'd need an at least two-dimensional array of antennas. For the resolution displayed, you'd need both high bandwidth (think in the order of 1GHz) as well as small wavelength (think of beyond 76 GHz). So, I'd say, without a millimeter-wave lab, this is probably not achievable.

I've not seen pictures like that myself (the article says "source: NXP" without elaborating any further, so this might just as well be "artist's interpretation of what NXP wants to achieve in 10 years).

Whilst this resolution might be impossible to achieve with commodity hardware, in principle, SDRs that are MIMO-capable do exist. If you can get your hand on a 4x4 MIMO phase-aligned SDR system, you could, behind/infront ambiguities notwithstanding, build a simple patch antenna array of 2 rows / 2 columns parallel patches and build a digital beamforming radar with that. You'd probably want to avoid the problems of isolating RX from TX by having a separate array for TX and RX.


Covariance Matrix Polar to Cartesian Let $R$ and $\Theta$  be zero mean independent Gaussian random variables with variances $\sigma_R^2$ and $\sigma_\Theta^2$ respectively. Hence their covariance matrix would be diagonal give by
$$C_{R,\Theta}=\begin{bmatrix}\sigma_R^2 & 0 \\ 0 & \sigma_\Theta^2\end{bmatrix}$$
Now, let's define transformations given by $X = R\cos (\Theta)$ and $Y = R \sin (\Theta)$.

The question is, what would be $C_{X,Y}$, the covariance matrix of transformed random variable $X$ and $Y$?

Thank you in advance.

The covariance matrix is given by

$$C_{X,Y}=\begin{bmatrix}E(XX)& E(XY) \\ E(YX )& E(YY) \end{bmatrix}$$

This can be written as below:

$$C_{X,Y}=\begin{bmatrix}E(R^2cos^2(\Theta) )& E(Rcos(\Theta)Rsin(\Theta)) \\ E(Rsin(\Theta)Rcos(\Theta) )& E(R^2sin^2(\Theta)) \end{bmatrix}$$

Since $R$ and $\Theta$ are independent the expectation will distribute. Now the task is to calculate expectation of $Cos^2(\Theta)$
, $Sin^2(\Theta)$ and $Cos(\Theta)Sin(\Theta)$ or equivalently $Sin(2\Theta)$. These can be found by by calculating the expectations of each of these functions for ex: for $Cos^2(\Theta)$ this will be $$E(cos^2(\Theta)) = \int_{ - \infty}^{\infty}\, Cos^2(\theta) f_\Theta(\theta) \, d\theta$$

Where $f_{\Theta}(\theta)$ is the PDF of random varible $\Theta$

Hint: you can also utilize symmetries while calculating these integrals, $Sin(2\Theta)$ is an odd function and it multiples a Gaussian with zero mean which is an even function. Similarly $Cos^2(\Theta)$ and $Sin^2(\Theta)$ are even functions

EDIT, Solving the integral:

$$\cos^2\theta=0.5+0.5\cos 2\theta=0.5+0.25(e^{2i\theta}+e^{-2i\theta})$$

$$I=\int_{-\infty}^{\infty} \cos^{2}(\theta) f_{\Theta}(\theta) d\theta=0.5+0.25\int_{-\infty}^\infty {1\over \sqrt{2\pi}} e^{-\theta^2\over 2}(e^{2i\theta}+e^{-2i\theta})d\theta$$

$$\int_{-\infty}^\infty e^{-\theta^2\over 2}e^{ki\theta}d\theta {= \int_{-\infty}^\infty e^{-\theta^2\over 2}e^{ki\theta}e^{-{(ki)^2\over 2}}e^{{(ki)^2\over 2}}d\theta \\= \int_{-\infty}^\infty e^{-(\theta-ki)^2\over 2}e^{-{k^2\over 2}}d\theta \\= \int_{-\infty}^\infty e^{-\theta^2\over 2}e^{-{k^2\over 2}}d\theta \\= \sqrt{2\pi}e^{-{k^2\over 2}} }$$ hence $$
I=0.5+0.5e^{-2}$$

Original answer to the integral  https://math.stackexchange.com/questions/3650453/integral-of-int-infty-infty-cos2-theta-f-theta-theta-d/3650494#3650494