Digital carrier recovery I was reading Chapter 2 from "Telecommunication Breakdown" by Johnson and Sethares and came across this paragraph:

  However, sometimes it is more cost effective to perform certain tasks
  in analog circuitry. For example, if the transmitter modulates to a
  very high frequency, then it may cost too much to sample fast enough.
  Currently, it is common practice to perform some frequency translation
  and some out-of-band signal reduction in the analog portion of the
  receiver. Sometimes the analog portion may translate the received
  signal all the way back to baseband. Other times, the analog portion
  translates to some intermediate frequency, and then the digital
  portion finishes the translation.

And so I was wondering, if I wanted to do carrier recovery digitally, where in the receiver chain would I perform such task? I guess I'd need to translate input RF signal to IF and then perform carrier recovery digitally but in order to do so, I need to know the frequency and phase of the RF signal to translate it to IF, don't I? And I can only do that by carrier sync/recovery.

You don't need to know the phase at all and you only need to know the frequency approximately. It is the job of your digital carrier recovery scheme to estimate and remove the phase and residual frequency offset.


Filtering and Fourier Transforming, does the order matter? I have a signal $x(t)$. I want to find the Fourier Transform of it, $X(f)$, and then extract a narrow frequency range from $X(f)$ by use of a Band Pass Filter (BPF) in frequency domain.

Can I instead filter $x(t)$ by using a BPF in time domain and then find the Fourier Transform of the filtered signal?

I believe these two are equivalent.

There may be some slight differences due to the band-pass filtering of FFT results in the frequency domain being a circular convolution (with some wrap-around artifacts) rather than a pure linear convolution.

If you filter first, starting in time before your FFT window, any windowing artifacts from any out-of-band spectrum will be reduced before leaking into the FFT filter's pass-band.  Or alternatively, you could zero-pad the FFT window by the length of your filter's impulse response.


Inverse Continuous Wavelet Transform derivation? Wiki writes iCWT as
$$
f(t) = C_{\psi}^{-1} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} W_f(a,b) \frac{1}{|a|^{1/2}} \tilde\psi \left(\frac{t - b}{a}\right) db \frac{da}{a^2}, \tag{1}
$$
where $\tilde\psi$ is the dual wavelet of $\psi$ in CWT, $W_f$ is the CWT of $f$, and $C_\psi$ is the admissible constant
$$
C_\psi = \int_{-\infty}^{\infty} \frac{\bar{\hat{\psi}}(\omega) \hat{\tilde{\psi}}(\omega)}{|\omega|} d\omega, \tag{2}
$$
Hat = $\mathcal{F}$, overbar = complex conjugate. Some wavelets are own duals, $\psi = \tilde\psi$; then the integrand is $|\hat\psi(\omega)|^2 / |\omega|$.

How does $(1)$ invert the CWT? And why is there a $a^{-2}$ normalizing factor? Above wavelets are L2-normed - do we normalize iCWT differently if they are L1-normed (i.e. $|a|^{-1}$)?
The normalization is particularly troubling; in forward-transform, we do $|a|^{-1/2}$, so in inversion one might expect $|a|^{1/2}$ for cancellation, but it only goes the other way: another $|a|^{-1/2}$, along $a^{-2}$. The role of $C_\psi ^{-1}$ is also unclear.

Summary: the dual wavelet's role is analogous to that of $e^{j\omega t}$; it undoes the wavelet's convolving with the signal (integrated inner product). The main intricacy's indeed in normalization; there's no exact answer without digging in math (done below) - but intuitively, the seemingly matter-worsening $a^{-2}$ is due partly to the dilating nature of the wavelet, i.e. $\psi(t/a)$, which upon inversion requires downscaling not only to properly preserve energy, but reconstruct correctly by undoing the inner product.
Below are more "proofs" rather than "derivations", as they work backwards to show equality with $f(t)$, but I term them latter as full-fledged proofs are beyond the scope of this answer.

Derivation 1: Convolution Theorem
Pasted for convenience the CWT and iCWT
$$
W_f(a, b) =   \int_{-\infty}^{\infty} f(t) \frac{1}{|a|^{1/2}} \psi \left(\frac{t - b}{a}\right) dt \tag{0}
$$
$$
f(t) = \frac{1}{C_{\psi}} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} W_f(a,b) \frac{1}{|a|^{1/2}} \psi \left(\frac{t - b}{a}\right) db \frac{da}{a^2}, \tag{1}
$$
Begin by rewriting the CWT as a convolution: $W_f(a, b) = f \star {\bar\psi_a}(b)$, where $\psi_a(t)=a^{-1/2}\psi(t/a)$. Inserting into $(1)$, and letting the right side equal $b(t)$,
$$
\begin{align}
b(t) &= \frac{1}{C_\psi} \int_0^{\infty} W_f(.,a) \star \psi_a(t) \frac{da}{a^2} \\
     &= \frac{1}{C_\psi} \int_0^{\infty} f \star {\bar\psi_a}(b) \star \psi_a(t) \frac{da}{a^2} \tag{2}
\end{align}
$$
where we also expressed the $db$ integral as a convolution, and $"."$ indicates the variable over which convolution is computed. We prove $b=f$ by showing their Fourier Transforms are equal. With $\hat{} = \mathcal{F}$,
$$
\begin{align}
{\hat b}(\omega) &= \frac{1}{C_\psi} \int_0^{\infty}
    {\hat f}(\omega) \sqrt{a} {\bar{\hat \psi}}(a\omega) \sqrt{a} {\hat\psi}(a\omega) \frac{da}{a^2} \tag{3a} \\
                 &= \frac{{\hat f}(\omega)}{C_\psi} \int_0^{\infty} \tag{3b}
    |{\hat \psi}(a\omega)|^2 \frac{da}{a}
\end{align}
$$
With change of variables $\xi = a\omega$, we thus 'prove'
$$
{\hat b}(\omega) = {\hat f}(\omega) \left(\frac{1}{C_\psi} \int_0^\infty \frac{|{\hat \psi}(\xi)|^2}{\xi} d\xi \right) = {\hat f}(\omega). \tag{4}
$$
Interpretation:
It's tempting to attribute the $a^{-2}$ to the wavelet dilations, $\psi(t/a)$, as they manifest as rescalings of Frequency-domain components ($\cdot \sqrt{a}$), and while that's true, it's not the crux.
Recall when I said it's not just about energy? Consider what happens if the normalization is anything other than $a^{-2}$. To arrive at $(4)$, we relied on $a^1$ in denominator in $(3b)$, else there is no cancellation with $C_\psi$, so no reconstruction. Thus, one can conclude that $a^{-2}$ is set in retrospect - i.e. "what do we make it to cancel $C_\psi$".
Note that we can't go the other way, i.e. redefine $C_\psi$ such that $a^{-2}$ is no longer needed - at least not entirely; per above derivation, it's clear that $\psi$ must somehow drop, thus $C_\psi$ must include it. The usage of exactly $C_\psi$ might be convention-based, since $C_\psi$'s applicability spans much outside the iCWT as the fundamentally defining criterion of wavelets (Farge, 1992).
So how to invert L1-normed CWT? Obvious from above.
Notes:

Proof followed Mallat, Theorem 4.4
Split integrals based on identity
Assumed $f \in \mathbf{L^2}(\mathbb{R})$ (square-integrable)
$\psi$ was also assumed real, and note the tilde is missing. I'm unsure why such a restriction was made, or whether Mallat implies that all real $\psi$ are self-dual, but the rescaling identities (e.g. $a \Rightarrow a^{-1}$ via $\mathcal{F}$) should hold regardless. Derivation 2 does not restrict $\psi$ to real.

Derivation 2: Resolution of Identity
For all $f, g \in \mathbf{L^2}(\mathbb{R})$ (real or complex), below holds:
$$
\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{da db}{a^2} W_f(a, b) \overline{ W_g (a, b)} = C_\psi \left< f, g \right> \tag{5}
$$
Begin proof by transforming to frequency domain:
$$
\begin{align}
(5) = \int \int \frac{da db}{a^2} & \left[
      \int d\xi {\hat f}(\xi) |a|^{1/2} e^{-jb\xi} \overline{ {\hat \psi}(a\xi)} \right] \\
      \cdot & \left[ \int d\xi' \overline{ {\hat g}(\xi')} |a|^{1/2} e^{j b\xi'} {\hat \psi}(a\xi')
      \right] \tag{6}
\end{align}
$$
This is similar to derivation 1, except $\psi$ replaced with another CWT of some different function $g$. In the brackets we have the Fourier transform of $F_a(\xi) = |a|^{1/2} {\hat f(\xi)} \overline{\hat \psi (a\xi)}$ times the Fourier transform of complex conjugate of $G_a(\xi) = |a|^{1/2} {\hat g(\xi)} \overline{\hat \psi (a\xi)}$. By the unitarity of the Fourier transform, we have
$$
\begin{align}
(6) &= 2\pi \int \frac{da}{a^2} \int d\xi F_a(\xi) \overline{G_a(\xi)} \\
    &= 2\pi \int \frac{da}{|a|} \int d\xi {\hat f(\xi)} \overline{ {\hat g}(\xi)} |{\hat \psi}(a\xi)|^2 \\
    &= 2\pi \int d\xi {\hat f(\xi)} \overline{ {\hat g}(\xi)}  \int \frac{da}{|a|} |{\hat \psi}(a\xi)|^2 \ \ \ \ \text{(Fubini's theorem)}  \\
    &= C_\psi \left< f,g \right>,
\end{align}
$$
where in last step a similar change of variables was used, $\zeta = a\xi$. $(6)$ can be read as
$$
f = C_\psi^{-1} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{da db}{a^2} W_f(a, b) \psi^{a, b}, \tag{7}
$$
quoting the source,

with convergence of the integral "in the weak sense", i.e., taking the inner product of both sides of $(7)$ with any $g \in \mathbf{L^2}(\mathbb{R})$, and commuting the inner product with the integral over $a, b$ in the right-hand side, leads to a true formula.

This seems to say that the original derivation was for the inner product $\left< f, g \right>$, whereas $(7)$ reads only for $(6)$, which will have somewhat different properties. The "original" is recovered by convolving with a $g$. -- I'm also unsure what exactly "in the weak sense" means; this is the most I found (clarifications welcome).
Once again, it is evident that $a^{-2}$ was set in retrospect, depending on how the forward transform $W$ was normalized.
Source: Daubechies, Ten Lectures on Wavelets, Chapter 2


Circular Convolution as Cyclic Shift Operator Given the following signal vectors:
$$ γ=[ψ_0,0,ψ_1,0,ψ_2,0,…,ψ_{N-1},0]^T\in \mathbb{R}^{2N}, ϕ=[1,\frac{1}{2},0,…,0,\frac{1}{2}]^T \in \mathbb{R}^{2N}$$
I want to show that the convolution of $γ$ and $ϕ$ is actually a cyclical-shift operator. i.e.:
$$γ*ϕ = γ + \frac{1}{2}T_1(γ) + \frac{1}{2}T_{2N-1}(γ) $$
Where $T_{t_0}$ is cyclical-shift operator with offset of $t_0$ places.
I tried to develop it according to the definition of (circular) convolution:
$$h_l=∑_{n=0}^{2N-1} γ_n ϕ_{l-n (mod \space 2N) } $$
However, I couldn't come up with something helpful...
How can I prove this identity? Thanks.

You may solve it by 3 steps:

Show yourself that a Cyclic Convolution with a vector $ \boldsymbol{e}_{i}^{N} $ is a Cyclic Shift Operator $ {T}_{i - 1} \left( \cdot \right) $. Where $ \boldsymbol{e}_{i}^{N} $ is defined as a vector of length $ N $ which all its elements is zero but the $ i $- th element which is 1.
Decompose $ \phi $ into $ 3 $ vectors of type $ \boldsymbol{e}_{i}^{2 N} $.
Use linearity of the Cyclic Convolution to prove your assertion.


Bluetooth locating with RSSI and only one receiver Is there a method to find the relative direction of a Bluetooth device to the receiver?

A quick Google search on "Bluetooth positioning" returns several methods using triangulation, but I don't want to use multiple receivers. Also, I don't need to know the distance, only the relative angle to the receiver is enough.

If the environment does not contain a lot of 2.4 GHz radio reflecting objects, you can try shielding one side of the BTLE receiver's antenna (it might be tiny) with a large conductive object, the rotating the shield and antenna configuration.   The difference between RSSI's will be greatest when the shield is between the source and receiver antenna compared to the opposite orientation.


Audio localization accuracy with two microphones I am attempting to find the location of sound using two microphones, using cross-correlation methods. I am not getting very consistent and accurate results, there are instances where there seem to be some kind of noise signals that interfere with the signal of interest and generate rubbish values of angles of arrival.

I have tried looking at various methods to make it more accurate, but now I'm wondering whether I should be using more than two microphones to localize the sound. Is it possible that my decision to use two microphones only is causing the poor accuracy of localization?

Each ear have a Frequency Response $H_L(f)$ and $H_R(f)$ and transmission delays $\tau_L$ and $\tau_R$ according the 3D angle from where the sound source is located relative to each ear.

Hence, for a fixed magnitude and frequency sound, $s(t)$, you have the response $H_L(f)s(t-\tau_L)$ and $H_R(f)s(t-\tau_R)$.

Actually, you cannot solve this for every position, or equivalently, you need a head movement (or ear movement in animals) for retrieving the proper angle. Besides, there exist some cone of confusion which means there is actually not a well defined 3D angle, mainly you only can locate 2D cones instead of location.

I will not solve it for you, but from that information you can trivially determine a location considering head movement onto some estimated direction, or an enough rich movement trajectory.


Cross-correlation - Why is there a complex part after FFT/IFFT when deleting the mirrored part of the spectrum? I want to detect a special shape in my time series and I apply an matched filter using cross-correlation to increase the SNR.

e = conj(fft(shape,2048));
f = fft(signal,2048);
g1 = real(ifft(e.*f));

However, when deleting the second half of the FFT results, because of periodicity, then my g2 still has a small complex part and the result is less good in terms of SNR.

_e = e(1:length(e)/2);
_f = f(1:length(f)/2);
g2 = real(ifft(_e.*_f));

Could you explain me why there remains a complex part in g2?
And is there any easy method to improve the signal again after deleting half of the spectrum? Or is not advisable at all to discard the second half?

If the input signal is real-valued, then the DFT returns a sequence of length $N$ that satisfies $$X[N-n]=X^*[n]$$

As one usually plots the magnitude of the spectrum, then the relationship of interest is:

$$|X[N-n]|=|X^*[n]|\implies|X[N-n]|=|X[n]|$$

When you work with real-valued signals, one tends to "ignore" the second half of the FFT spectrum due to this symmetry property.

However, if the signal is complex, then both sides of the spectrum are fundamental, as the symmetry is not accomplished anymore. When you call the function ifft after discarding one side of the spectrum, the function is assuming that the original FFT (the input) is non-symmetrical (as you deleted one half of it), so the time-domain sequence from which it came from must be complex. That's why you get a complex-valued signal when coming back to the time-domain.


FFT of a signal with 0 added between samples I'm having a rough time doing the following question from a HW. Given a sequence of your choice, what happens to the DFT if we add a 0 between every sample? I.e. if we have x[n] = [A,B,C], we turn it into x[n] = [A,0,B,0,C,0].

My initial thought was that the sampling frequency would double, similar to when we zero pad at the end of the signal. To do some analysis, I used the following MATLAB code.

Fs = 300;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

S = 0.7*sin(2*pi*50*t); %Signal - sine 50 Hz
z = zeros(1,3000);
z(1:2:end) = S; %Zero interp
z_fft = fft(z);
z_abs = abs(z_fft/(L*2));

P = z_abs(1:L*2/2+1);

P(2:end-1) = 2*P(2:end-1);
f = Fs*(0:(L*2/2))/(L*2);
plot(f,P);

Which yields the following plot:

That is different from the result when I add the same number of 0's but to the end of the sequence:

So I got confused, what is the result of adding the 0 between every sample? Does that doubles the frequency of my signal? What about that peak on f = 125Hz? Or am I making something wrong on my simulation?

Thanks for your time!

You can use expansion analysis to deduce what happens when a sequence is zero staffed in between its samples.

Let $x[n]$ be your original length $N$ sequence. Then its expansion by $2$ yields the new sequence $y[n]$ of length $2N$ :

$$ y[n] = \begin{cases} { x[n/2] ~~~,~~~n = 2m \\ ~~~~~ 0 ~~~,~~~ \text{otherwise} } \end{cases} $$

Then the DTFT relation between $Y(e^{j\omega})$ and $X(e^{j \omega})$ will be : $$Y(e^{j\omega}) = X(e^{j 2 \omega})$$

Then the DFT relation between N-point $X[k]$ and 2N-point $Y[k]$ will be

$$Y[k] = Y(e^{j \frac{2\pi}{2N} k}) = X(e^{ j 2 \frac{2\pi}{2N} k}) =  X[k] ~~~,~~~ k = 0,1,...,2N-1$$

In other words, $Y[k]$ will be two copies of $X[k]$


Understanding the meaning of amplitude in FFT I am recording data with a magnetometer of the background magnetic field in a building. I have applied the FFT algorithm to the data in order to look for the frequencies that appear in it. I would like to use this in order to identify (or at least make an educated guess) of the sources of the disturbances that I observe.

My question is: What is the meaning that I can attribute to the amplitude that I obtain from the FFT algorithm? Is there some unit that can be ascribed to it?

Looking at the formula for the continuous fourier transform (which I took from Wolfram Mathworld) :
\begin{align}
f(\nu)=\int\limits_{-\infty}^{+\infty}f(t)e^{-2\pi i \nu t} \mathrm{d}t
\end{align}
I do not really know how to accomodate the dimension of Tesla in there.

Thank you

The continuous-time Fourier transform of a function $f(t)$ is in essence an integration of $f(t)$ multiplied with a complex exponential kernel:
$$ F(\omega)=\int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt \tag{1}$$
Since the exponential function is unitless, the unit of the Fourier integral will be the multiplication of the units of the function $f(t)$ and the differential $dt$.
Assuming that the function $f(t)$ had a unit of micro Tesla, and its argument $t$ is time (in seconds), then the unit of $dt$ will be seconds. As a consequence, the unit of the Fourier transform, $F(\omega)$, will be  micro Tesla  second  $$\mu T \cdot s \tag{2}$$
However, what you actually compute is the discrete-time Fourier transform, $F(e^{j\omega})$, of the samples $f[n]= f(nT_s)$ of the function $f(t)$, via the summation: $$F(e^{j\omega}) = \sum_{n=-\infty}^{\infty} f[n] e^{-j\omega n} \tag{3}$$ where $T_s$ is the sampling period in seconds.
Furthermore, instead of the continuous-argument function $F(e^{j\omega})$, you will compute its samples $F[k]$
$$ F[k] = \sum_{n=0}^{N-1} f[n] e^{-j \frac{2\pi}{N} n k} \tag{4}$$
through a DFT (discrete Fourier transform) of the samples $f[n]$ of length $N$, possibly implemented with an FFT algorithm.
The unit of the samples $f[n]$ is the same as the unit of $f(t)$, making the unit of $F(e^{j\omega})$  and $F[k]$ as micro Tesla. Therefore the unit of the FFT samples $F[k]$ of $F(e^{j\omega})$ will be micro Tesla.
Note that there is an (implicit) amplitude scaling by $1/T_s$ in the computed DFT samples $F[k]$, and when you want to display the continuous-time Fourier transform $F(\omega)$ from the samples $F[k]$, you multiply them with $T_s$, which corrects not only the amplitude scaling, but also the unit of it, by making it micro Tesla second as in (2).


Difference between 2D DFT's and 1D DFT's of Linearized Matrices I have recently left the safe and easy MATLAB environment and begun to use CUDA-C/C++ for image processing.   Since CUDA doesn't allow 2D arrays to be passed into kernels I am now used to linearizing my gray-scale images (2D matrices) into 1D arrays before processing them.  I seem to have no problem applying 1D-DFT's to my linearzied arrays.  This leads me to wonder why people use 2D-DFT's at all?   Why not just linearize your matrices into 1D-arrays and apply 1D-DFT's?

Can 1D-DFT's always operate on 2D-data by simply laying the matrix out into a 1D array?

The DFT is "Separable Operator" (Also the classic Fourier Transform) and hence can be applied on the Rows and Columns of the image separately (It can be generalized to N dimension and not only 2).
Yet, if you create  1D signal from your image (Let's say by Column Stack) and apply 1D DFT you don't get the information you would by using 2D DFT (By going on the Row and them Columns).

Remember, Fourier Transform is all about synthesizing the signal using different functions. In this case if it is 2D signal you want to build it using 2D Signals.

You need to play with the indices to first apply it on each row and later on each column.


Does perceived ADC resolution, range and precision depend on fixed point format? What is the best way to deal with this for signal processing? I am working on a system that is implementing a digital filter that was designed for another system that had an FPU and a 12bit ADC. The system that I am working on is using the same filter without an FPU using a 12 ADC and fixed point arithmetic.

The ADC resolution is 3.3v / 2^12 = 0.000805664062v

With an FPU each ADC integer value could be scaled with the above resolution.

On the system using fixed point where

m2^-e, e = 10

The resolution would be

2^-10 = 0.0009765625

This leads me to think that the use of fixed point is giving me slightly lower resolution.

Looking at it differently and representing a delta between ADC integers with 14bits for the fraction the resolution could be

2^(-11) + 2^(-12) + 2^(-14) = 0.000793457031

but then the range would be perceived as less than 3.3V

Is this thinking correct?

it may turn out to be very inconvenient, but you can do essentially anything in fixed-point that you can do in floating.  as long as you got "Add-with-carry" as an instruction and enough MIPS to do all of the multiword instructions you may need.  the 12 bits that come from your ADC don't care if it's fixed or floating and your ADC resolution is the same regardless.  can't promise that the quality of code coming out of a C compiler will be any good.


Realtime filtering on a microcontroller, with dynamic cut off frequencies Is it possible to dynamically filter a signal in realtime on a microcontroller?

What I mean by dynamically is to change the cut off frequency of the filter on the fly.

I have some DSP knowledge and realize that I probably do not have the computational power to do an FFT on the fly. Instead I should be looking at creating an IIR or FIR filter that I can continually convolve with my signal.

The issue is now whether or not I have enough computational power to generate the filter coeffecients quickly and then convolve with the signal. I'm guessing that because of this I'm going to have to be very lenient on the response of the filter which is ok for my application.

I have looked at algorithms such as the Remez Exchange algorithm and believe this may be able to work but am not sure. I wanted to check before I start the lengthy process of implementing whether or not I am on the right track.

You can pre-compute the coefficients for a bunch (say $N=1024$) different cutoff frequencies. If your application doesn't need to provide fine control over the cutoff frequency, just read from your table and you're done. This works best if you have a lot of ROM and/or you use short FIR or IIR filters.
If you need finer control than those $N$ values, and if the relationship between your coefficients and cutoff frequency is well-behaved (ideally, monotonic), you can use linear interpolation within the pre-computed table of coefficients. This works best with well-behaved structures like state-variable filters - for which the coefficients are monotonic functions of the cutoff frequencies, and which do not get unstable when modulating them at fast rates. This is a very common problem in musical audio applications (equalizers, synthesizer analog filter emulation...). If you use other filter topologies, make sure that they cannot get unstable when their coefficients are changed "in mid-flight" - to me, this is a harder problem than actually computing them.


Comparing 2 sinusoids: Discrete case vs Continuous case I'm doing MIT's OCW course with Allan Oppenheim - "Signals and Systems".

I'm given 2 continuous signals$$x(t)=\cos[\omega_x(t+r_x)+\theta_x]\tag1$$ and $$y(t)=\cos[\omega_y(t+r_y)+\theta_y]\tag2$$

Assuming their phases are the same, the notes tell me that the signals are identical if $$\omega_x=\omega_y$$which is quite obvious.

Now, the discrete case:

Again, I'm given 2 signals $$x[n]=\cos[\Omega_x(n+P_x)+\theta_x]\tag3$$ and $$y[n]=\cos[\Omega_y(n+P_y)+\theta_y]\tag4$$

This time, assuming that the phases are the same, the notes say that the signals are identical if $$\Omega_x=\Omega_y+2\pi k \tag5$$

Why is the $2\pi k$ present in equation 5?

Because of the $2\pi$-periodicity of the cosine and the discrete nature of $n$ ($n,k\in \mathbb{N}$):
$$
x[n] = \cos[\Omega_\mathrm{x}] = \cos[(\Omega_\mathrm{y}+2\pi k)n] = \cos[\Omega_\mathrm{y}n+2\pi kn] = \cos[\Omega_\mathrm{y}n] = y[n]
$$
This important "phenomenon" is called aliasing and only occurs for discrete-time signals. The normalized frequency $\Omega$ is therefore usually restricted to the interval $[-\pi, \pi]$.


Terminology: Camera Tracking, Camera Relocalization I am reading the paper RGBD Datasets: Past, Present and Future available here, and two terms that I am not familiar with are presented.
On page three, middle left side of the page, the author refers to

Camera Tracking
Camera Relocalization

From what I understand, camera tracking is when you have a sequence of 2D frames, and you track the position of an object within those frames to determine the camera's pose as it moves in space. I am wondering if this is correct?
Further, I am not sure what camera relocalization means, and some Google searching did not yield a conclusive answer. Any insights appreciated.

In the context of SLAM papers, camera tracking means "local relative pose estimation of different views of the scene viewed by a given camera". It's also called visual odometry. What is important in this task is to get accurate estimates of from-to-frame (or in a small time window) camera displacements. It's to make a difference with mapping where you have global map optimization applied (think loop closures or global bundle adjustment). I refer you to the PTAM paper for more details as to why these 2 tasks are considered.
Relocalization means your robot has a map of a known environment, then you kidnap and blind it, and when you uncover its sensors you want it to recognize what is its new position with respect to the map it already knew (you did not bring it to an unknown place).


I'm looking to do my senior capstone project in the field of medical tomography, what is a good book to get acquainted with the field? I'm a junior currently pursuing a mathematics degree with aspirations toward medical school. I have a senior project I have my entire last year to devote time to, and I've decided that I want to tackle the challenge of medical signal processing. I'd like to write a program that reconstructs a 3d-model of the object from the cross-sectional images provided to it, very much like a CT scan. However, I have no idea where to start. You can assume that I have a good level of mathematical knowledge or the drive to acquire what is needed, so don't be afraid to throw some academic papers my way! I'd really like to dive in to my research as soon as possible, but I can't quite figure out where to start. What materials would you recommend to someone starting out, such as myself?

What about for a start:

Buzug, Computed Tomography
Hermann, Fundamentals of Computerized Tomography


Matlab applying an FIR filter to a signal, generated in Matlab I want to create an FIR in Matlab and apply it to a signal (also generated in Matlab). You will find the code below. I will put my questions in order from here.

1) The filter coefficients vector "b" should contain coefficients for a bandpass with a filter order of 100 and a passband of 1000Hz to 2000Hz. My sampling frequency is 50.000Hz. Matlab uses normalized frequency. As far as I have understood, my formula to transform my frequency to a normalized frequency is f_norm = (2*f)/f_sample. So in my example I would get 0.04 for 1000Hz and 0.08 for 2000Hz. Is that correct?

2)I am creating a test signal sinuses with the frequency 200Hz, 2500Hz and 5500Hz and compute the two sided fft of the amplitude. The graph is attached

This looks good, since the overall amplitude should be about 0.33 for every signal part. Now, when I apply my filter to the signal and plot the two sided fft of the result, it shows the next graphFFT Signal after filter

I mean I can see the 5500Hz being filtered out, but the amplitude of the 2500Hz signal drops by a factor of 10 and the 200Hz signal is also not filtered out completely. What are some things that I can improve here? 3) Also, when I only raise the filter order, I am getting worse results with more suppression on the 2500Hz signal and less suppression on the signals that I want to have filtered out.

%fs = 50000
b = fir1(100, [0.04 0.08], 'bandpass');    %passband 1000Hz bis 2000Hz
figure(1);
freqz(b,1,2001);

f1 = 100;                  %Hz
f2 = 5500;
f3 = 2500;
A = 1.00;                  %Amplitude
sampleRate = 50000;
sampleTime = (1/sampleRate);
t = 0:sampleTime:1;

s1 = A/3 * sin(2 * pi * f1 * t);
s2 = A/3 * sin(2 * pi * f2 * t);
s3 = A/3 * sin(2 * pi * f3 * t);

signal = s1 + s2 + s3;

Y = abs(fft(signal));
samplePoints = length(Y);
Y = Y / samplePoints;
F = linspace(0,50000,50001);

figure(2);
plot(F,Y);
grid on;

y = filter(b,1,signal);

Y = abs(fft(y));
samplePoints = length(Y);
Y = Y / samplePoints;
F = linspace(0,50000,50001);

figure(3);
plot(F,Y);
grid on;

In general your result match expectations, I don't see anything in there that's unusual or unexpected.

Filter design is a complicated trade off between complexity/cost, stopband attenuation, pass band ripple, transition steepness, phase distortion, time domain ringing, non-causality, latency, etc. The best way to go about is to clearly articulate the requirements for your specific application and then evaluate any design against these requirements line-by-line.

I'd also be careful with spectral analysis: direct FFT works sometimes, but something like pwelch() is typically much more robust.


Why the calculated value of power is not match? I am trying to calculate the power of a signal by using following two methods one is a power of a signal and another one is using the periodogram in Matlab as shown in the code.But why I am getting a different answer for the power?

x=[1 2 3 4 ];

s=0;

for i=1:length(x)
    s=s+x(i).^2;
end

p=s/length(x);

D=fft(x);

P= sum((abs(D).^2)/length(x));

p=7.5,
P=30

These are the formulas for both methods.

Discrete signals are by definition periodic and you need to calculate the power in both domains over one period with proper normalization.

In practice this means, you need to divide the spectrum by the FFT length. See https://en.wikipedia.org/wiki/Parseval%27s_theorem


What will be the code for this following signal in MATLAB? I was trying to generate a signal like this:

I have done this so far like this:

clc;
clear all;
close all;
t = 0:.1:10;
f=.5;
y=.2*sin(2*pi*f*t);
plot(t,abs(y));
axis([0 10 -2 2]);
grid on;

What can I do to get desired output?

Well, so far using your code, you are getting the following:

clc;
clear all;
close all;
t = 0:.1:10;
f=.5;
y=.2*sin(2*pi*f*t);
plot(t,abs(y));
axis([0 10 -2 2]);
grid on;

So now, the only missing part is to zero-out all samples using PWM signal with 25% duty cycle.
For example, this can be done using the square wave. Its values are between +1 and -1, so simply add 1 and use the signum function (or divide it by 2):

Putting all that together:

clc;
clear all;
close all;
t = 0:.1:10;
f=.5;
y=.2*sin(2*pi*f*t);
s=sign(square(t/2*pi, 25)+1);
x=y.*s;
plot(t,x);
axis([0 10 -2 2]);
grid on;


Finding the Harmonics contained in an Image? I am beginer in image processing. Any signal whether it is 1D,2D or any multidimensional signal can be represented using combination of number of sine and cosine wavesforms (harmonics). Similerly any image can be termed as a function of sinusoidal signals.

I want to see individual pattern for the number of waveforms/harmonics present in an image.for example the following  link contains few of them http://www.cs.toronto.edu/~rfm/factored/filters_out.png . So is there any code in Matlab/ OpenCV  or results to understand easily these components(sine and cosine functions(harmonics)) present in any image?

this questions is asked fairly often. If you don't have a grasp of fft in 1d, higher dimensions may be difficult to grasph. But logically it makes a lot of sense once you see it What does frequency domain denote in case of images?

basically sharp changes (high contast noise) usually has high frequencies. and smoother gradients have low frequencies. just as a quick crash course

edit

sorry i thought the other post had code, its just code agnostic conceptual stuff. anyway ill tell you matlab (since your tag says matlab in it) and it is more familiar to me than openCV

function fft_im = createFFTIm(grayscale_im)
    fft_im = fft2(grayscale_im);

    figure()
    subplot(2,2,1:2);imshow(grayscale_im); title('original');
    %fftshift puts the 0 frequency in the middle of the plotting area, normally
    %it would be split, part of it in the 4 corners of the image
    %we add 1 because log(0) is -inf, which we dont want
    %we scale it jsut to show the colors better
    subplot(2,2,3);imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
    title('magnitude spectrum');

    subplot(2,2,4);imagesc(angle(fft_im));  colormap(gray);
    title('phase spectrum');
end

edit #2 i slightly modified the code

and here are two examples

notice with the chess example, there are visible peaks in the magnitude spectrum. Just like the regular fft, these peaks correspond to a sing wave of a specific frequency. Since out image is periodic there is are strong sines at the corresponding frequencies. You also notice since the image is near symmetric, the phase is also periodic.
The phase is responsible for modeling the fact that an image is not symmetric

edit #3
this decomposes the figure into individual elements. In the code I limited it to a 4x4 subblock of pixels (the upper left portion of the figure) but you can change this all the way up to the size of the image (though for any large image this would be incredibly slow)

As in my comment, I processed the image one pixel at a time, and did an FFT2 on each of those 1 pixel images. The images themselves are still the same size as the original, it just has only one pixel value populated

function fft_im = decomposedcreateFFTIm(grayscale_im)
    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper left corner of the image
    max_rows_to_process = 4;
    max_cols_to_process = 4;

    [m,n] = size(grayscale_im);

    %this is used to isolate pixels, for individual analysis
    zer_array = zeros(m,n);

    figure(1);title('magnitude')
    figure(2);title('phase')

    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper right corner of the image
    m=4;n=4;

    %for all rows
    for ii=1:1:max_rows_to_process
        %for all columns
        for jj=1:1:max_cols_to_process
            curr_idx = (ii-1)*max_cols_to_process + jj;

            %creates an image with only one pixel
            indiv_pix_pic = zer_array;
            indiv_pix_pic(ii,jj) = grayscale_im(ii,jj);

            %does fft
            fft_im = fft2(indiv_pix_pic);

            %fftshift puts the 0 frequency in the middle of the plotting area, normally
            %it would be split, part of it in the 4 corners of the image
            %we add 1 because log(0) is -inf, which we dont want
            %we scale it jsut to show the colors better
            figure(2)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(angle(fft_im));  colormap(gray);

            figure(1)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
        end
    end
end

magnitude

phase


How can the Poles of the Root Locus be negative? My understanding of drawing a root locus diagram is that stability requires all roots of the characteristic polynomial of the open loop transfer function to lie in the negative real part of the plane. In other words, it must have roots $(s+a)(s+b)...(s+n)$ where $s$ is a negative number.

My question is that the Laplace Transform requires that the real part of $s$ be positive, so wouldn't the possibility of negative real parts of $s$ be a contradiction?

Thanks.

If the denominator polynomial is given by

$$Q(s)=(s-s_0)(s-s_1)\ldots (s-s_{n-1})\tag{1}$$

then the corresponding time domain function is made up of a sum of scaled exponentials of the form $e^{s_it}u(t)$, where $u(t)$ is the unit step function. These terms decay only if $\textrm{Re}\{{s_i}\}<0$, so for stability the zeros of the denominator polynomial (i.e., the poles of the system function) must lie in the left half-plane.

Let $h(t)$ be the impulse response of a causal LTI system. The transfer function is the Laplace transform of $h(t)$:

$$H(s)=\int_{0}^{\infty}h(t)e^{-st}dt\tag{2}$$

This integral only converges in a region $\textrm{Re}\{s\}>c$ ("region of convergence"), with some real-valued constant $c$. That constant depends on the properties of $h(t)$. If $h(t)$ is a weighted sum of exponentials $e^{s_it}u(t)$ then $c$ equals the most positive real part of the complex numbers $s_i$:

$$c=\max_i\big\{\textrm{Re}\{s_i\}\big\}\tag{3}$$

This makes sure that the damping by the term $e^{-st}$ in $(2)$ is enough to make $h(t)e^{-st}$ decay sufficiently such that the integral converges.

For stability we require $c<0$, i.e., all poles lie in the left half-plane, and the imaginary axis of the $s$-plane $s=j\omega$ lies inside the region of convergence.


Discontinuity in the angle of a complex exponential signal Plotting this complex expontential:

$$z[n]=100e^{j0.1n},\quad\text{we have a discontinuity at}\quad \arg[n]=10\pi$$

I guess this is related to the angle formula:

$$\theta=\arctan\left(\frac{\Im\{z[n]\}}{\Re\{z[n]\}}\right)$$

But I do not quite see the relationship. Can anyone tell me exactly why we have a discontinuity at $10\pi$?

That's not really a discontinuity. On a circle the two points $-\pi$ and $+\pi$ are identified: They are the same point. That is true for all $x$ and $x+n 2\pi$ for integer $n$.

If you would like to keep the continuity obvious, you might want to switch to a phasor representation of angles. Those are complex numbers on the unit circle $\exp(i\phi)$ and they form a group under multiplication. They nicely absorb the non-obvious topology of the phase angle.


How do I apply a binary mask and STFT to produce an audio file? So here's the idea: you can generate a spectrogram from an audio file using shorttime Fourier transform (stft). Then some people have generated something called a "binary mask" to generate different audio (ie. with background noise removed etc.) from the inverse stft.

Here's what I understand:

stft is a simple equation that is applied to the audio file, which generates the information that can easily be displayed a spectrogram. By taking the inverse of the stft matrix, and multiplying it by a matrix of the same size (the binary matrix) you can create a new matrix with information to generate an audio file with the masked sound.

Once I do the matrix multiplication, how is the new audio file created?

It's not much but here's what I've got in terms of code:

from librosa import load
from librosa.core import stft, istft
y, sample_rate = load('1.wav')
spectrum = stft(y)
back_y = istft(spectrum)

Thank you, and here are some slides that got me this far. I'd appreciate it if you could give me an example/demo in python

The idea of the short-time Fourier transform (STFT) in this case is to compute a representation of the input signal where we see how the frequency content of the signal evolves over time. If you aren't familiar with the STFT, I would suggest you read up on it a bit first. The Wikipedia article on the subject should get you started.

To your question. The STFT function in your code gives you a complex-valued matrix in the output. When you move up or down a column in the matrix you are moving in the frequency direction. When you move along a row in the matrix you are moving along the time direction. The magnitude of a single complex-valued element in the matrix gives you an estimate of the energy in a certain frequency band at a certain time position in the input signal.

I recommend you try to understand the parameters that are given to you in the function call. You should try to understand how the window size, hop size, number of FFT bins, and the windowing function used affect the output of the function.

If we denote your complex valued STFT with $$X = \textrm{STFT\{x(n)\}},$$ and the binary mask with $M$, the output signal in the time-frequency domain is given by:

$$Y = X \odot M,$$
where $\odot$ refers to element-wise multiplication. Then you want to convert your modified time-frequency representation $Y$ back into to a time signal. This is done by doing the inverse STFT:

$$y(n) = \textrm{ISTFT}\{Y\}.$$

This should get you going, but time-frequency processing is a complex topic and I highly recommend trying to find resources to help you really understand what is going on when you are doing things such as applying a binary mask in the STFT-domain.


Spectrum of FM signal? If a carrier ($f_c$) is FM modulated by a single sine wave ($f_m$), its spectrum is composed of frequencies at $f_c + kf_m$ ($k \in \mathbb{Z}$) whose amplitudes are weighted by $J_n(\beta)$, the $k$-order Bessel function of modulation index $\beta$. We could get such spectra :

Now, what would be the spectrum of the carrier if it was modulated by a more complex signal like speech or music ?
Short term Fourier transform of the modulating signal would be :

Harmonic (e.g. : vowel in speech, or a violin note)
Non harmonic (e.g. : some consonants in speech, some percussive instrument)

In both case, modulating signal will have a more richer frequency content than a single sinusoid.
Will these rich spectra just repeat regularly around $f_c$ as in the case of a single tone modulating signal ? Why ? What would happen in the simple example of a 2 tones modulating signal ?
If a carrier is AM modulated by complex signal like speech, we find the spectrum of modulating signal shifted around the carrier (as seen in the spectrogram below) : what would happen with FM modulation for such rich signals ?

Because frequency modulation is a nonlinear phenomenon, there's no good general-purpose answer to your question.  You can take any specific baseband signal and figure out its modulated spectrum, but that's it.
For lack of a better technical term, FM modulation takes the frequencies in the parent signal and munges them together, spreading them out by more than the deviation frequency.  For a sine wave this is that Bessel function.  For a pair of sine waves, it'd be all the harmonics of each sine wave plus all possible combinations of intermodulation ($n f_1 + m f_2\, \forall\, n, m \in \mathbb{Z}$) -- so for the two-sine case there may be some discernible, identifiable lines in the output.  For sound with a continuous spectrum, you'll get a continuous spectrum out, but nothing you could just recognize on a spectrum analyzer.
For practical purposes, you're left with rules of thumb, inspired by the above results but limited by our inability to solve the nonlinear differential equations in any general way.  If you have a signal with bandwidth $B$, and a maximum frequency deviation of $f_D$, then the bandwidth of the transmitted signal can be limited to $2(B + f_D)$.  This isn't satisfying from a theoretical standpoint, but it's backed by over 70 years of practical experience with broadcast and voice communications systems.


FFT of High Frequency and Low Frequency? When we plot FFT in a matlab it doesn't ask for Frequency , all we need to input the sequence and it simply shows the output using stem command.
If i want to compute two different signal using FFT, one with High and the other with Low frequency what should i suppose to do ?
and what will be the difference between the two of them ?

Here is simple code for the FFT calculation ,

x=[1 2 1 0]
n=0:3
y=fft(x)
subplot(2,2,1)
stem(n,y)

The frequency span of an FFT is determined by the number of samples fed into it and the sample rate of those samples.  Basically, the maximum frequency analyzed by the FFT will be half the sample rate.  The lowest frequency will be DC, and the remaining output "buckets" will be evenly spaced in frequency up to the (last) max frequency bucket.  The number of output "buckets", in the standard implementation, is half the number of input samples.

(Actually, the total number of outputs is equal to the number of input samples, but half the outputs are the imaginary half of complex numbers ... and at about this point my head starts to hurt.)


How should calibration points be distributed? I do camera calibration for a single camera with a chessboard pattern in opencv for the first time.
The focus is fixed.
I'm not sure how to do it, there seems to be many questions:

How many images should be used for calibration?
What pattern should be used? number of points, size of squares
What transformations should be applied to the calibration target? rotated, tilted, shifted...

I came to the conclusion that none of the questions above can be answered on their own, because they all depend on each other. Say for example taking a lot of images that are very similar to each other are intuitively not as good as having less images that are better distributed in the frame.

I took a step back and thought about what's going on. (I tried)
Basically, there's a mathematical model of the camera and its aberrations (focal length, distortion, etc...) and the calibration is the process of collecting data to fit the model (its values) to the hardware at hand. This appears to be mostly statistics. In order to fit the model to the entire image well, sample data should be collected from the entire image.

I conclude that all questions above boil down to how should the detectable points of the pattern be distributed in the 3D space.

I'd like to know if this thought makes any sense.

I'm tired of reading handwavy suggestions like "the more images, the better", "at least arbitrary number of images should be used", etc. that do not provide a quantifiable measure for the quality of a calibration. Sure, there's a reprojection error, that tells how good the calibration fits to the provided calibration data, but not how good that data is suited for calibration in the first place.

I came to the conclusion that it makes more sense to evaluate the quality of a calibration by the distribution of the sample data used to create it. Makes sense?

Questions 1 and 2 above are thus combined to : What should be the overall density (or density distribution) of detected calibration points among all calibration images?

Say I take all detected points from all calibration images and place them in one image, I could see how many calibration points there are per pixel area $\text{px}^{2}$.If I have a lot of images with the pattern in the center of the image, I get a higher density there.

What would be a good minimum (baseline) value for the density of calibration points per pixel area?
$0.1 \ \text{px}^{-2}$ ? $0.5 \ \text{px}^{-2}$ ?

Should the distribution of calibration points per pixel area be constant over the entire image? Or do I need a higher density (read: more sample values) at the edges for example?

Question 3 is asking for the distribution in the 3rd dimension. The calibration target with the pattern can be moved back and forth (along the optical axis) or tilted. Both transformations lead to calibration points in front or behind the plane that's in focus1. Depending on the depth of field, corner points can be detected up to some distance from the focal plane. This adds a 3rd dimension to the space in which calibration points can be detected.
How should calibration points be distributed with respect to distance from the plane of focus1? Should most points be around the plane of focus1 given that detection of points further away from it is not as good?

tl, dr;

In front of a camera, there's a 3D space enclosed by a frustum of a pyramid/cone in which calibration points can be detected. How should calibration points be distributed in that space in order to get a good calibration?

1 With that I mean the plane of points that will be "in focus" in the image.

For that I have some guidelines, which I'd like to share. Even though some of these might not be valid for certain calibration software, for typical cases, they hold:

1) The entire pattern should be visible in all sequences.
2) At least 1/4th of the image should be occupied by the pattern.
3) Preferably, the lighting should be uniform everywhere on the pattern.
(As little shadows as possible)
4) The light should reflect from the pattern to the camera (as much as possible).
Thus, objects made up of metal are not suggested.
5) There should be minimal blur. If focus/depth-of-field is an issue,
use a smaller plate or less tilts.
6) If you can't fix the camera, it shouldn't move at all.
For extrinsic calibration, moving the camera after calibrating will ruin the entire
calibration. Be careful about this as even small shakes or movements distort the accuracy.
7) Don't pre-process the calibration images.
Even if you do very smart things, you might perturb the accuracy.
8) Ideally, have 15-50 images fitting the above criterion.
This should give you 450-1500 points to optimize in total.
9) If you are printing the pattern on a piece of paper, then make sure that
contrast between black-white squares is significant and no strange shape is appears on the paper.
10) Don't bend the calibration pattern. Make sure that it is as rigid (and planar) as possible
11) If you are printing on an A4 paper, make sure you have some white space around the pattern,
remember you will want to hold it somehow.
12) If you have distortion, use small circles, they are less affected from distortion.
If you have low resolution but less distortion use larger circles.
13) Start with low tilts and fronto-parallel images.
Move to the tilts and rotations in all axes, gradually, covering the entire visible space.
If you have bias in the space coverage or distribution of points,
then it means that some regions of the space will be better calibrated.
Keep this as equal as possible. Some variation can be tolerated.

I hope these help.


Converting a lowpass filter to a highpass filter. FIR filter-type 1 I was told that in the FIR filter type I, the impulse response
of a highpass filter $h_H[n] $can be obtained as follows,
$$h_H[n] = \delta[n − m] − h_L[n]$$
where $h_L[n]$ is the impulse response of a lowpass filter, $\delta$ is a Dirac pulse (Dirac delta function), which delayed by an appropriate value of $m$.
I've computed the impulse response according to the formula in MATLAB and observed the frequency response. I noticed that when $m=\frac{M+1}{2}$, where $M$ is the number of filter taps, the formula is indeed correct. However, I don't understand Why and How it possible.

If you have a real-valued amplitude function $A_{LP}(\omega)$ of a low pass filter with unity gain in the passband, then it is easy to see that
$$A_{HP}(\omega)=1-A_{LP}(\omega)\tag{1}$$
is the (real-valued) amplitude function of a high pass filter.
A type I length $N$ linear phase low pass filter's frequency response is given by
$$H_{LP}(\omega)=A_{LP}(\omega)e^{-j\omega(N-1)/2}\tag{2}$$
Consequently, according to $(1)$, the desired high pass amplitude function is obtained by computing
$$\begin{align}H_{HP}(\omega)&=\big(1-A_{LP}(\omega)\big) e^{-j\omega(N-1)/2}\\&=e^{-j\omega(N-1)/2}-H_{LP}(\omega)\tag{3}\end{align}$$
The inverse discrete-time Fourier transform of $(3)$ is
$$h_{HP}[n]=\delta\left[n-\frac{N-1}{2}\right]-h_{LP}[n]\tag{4}$$
Note that for type I linear phase filters the filter length $N$ is odd, so $(N-1)/2$ is an integer shift. Eq. $(4)$ shows that apart from sign inversion we only modify the center bin of the low pass filter in order to obtain a type I linear phase high pass filter.


This is a major bug in this library. see http://github.com/jaidevd/pyhht/issues/43#issuecomment-398388106 I think it's better to use more robust libraries that are well-tested.


Generation of noise floor file for estimating signal-to-noise ratio of a wave file I am trying to estimate the signal-to-noise ratio of a wave file (The wave file is a recording of a conversation in a room). I found examples where the Matlab scripts compare the wave file data to a noisy wave file. This noisy wave file is generated by either mixing Gaussian noise into the actual recording or they have mixed other sounds into the wave file and estimate the SNR.

How do I determine the noise floor for calculating SNR of a wave file? The reason for estimating the SNR of the recorded wave files is to determine if the speech content recorded in the file is useful or rendered by useless due to noise in the background.

One method that is reasonable for white noise is to use a spectrogram. For each spectrum, compute a quantile value of the power bins and use that as the estimated noise floor. The quantile should be less than the maximum spectral occupancy of the speech. For example if the speech never fills more than 50% of the spectrum then the median (quantile 0.5) might be reasonable. If the speech fills more than half the spectrum then use a smaller quantile.


Find the area of all triangles from an output of Delaunay Triangulation I am detecting roadline pattern in an image. I have applied delaunay triangulation.
I get 160 triangles in that. Now I need to find the area for each triangle. Can you kindly help me. Only after this I should proceed with my project. Please help me...

Why not to use following matrix form, while given vertices pairs $(x_i, y_i) $:

$A = \dfrac{1}{6} \det \begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1\\x_3 & y_3 & 1 \end{bmatrix}  = \dfrac{1}{2}( x_3y_1 + x_1y_2+x_2y_3-x_3y_2-x_1y_3-x_2y_1)$


Energy definition for Autocorrelation lag 0 and lag 1 for complex signals I am studying the role of an auto-correlation matrix for random signals and the difference of energy between a lag 0 and lag 1 matrix.

Consider a complex input signal $x(k)=[x1,x2]^T$ and $x(k-1)=[x0,x1]^T$, as column vectors with auto-correlation matrix of lag 0, $R(0)$, with diagonal $r(0)$ and auto-correlation matrix of lag 1 $R(1)$, with diagonal $r(1)$.

where $R(0)=E[x(k)x^H(k)]$ and $R(1)=E[x(k)x^H(k-1)]$

what is the relationship between the energy resulting from $R(0)$ and $R(1)$? can we thus say that eigenvalue max of $R(0)$ is greater then that of $R(1)$
and the difference in diagonal entries $r(0)$ and $r(1)$?

given that $r(0)$ is always positive real and $r(1)$ can be complex.

What can we comment of the real and imaginary part of the $r(1)$ entry?
and where does the energy resides i.e. real portion or imaginary?
Additionally to a wide sense stationary signal $x(k)$ can we assume that $r(1)$ is approximately equal to $r(0)$?

Thus if we have $B = R(0) + R(1)$, we know that $R(0)$ is positive definite, symmetric and hermitian however $R(1)$ is not. Can we say that eigenvalue max of $B$ is eigenvalue max of $R(0)$ + the norm of that of $R(1)$ assuming its complex?

My apologies for the successive questions however I am trying to understand the differences between the auto-correlation matrix at lag 0 and lag 1 in terms of eigenvalues and energy because I have to work with a matrix $B = R(0) + R(1)$.

I had some problems understanding your questions, but from what you wrote I'm assuming you meant $\mathbf x(k) = [x(k), x(k+1)]^T$ so that  $\mathbf x(k-1) = [x(k-1), x(k)]^T$ (using bold faced letters to distinguish scalars from vectors). Is that right?

In this case, expanding your equations, we get the following: $$ R(0) = \begin{bmatrix} r(k,k) & r(k,k+1) \\ r(k+1,k) & r(k+1,k+1) \end{bmatrix}, \quad R(1) = \begin{bmatrix} r(k,k-1) & r(k,k) \\ r(k+1,k-1) & r(k+1,k) \end{bmatrix}.$$

What we can see is that $R(0)$ is Hermitian symmetric (since $r(k,k+1) = r(k+1,k)^*$) and positive definite. Otherwise, unless further assumptions are made, we cannot argue much.

If you assume that your sequence is white sense stationary, we have $r(k,\ell) = r_{\ell-k}$ and we obtain $$ R(0) = \begin{bmatrix} r_0 & r_1 \\ r_1^* & r_0 \end{bmatrix}, \quad R(1) = \begin{bmatrix} r_1^* & r_0 \\ r_2^* & r_1^* \end{bmatrix}.$$ Moreover $r_0$ represents the signal's variance ("energy") and you have $|r_1|\leq r_0$ and $|r_2|\leq r_0$ from Schwartz' inequality. In particular, you can say that $r_1 = r_0 \rho_1$ where $\rho_1$ is the normalized (Pearson) correlation coefficient between adjacent lags of your sequence. Its magnitude is in $[0,1]$ (0 for an uncorrelated sequence, 1 for a constant [or alternating] sequence) and its phase can be anything. This answers why $r_1$ is complex and how to interpret real and imaginary part. It's only close to $r_0$ if the temporal correlation is high. For uncorrelated sequences its zero and hence nowhere close to $r_0$.

Now, $R(0)$ will have eigenvalues $r_0 \pm | r_1|$, so the max eigenvalue is between $r_0$ (for an uncorrelated sequence) and $2r_0$ (for a constant sequence). The eigenvalues of $R(1)$ will be $\left(r_1 \pm \sqrt{r_0 r_2}\right)^*$. From this, one can see that the dominant eigenvalue of $R(0)$ will be greater than or equal to the dominant eigenvalue of $R(1)$ with equality only for a constant sequence.

Finally, the eigenvalues of $R(0)+R(1)$ can be shown to be $ r_0 + r_1^* \pm \sqrt{(r_1^* + r_2^*) (r_1 + r_0)}$, so it's not just the dominant eigenvalue of $R(0)$ plus something.


MRI K-space to image: why track frequencies in two dimensions? While trying to better understand how an MRI goes from k-space to an image, I came across this wonderful website that explains how you would represent an image as a collection of rows of pixels  where each row of pixels is represented as a sum of waves of greyscale intensity.
First you select a single row from an image.

Then you apply the Fourier transform to represent the row as a collection of waves.

What I don't understand is, if we can get the frequencies for each row and accurately reconstruct what a row of the image looks like via frequency encoding, and if we already know where each row is in relation to the other rows (position of the blue line), why do we need to do the same basic thing along a different dimension (the y-axis using phase encoding)? Isn't it sufficient to simply stack the rows to get the full image?

in MRI, you don't start with the image as pixels in spatial domain, but with the k-space data!
So, you don't "construct k-space from the image", you "construct the image from k-space", which takes a 2D inverse Fourier transform.
And that's not "do it row-by-row"; that is "do it on all rows, then all columns of the result".


Bringing the signal into the same level I have a signal, shown in the attachment (blue one). As you can see that the latter part of my signal has a higher dc offset. I want to bring the signal into the same level. I have tried eliminate the low frequency signal by using FFT (I FFT the signal, then set the low frequency element to be 0, then returned it by IFFT), but the result that I get is shown in red color.
I am open to any suggestion and any python code is very appreciated.

Additional info:
sampling frequency = 100 Hz
Low cutoff frequency = 0.05 Hz

First part for an answer. There are several unknowns in you question. The fist is of course:

why do you really want to bring the signal to the same level?

and the second:

what is this "same level"?

Assuming that the baseline is the lower  levels, do those simulated signals approximately match your expectations?


Multiplication and Convolution in Wireless Communications In Wireless Communication, sometimes the received signal is expressed (in time-domain) as follows:
$$r=hx+n$$
where $h$ is called a fading coefficient
and sometimes it is expressed as:
$$r=h*x+n$$
where * is the convolution operator and $h$ is called the channel impulse response.

My question is: Are the fading coefficient in the first equation and the impulse response in the second equation the same thing? If so, why we use the multiplication instead of convolution?

The short answer is that in the first case the channel delay spread spans one symbol time duration or less (which is called frequency flat fading channel), and can be approximated by a complex number according to the central limit theorem, while in the second, the channel delay spread spans several symbol duration (which is called ISI or frequency selective channel).

EDIT1: Mathematically, the flat fading channel can be represented as $$h(t)=h_k,\,\,\,\text{ for } kT_s\leq t<(k+1)T_s$$ while in the frequency selective case it can be represented as $$h(t)=\sum_{n=0}^{N-1}h_n\,\delta(t-nT_s)$$ where $N$ is the number of symbols the delay spread of the channel spans, and $T_s$ is the symbol time duration.

Now suppose the transmitted signal is $s(t)$. The noise-free received signal for the first channel is $$r(t)=h(t)\star s(t)$$ which when sampled at time $kT_s$ for the $k$th symbol equals to $$r_k=h_k\,s_k$$ while in the second case it's $$r(t)=\sum_{n=0}^{N-1}h_n\,s(t-nT_s)$$ which when sampled at $t=kT_s$ equals to $$r_k=\sum_{n=0}^{N-1}h_n s_{k-n}$$ where for brevity I used the notation $s(kT_s)=s_k$.


FIR Filter Design - input parameters conversion I am new to DSP and filter design. I have developed a code in C++ to calculate FIR coefficients using Parks-McClellan algorithm. The inputs to calculations are:

Filter type (Low-Pass, High-Pass)
Passband Freq [Hz]
Stopband Freq [Hz]
Passband Ripple [dB]
Stopband Attenuation [dB].

I got a new requirement for inputs, namely to use:
1. Cutoff freq (6dB@1.6kHz)
2. Slope [24dB/octave].

Is there a way to convert these inputs into the inputs I used in my C++ code, as listed above?

Or, more general question, how to design FIR filter with custom slope?

Thank you.

As I've mentioned in a comment, the Parks McClellan algorithm is usually used to design frequency selective filters with a fixed maximum stopband error, which results in an equiripple behaviour in the stopband. Note that the algorithm can in principle approximate any desired frequency response shape. However, many implementations just allow for piecewise constant desired responses (which includes all standard frequency selective filters, such as low pass, high pass, etc.).

If you want a frequency response that decays monotonically from the cut-off frequency, you probably want something like a Butterworth filter. This is an IIR filter which is usually designed via the bilinear transform from an analog prototype filter. If you want an FIR filter, then one option would be to truncate (or window) the infinite impulse response of the Butterworth filter.

In Matlab/Octave this can be easily done as follows. First design a digital Butterworth filter with the given specifications, compute its impulse response, and truncate it. Where to truncate depends on the error that you allow. For your example a third order Butterworth filter is sufficient. In the following example I assume a sampling frequency of $8\,\text{kHz}$:

[b,a] = butter(3,1.6/4);    % Butterworth IIR filter
h = impz(b,a,30);           % impulse response (30 coeffs)
f = logspace(0,3.6,400);    % log. frequency grid [1,4000]Hz
H = freqz(b,a,f*pi/4000);   % Butterworth frequency response
Hf = freqz(h,1,f*pi/4000);  % FIR frequency response
semilogx(f,20*log10(abs(H)),f,20*log10(abs(Hf)))
axis([1,4000,-80,5]), grid on

You see just a single curve because both filters (IIR and FIR) have virtually the same magnitude response.


Why different smartphones with same Megapixel camera have different quality of image? I am image processing student. Nowdays there are different brands of smartphones available. But the common thing is two different big brand smartphones with same Megapixel camera have different quality of image.

Moreover,the Digital cameras with same megapixel resolution have altogether best quality. Why is it so?

I am not sure this is the right place to ask (Perhaps photo.stackexchange.com?) but this is my view on the subject. The term 'Megapixel' has little if any to do with image quality nowadays, since all the popular devices have sufficient amount of pixels. Way more important is the image processing algorithms, optics and sensor size. All digital cameras have huge advantage on the last two while smartphones usually 'make up' for this with image processing algorithms, but not enough to achieve similar quality.


bode plot, zero at 0Hz I'm learning to draw magnitude bode plots from transfer functions:

$$ H(s) = \frac{sRC}{sRC+1} $$

So I can see that I have a zero at s = 0, which corresponds to 0 Hz. I know this will add a +20dB/dec slope for all frequencies in a logarithmic plot. I'm struggling to figure where the gain will start at the left side of my plot. I can kind of cheat a little bit because I know that this is a passive high pass filter, so I can start at the right side of the plot and call that 0dB, move to the left, and then apply my 20dB/dec slope once I pass the break frequency of the pole. But if the transfer function was more complicated and I didn't already know what to expect, what would I do?

$$H(s) = \frac{R C s}{R C s + 1}$$

If $R C s$ is much lower than 1 (i.e., $R C s \ll 1$, in Math), then you can make an approximation: $$H(s) \simeq \frac{R C s}{1}$$

Basically, the effect of the $R C s$ in the denominator becomes insignificant.

So, just pick a radian frequency well below $1 / R C$ and compute the magnitude of $H(s)$ at that point -- mark that on your plot, and draw your line rising at 20dB/decade from there.


why GMSK use a gaussian filter instead of raised cosine filter? I'm wondering why GMSK use a gaussian filter instead of raised cosine filter as a pulse shaping filter? Is there special reason?

(Assuming the context is FSK signals, where the pulse shape is applied to the phase, instead of the amplitude as in linear modulation).
The main reason is that the sideband and out-of-band emissions of GMSK are lower than those produced with a raised cosine filter. See for example this plot from Wikipedia: https://en.wikipedia.org/wiki/File:GMSK_PSD.png
This property is especially important for mobile communications. Since spectrum is so limited, it has to be used very efficiently. Users may use channels that are adjacent to each other, with little or no guard band between them. Any power emitted outside the band assigned to a channel resuls in degraded SINR for other users in the cell.
GMSK requires a more complex receiver, because it reduces bandwdith at the cost of introducing ISI. However, in many scenarios this is an acceptable tradeoff, since bandwidth is more valuable than receiver simplicity.


Object Tracking with Improved Detector of Objects Similar to Target I read this paper and need  your help to understand it.

Object Tracking with Improved Detector of Objects Similar to Target

In paper,  particle at the frame $t$ is represented by:
$$c^{(i)}_t(i=1,\ldots, N)$$
and after that say let $c^{(i)}_t$ used by the proposed method be defined as

$$(x^{(i)}_t, y^{(i)}_t, u^{(i)}_t, v^{(i)}_t, w^{(i)}_t, h^{(i)}_t)$$

where $(x^{(i)}_t, y^{(i)}_t)$ represents the center position of the target object, $(u^{(i)}_t, v^{(i)}_t)$ represents the velocity of it and $(w^{(i)}_t, h^{(i)}_t)$ represent the width and the height of the rectangle which involves the target object.

I can't understand if $c^{(i)}_t$ is particle and a particle is a pixel. How a particle in Particle Filter can have width and height?

How to set $x$,$y$ of the target for a particle?

I didn't read the paper but let me provide some intuition about object detection and tracking. When you try to track a target in a video, object detection algorithms might not be enough and you need to support the algorithm with target tracking methods in bayesian framework. It could be Kalman filter or particle filter.
In particle filter or any Bayesian filter you first define your state which is composed of positions, velocities and dimensions in your case. You create random samples of states using importance distribution and use them to estimate the density of your state.
You are kind of correct by saying that these samples are points of pixels in the frame. However, there exist more information in the samples. There are also velocities and dimensions inside those samples. The algorithm eliminate less possible samples and move forward with more possible ones. In this way, you obtain most probable pixels that object could exist. You also obtain most probable velocity and dimension for the object.


Overlap Add - which length to use? I am doing realtime audio. Let's say I have a signal $x(n)$ which is passed in blocks of 512 (siobuffer) to the program filter/FFT method. Within this block I want to convolve the signal with the filter $u(n)$ (impulse response) of length 8192K.
So I zero pad the asio block from 513 to 8192K, do FFT of both arrays and multiply the results. After IFFT, I get the convolved signal $y(n)$ of length 8192K which I want to pass to the output asio buffer (lets say of length 512, too).

Now comes my problem: how much overlap-add to use? From my understanding, I must overlapp add the whole 8192K to the last byte, no? This means I divide $y(n)$ of 8192 by 512 into 16 blocks and subsequently overlapp-add all those 512 blocks to the next processed input blocks. When I read about overlap add in general, there are always examples where people have outputs $y(n)$ just a little longer than the input, so they normally overlap add approximately 1.5 of a imputblock.

Is there a mistake in my thoughts? Do I just need to overlapp add parts of $y(n)$, or must I overlap add all my 16 blocks each time?

You have a bunch of options here

You can zero pad both impulse response and signal to 16k, FFT and multiply. That's very inefficient though. Keep in mind that the sum of length of filter and signal need to be equal or smaller then the FFT size, hence an 8k FFT doesn't work. This gets a little better if you can truncate the 8k impulse response down to 7.5k. Then you can use an 8k FFT since 7.5k+512 <= 8k
You can keep accumulating 512k buffers until you have 8k of signal. Then zero pad both to 16k, FFT and multiply. That's a lot cheaper but also increases the latency substantially and you need to do the processing in a different thread.

Typically the best option here is a Segmented Overlap Add or Block Convolver. This works roughly like this

Break your 8k impulse response into 16 chunks of 512 each. Zero pad to 1024 and FFT. Now you have 16 frequency domain transfer function vectors $H_k(z), k= 0,1,2...15$
At every frame n, you get a new signal block $x_n(t)$ of length 512. Zero pad to 1024 and FFT. Now have the frequency domain signal vector $X_n(z)$. Make sure you keep the last 15 frames around as well, so you have $X_n(z), X_{n-1}(z) ... X_{n-15}(z)$
Now multiply the signal vectors with the transfer function vectors and sum up all the results. $Y_n(z) = X_n(z) \cdot H_0(z) + X_{n-1}(z) \cdot H_1(z)+ ... +X_{n-15}(z) \cdot H_{15}(z)$
Inverse FFT, you get 1024 time domain samples, $y_n(t)$
Manage the overlap. Create the output as the first 512 samples of $y_n(t)$ of plus the last 512 from the previous frame $y_{n-1}(t)$. Keep the last 512 samples from the current frame, $y_n(t)$, as overlap for the next frame.

This method is very efficient since you only have to do two 1k FFTs per frame and it also keeps the latency low. There are ways to improve this further by leveraging that the input data is real and not complex, but that should probably be a different question.

If you don't want to manage the overlap manually, use overlap-save instead.
https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method


Please help me analyse spectrogram I am working on noise reduction and I need to learn how to analyze spectrograms

I mixed a pure speech file with violet noise in audacity and got the following spectrograms:

Matlab version

Sonogram version: (enter link description here)

Spectrogram view in Audacity:

What do all the colours mean and why do they look different? Shouldn't everything have the same colours as in the Matlab version? Or are the exact colours not that important, since the darker pattern at the bottom of each figure still look similar in all of them?

--- EDIT :: Further additions: ---

Here is a screenshot of the spectrogram of noisy speech as well as output of my noise reduced sound file.

The one on the left is noisy signal while one on the right is the one passed through the noise reduction algorithm. Both figures are plotted through the spectrogram function of matlab, called as follows:

spectrogram(data, hanning(128), 64, 128, 16000)

The number 16000 is the Fs value (sampling frequency) returned by the wavread function used to read in the original noisy speech file into matlab. The sound file is a male adult speaking, mixed with violet noise in audacity

--- EDIT 2: ---

Also, here is the same thing in gray-scale, if it helps

The spectrograms look different because they use:

Different parameters for the FFT size and hop size (window overlap ratio).
They simply use different palettes. The spectrogram gives you an array of numbers, which are scaled and map to a color palette to produce a color image. Matlab's colormaps are listed here - by default it uses "jet" ; so the time-frequency bins of lowest energy are colored in blue, and those of highest energy are colored red.


Preserving formants using Cepstrum I've implemented an harmonizer, using Stefan Bernsee's code, which works great. I figured out the logic of finding the root note, and the right harmony in the scale.

The only problem is that, my upper voices sound like someone who has taken too much helium !

After some research, I found that, that was because I didn't preserve formants.
The problem is I don't know how to do that ! I've heard that I must use a Cepstrum, but that's pretty much it.

Stephan's code, modified by me is :

/* ***************** PROCESSING ******************* */
        /* this does the actual pitch shifting */
        memset(mSynMagn, 0, mParent.mFftFrameSize * sizeof(float));
        memset(mSynFreq, 0, mParent.mFftFrameSize * sizeof(float));

        for (k = 0; k <= mParent.mFftFrameSize2; k++)
        {
            index = k * pitchShift;
            if (index <= mParent.mFftFrameSize2)
            {
                mSynMagn[index] = mParent.mAnaMagn[k];
                mSynFreq[index] = mParent.mAnaFreq[k] * pitchShift;
            }
        }

I think I must mess with the SynMagn array but I don't know how.
I think calculating the Cepstrum, is just going through mAnaMag and taking the log ? (which base ?)
Then when I have that Cepstrum, what should I do with it ?
What should those line become ?

    mSynMagn[index] = mParent.mAnaMagn[k];
    mSynFreq[index] = mParent.mAnaFreq[k] * pitchShift;

Do you need extract the spectral envelope from original and pitch shifted signal and then warp the shifted spectral to match with the original!

You can use cepstrum, LPC or a channel vocoder to extract the envelope!

The real cepstrum can be computed using this equation:

Rcepstrum=real(ifft(log(abs(fft(windowed_signal)))))

Now do you need low pass Rcepstrum

And finally apply FFT again to get the Spectral Envelope, you will get a smoothed version of the spectrum!

At this point you will need get the diff spectral shapes from the (Shifted-Original), this is your correction factor, now apply the correction factor in the shifted spectral.

Edit

Let me try show how it is done:

FFT1 = Original
FFT2 = Pitch shifted

Get the Real Ceptrum from both signals:

cepstrum1 = real(ifft(log(abs(FFT1))));
cepstrum2 = real(ifft(log(abs(FFT2))));

Low pass window from cepstrum:

Choose some index value to cut the cepstrum (50, 100, 150), get just the firsts indexes (0 to 50 for example) and put zero at all.

Extract the envelope from the slashed window

envelope1 = exp(2*real(fft(cep_cut1))); % spectral shape
envelope2 = exp(2*real(fft(cep_cut2))); % spectral shape

Find the correct factor from envelopes:

warp_factor=(envelope1/envelope2)

Apply the correction:

new=(FFT2.*warp_factor)

Back you signal to time domain again (not forget, do you need overlap and add)

New_Formant_signal = real(ifft(new))

My small tests some times give me poor results, maybe you can use some peak picking in the FFT and interpolate the results to give a nice envelope, I will try one day, bu I never have time lol

Of course if you are pitch shifting monophonic singnals you can use PSOLA to keep the formants Intact...


Log Distance-Normal Shadow Model and Friis Since I'm founding several but different results about this model, I would ask you if all it is OK with this formulas:

1) Friis Model --> $Pr_0 = Pt + Gt + Gr + 20\log_{10}\frac{\lambda}{4\pi d}$

2) Log Normal Shadow --> $Pr_L = Pr_0 - 10\eta\log_{10} d + \chi_{\sigma}$

where:

$Pr_0$ is the power received in the Friis Model as depicted above
$Pr_L$ is the power received in the Log Normal Shadow Model as depicted above
$Pt$ is the power transmission
$Gt$ and $Gr$ are the respective transmitter and receiver gains
$\lambda = c/f$, where $c$ is the light speed and $f$ is the frequency
$\eta$ is the Path Loss exponent
$\chi_{\sigma}$ is a Gaussian zero mean random variable with standard deviation $\sigma$.

I am a bit concerned about the signs and the $\log$ operations.

I verified all the steps and they are all correct.


Using Fast Convolution to apply HRIR to soundwave I have a dozen source sound waves; each is about one second duration. I want to put each one through a digital filter (I'm truncating it to 128 taps).

Currently I'm doing this with a bog standard convolution. It is chewing up a lot of CPU.

I would like to switch to fast convolution.

Could someone outline the steps?

I'm aware that the basic technique is iFFT ( FFT(src) * FFT(filter) ), but the source and the filter are different lengths.

So what I'm imagining doing is this:

For each (source, filter) pair:

Truncate each filter at 512 taps  and  FFT it to give 256 frequency bins.
Now step through the source at 128-sample increments using a window size of 512.

For each step,
    * multiply by a (Hanning?) window,
    * FFT
    * multiply FFT with FFT of filter (complex multiplication on each bin)
    * add to destination buffer
Divide amplitude of destination buffer by 4, as we did 4 x overlap

Questions I have at this point are:
1. Is the basic method right?
2. Do I need to zero-pad my filter?
3. Is 4x overlap appropriate / optimal?
4. is that windowing correct?  Is there any way to optimise out the windowing step?
5. how good will this output be compared with standard convolution?

π

EDIT: I've just noticed that dividing by 4 at the end is probably wrong, shouldn't it be using the area under the windowing function A_w?  So  multiplying everything by A_w / 4.  Would that be correct?

Google "overlap add". Quick outline

zero pad filter to 256 taps, do an FFT
Initialize and "overlap" buffer of length 128 taps to all zeros
Initialize an input pointer to 0
Loop starts here
Take 128 samples from your input starting at input pointer. Zero pad to 256, FFT
Multiply the FFTs
Inverse FFT, add overlap buffer to the first 128 samples, output these samples. Put the second 128 samples into the overlap buffer
Advance your input pointer by 128 samples and go back to step 5. Repeat until the input is all done

Some quick tips:

Overlap save is slightly faster but a little less intuitive.
If filters and signals are all real valued you can speed things up further by either implementing a real valued FFT that's based on an N/2 complex FFT or by combining filters and/or signals in pairs
128 is probably close to the break even point between FFT based and direct FIR filtering. I wouldn't expect huge gains. Depending hardware and level of code optimization it may actually be slower.


Correlation between 2 functions I have a periodical function $f(x)$ that is sampled in time, at a rate of 1 millisecond.

Together with this, I also have a second function $t(x)$.

I want to find the correlations between periods of the first function and the second function $t(x)$.

How can I approach this issue?

If you want to merely compare the periods of the functions $f_1$ and $f_2$, you could use the following "metric" (it might not technically be one if it disobeys the triangle inequality, but someone else will have to check for that):

$\Delta (f_1, f_2) \equiv \dfrac{\left| T_1 - T_2 \right|}{\max {T_1, T_2}}$

Then the similarity can simply be defined as $s(f_1, f_2) \equiv 1-\Delta(f_1, f_2)$


FFT Overlap in Matlab I am getting incorrect amplitude of my signal post FFT (90% error) and it does not improve after applying any of the window functions prior to FFT. After reading a number of articles, I believe FFT overlapping might improve the accuracy. However, I am not sure how to do it in Matlab. Please could anyone shed some light on that. I tried using spectrogram but I keep getting error on input X to be double precision vector.

Here is an example without performing FFT overlapping:

L = 1024; %signal length
Fs = 1000; %sampling frequency
t = (0:L-1)'/Fs;

NFFT = 2^nextpow2(L);
x = 11*sin(2*pi*250*t);

J = flattopwin(L);
JX = x(:).*J(:);

P= fft(JX,NFFT)/L;
M= 2*abs(P(1:NFFT/2+1));

f = Fs/2*linspace(0,1,NFFT/2+1);

plot(f,M)

The amplitude observed from the plot is approximately 2.4 at 250Hz which is much lesser than the actual amplitude, 11. Is there any reason for the significant reduction of the amplitude on FFT of a windowed signal?

Thank you.

Matlab normalizes in the inverse Fourier transform, so if you want to calculate the fft accurately, you have to normalize by the fft vector by the lenght of the fft itself.

This matlab script will demonstrate the normalization. There is a 1 seconds long 1000 Hz
pure sine wave with amplitude 1. We will sample it with fs = 6000 Hz, and calculate the spectrum to it in 60 points.

After plotting, we expect two pure spikes at 1000 and -1000 Hz with amplitude 0.5.

fs = 6e3;
T = 1;
f = 1000;
x = sin(2*pi*f*(0:1/fs:T-1/fs));
Nfft = 60;
X = fft(x,Nfft)./Nfft;
stem((0:Nfft-1)*fs/Nfft-fs/2,abs(fftshift(X)))
title('Normalized spectrum')
xlabel('frequency')


Forward and Backward Filtering (MATLAB's `filtfilt()`) I have an input sequence $x(n)$ and impulse sequence

\begin{equation}
h(n) = \begin{cases} 0 & \quad \text{when } n = 0 \\ .239052 & \quad \text{when } n = 1 \\1.3357\exp(-0.650286n) & \quad \text{when } n > 1\\ \end{cases}
\end{equation}

I did a convolution and obtained $y(n)=x(n)*h(n)$. I trimmed some of the data  in the end from $y(n)$ and then I flipped $h(n)$ and I did reverse filtering $w(n)=y(n)*h(-n)$.
When I researched about this I found filtfilt in matlab is an easier option but I am trying to implement this in C.
So my question is without trimming the first output sequence  $y(n)$ I am getting a drifted version of the result (there is a delay).

Is it trimming really necessary or am I doing something wrong here.
One more question is I can see some transients both at the start and end  of $w(n)$ sequence when I use filtfilt I don't see that transients.The reason is matlab is using some initial states based on the below paper to get rid of the transients.

Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992.

To summarize :

Is it trimming really necessary ?
PS the figure attached without trimming the the resultant data has a lag..Trimming the $y(n)$ sequence before backward filtering aligns the result.
how to determine the initial states for getting rid of transients.

My opinion is that the trimming isn't necessary.

The method in the paper is based on a simple idea, if we filter the signal $ x [n] $ forward and then backward, the result should be equal to the result of doing it backward and then forward.

The paper shows how to chose the initial conditions which minimize the $ {\ell}_{2} $ norm between the the different approaches.


The complexity of such function run in Matlab The below function is representing an algorithm, so how can I get its complexity?  I don't mean the time of running by using the tic .. toc, I mean how many operation (Additions and multiplications) are performed in this loop.

for times=1:m;
    for col=1:N;
        product(col)=abs(T(:,col)'*r_n);
    end
    [val,pos]=max(product);
    Aug_t=[Aug_t,T(:,pos)];
    T(:,pos)=zeros(M,1);
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*Yy;
    r_n=Yy-Aug_t*aug_y;
    pos_array(times)=pos;
end

Size of parameters, m = 256 , N = 256, T= [256,256] and M = [256,1]

You should know what each operator (i.e., *) and each called function (i.e. product) does.  Then add up those operations.  For instance, I'm pretty sure that in Matlab, the way you're building up Aug_t means that Aug_t' * Aug_t generates a vector dot-product, so for each element in Aug_t there's a multiply-accumulate.


Group delay of moving average filter I have implemented a moving average filter (in Python) where I fill a list with values and average them. When new values arrive the oldest will be deleted.

Now I am wondering how much delay I get with this filtering.

I am sampling a sensor signal (voltage) with 1000 Hz and I have read often that the group delay (for my understanding that is the delay?) for a moving average filter for M values (window size) is (M-1)/2.

I can't really understand it because:

Lets say I have the a window size of 5 and following measured values arriving in my empty moving average filter:

6x 1 Volt
6x 3 Volt

|1|1|1|1|1|1|3|3|3|3|3|3|

First I am waiting until my moving average filter is filled (until the fifth measured value).

Then I get these average values

|1|1|1|1|1| 1 |3|3|3|3|3|3|   1   [1,1,1,1,1]
|1|1|1|1|1|1| 3 |3|3|3|3|3|  7/5  [1,1,1,1,3] <- 3 Volt first time measured
|1|1|1|1|1|1|3| 3 |3|3|3|3|  9/5  [1,1,1,3,3]
|1|1|1|1|1|1|3|3| 3 |3|3|3| 11/5  [1,1,3,3,3]
|1|1|1|1|1|1|3|3|3| 3 |3|3| 13/5  [1,3,3,3,3]
|1|1|1|1|1|1|3|3|3|3| 3 |3| 15/5  [3,3,3,3,3] <- 3 Volt average

I would expect that my delay would be (M-1)/2 = (5-1)/2 = 2 measurements
but I get my 3 volts in my moving average filter just after 4 more measurements (so delay M-1).

I would be very happy if someone could explain how I can determine my delay.

"Group delay" isn't the delay between the change on the input and the first effect; it's the delay that a packet of oscillations of different frequencies experience.

In the case of a linear phase filter (and your moving average, its impulse response being symmetric, is linear phase) that is the "center" of the effects of a single impulse. So, it's the center between you seeing the first effect of your input change, and the last effect.


How to Find Shift Between Two Signals with Cross Correlation? I have calculated the cross-correlation function between to waveforms. I would like to find the shift between the two waveforms from the peak of the cross-correlation function. I am not really sure how to do that precisely. In this example, I have shifted signal 2 by 1000 Angstroms with respect to signal 1. I want to recover this shift from the CCF.

Well, All in there is a simple play with the indices of the Auto Correlation Function.

Have a look at my solution for Finding Reference Audio Signal in Test Audio Signal and Cropping Accordingly.
You'll find there full MATLAB Code for the task you mention.

Enjoy...


Can STFT (Short-time Fourier Transform ) be more useful than FFT for analyzing stationary signals under some circumstances? the nature of STFT is to be applied on non-stationary signals. for stationary signals STFT and FFT sounds exactly same to me. However, I was wondering if STFT can lead to better results compared to FFT (when applied on stationary signals), especially when optimized the STFT window function or other possible parameters?

In the continuous version of the short-time Fourier transform, if you choose a uniform unit-height window, it just boils down to the continuous Fourier transform.

In practice, for the discrete case, you can find STFT parameters that just give you the same information as that of the FFT. Hence, STFT can always be as useful as the FFT.

Since real signals are often finite length, they are not stationary in the purest sense. Even if the noise is "stationary", we only have access to few realizations.

For instance, the estimation of the parameter of an additive noise from a few realizations would probably be more robust in the STFT space than from the FFT.


Speed from Depth/Pressure sensors I have a pressure based Depth sensor. For reference, MS5837-02BA.
All depth measurements I've been able to read are very accurate (±1cm)
I tried experimenting to find out the speed of the body by differentiating the measured depth value.
The measurements are taken at a 40ms time period. This is the time the sensor takes for each measurement. I've used this as the rate of measurement and calculated the speed but the values are not as correct i feel.
Here are the plots for both measured and calculated variables.

All the measurements are in meters.
Would like to hear what exactly I'm doing wrong, whether its the rate of measurement or something else.
Edit: accuracy  ±0.1cm to ±1cm.
I verified the depth measurements using a measuring tape and it was accurate down to the cm. I couldn't measure whilst the body was moving so that may a factor
Edit 2:
Response to Hilmar's comment:
I've gone through what you've stated and my original depth data. So i only took my raw Depth measurements this time and discarded the speed calculations my microcontroller has done. I've calculated the derivative by finding the difference between the current and previous measurement and dividing it by the sample period. $(x[i]−x[i−1])/dt$
Here are the plots again but with the sample number changed to time

As you can see the plots are similar as earlier. And about your comment on the rise at t = 2000, the rise is not instantaneous but its actually gradual over the period of 1.32s . Have a look at this pic, note that the period in question is between t=2025 and t=2058, which means its a total period of 33*0.04=1.32s .

And the change in depth is noticeable but not at a very high rate. Note that the measurement is actually in the order of 0.1m or 1 decimeter. Hence the body had moved 0.4719m in a period of 1.32s which is just 0.35 m/s which my object is capable of. Other than this i'm uncertain where exactly i've made a blunder here. Others suggest that the high frequency noise may be a cause of this but would a simple LPF solve this problem?

All depth measurements I've been able to read are very accurate (±0.1cm)

That's highly unlikely. It seems that something is quantizing your data at pretty rough interval: maybe 2cm or so plus some extra noise. There are is also a spike at t = 2000 units ?) where the depth jumps from 0.1m to about 0.5m in one sample. Your sample interval is 40 milliseconds. 0.5m in 40ms corresponds to a speed of 12.5 m/s. Can your object really move that fast ?
Overall your data looks quite bad (noisy, quantized & spikey) and not like a real physical object tends to behave.

find out the speed of the body by differentiating the measured depth value

The second curve doesn't look the derivative of the first curve at all. I would expect a massive velocity peak of 12.5 m/s at t = 2000, but your velocity curve just shows a minor blip. Another big peak should be at t= 4000. I have no idea what you did there, but it doesn't appear to be differentiation.

This is a simple calculation of the "velocity" from the depth measurements using Excel. It seems more reasonable than the plot in the question.
Note that the data used was from Web Page Digitizer, so it won't be quite the same as the original.

In the future please label the X and Y axis of your graph with proper units. I have no idea what the unit of your velocity are supposed to be but it's unlikely to be in m/s.


Square wave FFT results show frequency that seems too low? I have created an open source plotting package for a low cost oscilloscope - see link here: GitHub Plotting Software

But during my testing I have one result that may be wrong: a noisy 7 kHz square wave shows an FFT peak at 6.25 kHz. Can this be correct, or am I missing something?

Here is a python code extract for my FFT

# stack FFT plot in this figure vertically
plt.subplot(grid_stack[grid_row, :])

y_heading = self.y_heading

# Reduce signal_data array size by a factor of n (used to reduce the frequency range)
n = self.n
updated_df_data = self.reduce_sample_array(self.df_data, n)

# update sample rate and array size
sample_rate = float(self.info_dict['rate'][1])/n
array_size = updated_df_data[y_heading].size

# Calculate y-axis magnitude scaled to same units as y-axis in heading2Use data (ie volts)
yf = 2/array_size * fft(updated_df_data[y_heading].values)

# Calculate x axis as frequency in Hz
x = fftfreq(array_size, 1 / float(sample_rate))
x_half = x[:x.size//2]

freq_units, freq_multiplier = rescale_frequency(x_half)
x_half = rescale_data(x_half, 10 ** freq_multiplier)

y_half = abs(yf)[:array_size//2]
engr_power_v = calculate_scale(y_half)
y_half = rescale_data(y_half, 10**engr_power_v)

plt.plot(x_half, y_half, color=self.iplot_colors[-1])

Sample plot with clean square wave that seems to be ok

Sample plot with noisy square wave that seems to have wrong peak value (6.25 instead of 7 kHz)

Any comments or suggestions are welcome as this is my first python signal analysis software project.

The frequency resolution of your FFT is determined by its length and the sampling rate.

In your second plot you show Fs = 10e6 and 3200 samples. Assuming you don't zero-pad these samples when you pass them to the FFT this gives you a frequency resolution of 3125Hz:

$$ 3215 = \frac{10 \times 10^6}{3200} $$

So at lower frequencies you have frequency bins as 0Hz, 3125Hz, 6250Hz and 9375Hz. The largest peak in the FFT will be in the bin closest to the signal fundamental frequency.

If you want better resolution you either need to capture more samples, or append a load of zeros to the end of your time series (zero-padding).

I would recommend reading about windowing for FFT analysis also, especially if you go down the zero-padding route.

As an aside: the reduce_sample_array function seems to decimate the time series by simply discarding samples. You really should be doing proper down-sampling here, to avoid aliasing.


Different results applying DFT (formula only) and FFT on Matlab? I have this function in time-domain and in frequency domain after fourier transform:
$$s_1(t) = (t-2)e^{-t}u(t-2) $$
$$S_1(f)=\frac{e^{-2(1+j2\pi f)}}{(1+j2\pi f)^2} $$
First I create a time vector and a frequency vector ranging from [0s,8s) and [-50hz,50hz) respectively and evaluate them on my functions, using these lines:
s1=(t1-2).*exp(-t1).*heaviside(t1-2);  %s1(t)
S1f=(exp(-2*(1+(2*pi*f1*1i)))./((1+(2*pi*f1*1i)).^2));  %S1(f)

What I get when plotting magnitude/phase is this:

After that, I wanted to compare those results using FFT command in Matlab, so I did this:
S1k=fft(s1);
figure % creates a figure
subplot(2,1,1)     %creates a grid of 2 plots in one figure, selecting the   stem as the first plot
stem(k1,abs(fftshift(S1k,2)),'red') %plots magnitude of S1f
title('Magnitude vs Frequency')
subplot(2,1,2)      %selects the phase plot as the second one in the grid
plot(k1,angle(fftshift(S1k,1)),'blue') %plots magnitude of S1f

title('Phase vs Frequency')
And the result, to my surprise, is this:

As you can tell, there are plenty of differences in both plots, though "shape" its similar at least in the magnitude plot, but with different values in the y-axis.
What can be the problem? I'm sure the fourier transform I did by hand is good, but yet results are different.
Why? Isn't FFT and DFT similar except for the speed of calculation?
Any hints will be appreciated.

The Fast Fourier Transform (FFT) is a class of algorithms that efficiently implement the Discrete Fourier Transform (DFT). However, what you computed by hand is the Continuous-Time Fourier Transform (CTFT), which is substantially different from the DFT. The DFT is applied to finite length discrete-time sequences, whereas the CTFT is applied to continuous-time functions of possibly infinite duration:

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\quad&\text{(DFT)}\\
X(j\omega)&=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\quad&\text{(CTFT)}\end{align}$$

However, you can use the DFT to approximately compute the CTFT, but there will always be at least one of the following two errors (and usually both):

Truncation error: if the signal is too long (or has infinite duration as in your example), it must be truncated to fit inside the DFT window.
Aliasing error: the continuous-time signal must be sampled, and sampling will usually introduce aliasing unless the signal is band-limited and the sampling theorem is satisfied.

Since the signal in your example is neither band-limited nor time-limited you will get both errors mentioned above.

More information on approximating the CTFT by the DFT can be found in this answer.

The following Octave script shows how the approximation works for the example in your question.

% choose parameters (sampling frequency and time interval)
Fs = 50;        % sampling frequency
T1 = 2;         % lower and upper
T2 = 20;        % integration limits

Ts = 1/Fs;
DT = T2 - T1;
N = round(DT/Ts);   % DFT length

% sampled time-domain function on [T1,T2]
t = (0:N-1)*DT/(N-1) + T1;
xd = (t-2).*exp(-t);

% DFT approximation of CTFT
Xd = Ts * fft(xd);
fd = (0:N-1)*Fs/N;

% phase correction for offset T1
Xd = Xd .* exp(-1i*T1*2*pi*fd);

% exact expression for CTFT
N2 = round(N/2);
fc = fd(1:N2);
Hc = exp(-2*(1+1i*2*pi*fc))./(1+1i*2*pi*fc).^2;

% plot results
subplot(2,1,1), plot(fc,20*log10(abs(Xd(1:N2))),fc,20*log10(abs(Hc)))
    title('magnitudes (dB)'), xlabel('f'), legend('DFT','CTFT'),
    axis([0,Fs/2,-110,-10]), grid on
subplot(2,2,3), plot(fc,20*log10(abs(Xd(1:N2))) - 20*log10(abs(Hc)))
    title('magnitude error (dB)'), xlabel('f'), grid on
subplot(2,2,4), plot(fc,unwrap(angle(Xd(1:N2))) - unwrap(angle(Hc)))
    title('phase error (rad)'), xlabel('f'), grid on

The plot below shows the magnitudes of the CTFT and its approximation by the DFT, as well as the magnitude and phase approximation errors.