The Number of Sine and Cosine Waves in an $ N $ Point DFT This is bound to be an embarrassingly simple question, but here it goes...
I was reading the chapter on discrete Fourier transforms (DFT) of this really didactic online book, The Scientist and Engineer's Guide to Digital Signal Processing, by Steven W. Smith and I got stuck thinking about the last sentence in the caption to this illustration:

The parameter, $k$, determines the frequency of the wave. In an $N$
point DFT, $k$ takes on values between $0$ and
$\color{red}{\frac{N}{2}}.$

He goes on to say with regards to the $32$-point DFT,

Fig. 8-5

shows some of the $17$ sine and $17$ cosine waves used in an N = $32$ point DFT.

So it seems to indicate $16$ sine waves, or $(N/2)$, plus the constant component, and a parallel split for cosine waves.
However, the DFT is defined as:
$$X(k)\equiv \sum_{n=0}^{N-1} x(n)\; \exp\left(\mathbf-i\frac{2\pi}{N}kn\right), \quad k=0,1,\dots,\color{red}{N-1}$$
which is $\color{red}{N}$ values of $k$, not $\color{red}{\frac{N}{2}}.$
And given Euler's formula:
$$e^{-ix}=\cos x - i \sin x$$
it seems as though there should be an equal number of sine and cosine basis equations: one of each frequency parameter $k$ for a total of $N$ sine and $N$ cosine.
Of course, at some points, the sine or cosine components will be zero (as we move around the complex plane, and depending on $N$), but not necessarily $\frac{N}{2}$. For instance, the $4$-point DFT would indeed either fall on the $x$-axis (real - cosine), or $y$-axis (imaginary - sine), splitting the basis function in half. Yet this seems to be a very sweet and unique example.
What am I missing?

There is no contradiction.

Both cases you have information from $ -\frac{ Fs }{ 2 } $ to $ \frac{ Fs }{ 2 } $.

The representation using Cosine ans Sine functions is using Real Functions.
Those functions are symmetric as Real Valued functions.
Namely if you a coefficient for a Sine in a given frequency it is for sure have energy on the negative frequency as well (See the DFT / Fourier Transform of Sine and Cosine).
Namely using those functions all needed it to go from $ 0 $ to $ \frac{ Fs }{ 2 } $ as the negative value is implicit.

On the other end, for the Complex Exponent the energy is one sided.
Hence you must go from $ -\frac{ Fs }{ 2 } $ to $ \frac{ Fs }{ 2 } $.


Length of padding for multiple samples stored in a single MP3 file I have 20 short audio samples of a piano playing a single note. Due to internal limitations of my application, I need to concatenate and combine these into a single lossy file (either MP3 or AAC). This concatenated buffer will later be read, decoded to PCM, and split back into 20 buffers.

I realize that both VBR MP3 and AAC take advantage of various pyschoacoustic phenomena, including possible temporal masking. Is there a specific amount of silence that I must insert between the various samples in order to ensure that one does not affect the encoding of the other? Is this impossible to know due to the way encoders work?

I know that MP3 has a frame size of 1152 samples and AAC has a frame size of 1024 samples. Would 1152 samples of silence be enough?

The temporal masking effect is minimal. If your samples do not have extreme differences in loudness, one frame of silence should be enough. Depending on sampling rate one frame is >20ms, that should be plenty.


Region of Convergence for pole at z=-1 ? I am caught up in a following question:

I have a mixed (causal + anti-causal) input to a system whose output is causal, with system having transfer function with pole at z=-1 and zero at z=1. What is the ROC for the system (if any) ?

Thanks a bunch.

If the system is causal then the ROC is outside the pole with maximum magnitude. Since there is only one pole at $z=-1$, the ROC is $|z|>1$. Note that the system is not stable since the ROC does not include the unit circle.


Pesq numbers and clean speech comparison to processed speech I am trying to evaluate my processed speech for my DSP project. My project is an implementation of an MMSE algorithm proposed by Y.Ephraim and D. Malah, "Speech Enhancement Using a Minimum Mean Square Error Short-Time Spectral Amplitude Estimator," in the IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-32, No. 6, Dec. 1984.
As per project requirements, I ran PESQ, LLR, and segment SNR. I got numbers, but I am not sure what they mean. Is high number is good or bad? Also, what does it mean if I get a negative number?

Please help and thanks in advance

Well, I submitted my project. Thanks to the help of this group.

As for the PESQ number, I think the range is between 1 and 5, where 5 being excellent and 1 is low.

Anyhow, I did another trick to estimate my performance. Instead of feeding the functions (SESQ, LLR, and Segment SNR) the clean speech and the processed speech, I fed them clean speech with the clean speech. I got numbers and these were my reference (or maximum) numbers. Then I compared all my runs against these reference numbers


Pole-Zero Configuration of Glottal Signal I am currently studying for an exam in my Speech Processing course this Friday and I stumbled upon this question in one of the sample exams:

The answer to this question was given as:

Can someone explain to me why the pole is drawn in that location and why the reciprocal was taken?

Recall what a pole is? It's where the z-transform tends to infinity. What value(s) of $z$ would cause $G(z) =\frac{1}{(1-\beta z)^2}\rightarrow \infty$? (This will answer your reciprocal question).

Next, note that $\beta <1.$ So where would $\beta$ go relative to 1?
(This will answer your question about why the pole is drawn in the location shown.)


Convolution theorem with a kernel smaller than the image Suppose that I need to apply a convolution filter to an image. I want to do it by using the convolution theorem so I compute the kernel for the size of the input image to later calculate fft's and multiply. But I want the fft calculation to be less time consuming. The kernel loses its magnitude the further from its center, as usual. If I want to crop it, say, have a smaller kernel with only significant values, can I for instance compute its fft and later somehow extend the resulting frequency domain image of the kernel to fit the size of the input image, to later multiply them point-wise?

If you want to do the convolution in the frequency domain, you have to pad the kernel with zeros to make it the same size as the image.  There is really no way around it.  The upside is that you only need to take the FFT of the filter once, and then you can re-use it for multiple images.

If you are concerned with speed, you would have to do some analysis here. For small separable kernels doing the convolution in the spatial domain may be faster.


Ifft through Matrix multiplication I am still new to MATLAB, so apologies if I sound lazy to you. I am attempting to model a transformation as a set of matrix operations. I start with a vector, up-sample it by $U$ (up-sampling rate), take its Fourier transform, apply a Low-pass filter, take Inverse Fourier transform and finally down-sample at a rate $D$, where $U \neq D$.

I want to do all these transformations as a set of matrix multiplication, and not use MATLAB's built-in commands. Up-sampling / Down-sampling is easy. And so is Low Pass Filtering. But I am stuck with fft and ifft. I could use dftmtx to generate the matrix that would perform the fft, but using a vector $1000$ samples long, and up-sampling it by $1000$ makes a very big array, so that matrix generated by dftmtx becomes memory constrained. I have looked the code for dftmtx which is simply fft(eye(n)), where $n$ is length of the vector to be transformed. I tried replacing eye(n) with speye(n), so as to reduce the matrix values, but it would not work this way too.

The second problem is I do not have any idea how to do ifft using a matrix multiplication. I understand this matrix will also be a very large matrix and become memory constrained. However, no search on any online resource has helped me so far.

Whereas doing the above mentioned transformations is easily done by MATLAB, I need the matrix version for my application.

Any guidance is appreciated.

My code in simplest for is

n=1000;% length of original code
a=randn(n,1)+j*randn(n,1); %initialise
U=1001; %up-sampling rate
u=zeros(n*U,n); %Blank up-sample matrix
%find indices for up-sampling
for i=1:N;
    u(U*i-(U-1),i)=1;
end
ua=u*a; %up-sample a by rate of U
%Find Fourier transform
ft=dftmtx(length(ua));
fa=ft*ua;
%Designing LPF
lp=speye(length(fa)); %initialise
j=length(fa)/2-ceil(U/3):length(fa)/2+ceil(U/3); %define central region where high frequencies will be ignored
lp(j,j)=0; %ignore center values which have high frequency, render this region zero
lpa= lp*fa; %Low Pass Filter Up-sampled Fourier transformed sequence
%Inverse Fourier transform
ifa=ifft(lpa); %How to implement this as matrix multiplication so that we have ifa = ifft_matrix * lpa
%Down-sample
D=1000; %down-sample rate
j=length(ifa)/D;
d=zeros(j,length(ifa)); %Blank up-sample matrix %find indices for down-sampling code
for i=1:j
    d(i,D*i-(D-1))=1;
end
da=d*ifa; %calculate the down-sampled matrix da with length 1001 (original code stretched by 1 sample)

You can alternatively create a DFT matrix in matlab using this code:

exp(-1j*2*pi* ((0:N-1)/N).' * (0:N-1))

And the IDFT matrix thus:

1/N * exp(1j*2*pi* ((0:N-1)/N).' * (0:N-1))

As the only difference betweenm DFT and IDFT is the sign and a scaling factor.

You could alternatively just do:

ifft(eye(N))

But this doesn't get around needing the full DFT matrix, and honestly I don't see an obvious way of perfoming a DFT by matrix multiplication without actually having the DFT matrix. Perhaps if you were writing your own code you could leverage the symmetry of the DFT matrix so you only needed half of it but I don't think this idea will get you anywhere in Matlab.


Contributions to a DFT frequency bin So, a frequency bin is centered around an interesting frequency $f_0$. Assume the spectrum is such that there is no spectral leakage from other frequencies to that specific frequency $f_0$. Now. I center a frequency bin around f0. The value (magnitude) of this bin is then the value of EXACTLY this frequency $f_0$ and is not affected by non-zero frequency values (magnitudes) of frequencies that also lie in that bin's range. A frequency bin 'only' covers the value of its center frequency (+ potential spectral leaks to that center frequency).

Is that correct?

A frequency bin 'only' covers the value of its center frequency (+ potential spectral leaks to that center frequency). Is that correct?Is that correct?

No. All frequencies show up in all frequency bins EXCEPT frequencies that are an integer multiple of $f_{\Delta}=f_s/N_{FFT}$ sample rate divided by FFT length.

In other words a frequency bin is affected by ALL frequencies except the center frequencies of the other bins.

You need to take into account that the FFT assumes that the signal is periodic with N. For a periodic signal contains only integer multiples of the fundamental frequency.


DSP library compatible with QT? I have checked this and this Question in this site very much similar to my question. But I want to know are they Qt compatible? I want to start learn developing DSP application in Qt framework so I want some library compatible with Qt so that it can later be translated into Qtpython also? Can we use same libraries with Qt as it is? Will it make effect on changing data types to and from this library to QtObject? I am new to DSP and Qt. so may be the question seems silly.

Not quite the answer you requested, but I think the sanest thing(OOP wise)  to do in most of the situations is to have the backend of the application separated from the frontend.

You could create a standalone library(either static or dynamic) for your dsp code and separate it from the QT code.

This way your DSP code would be independent of the frontend you are using.
You could switch from QT to anything else(wxWidgets/GTK+ etc) and your DSP functions inside the lib would remain the same.


how can the channel fading model (Rayleigh, or Rice, or Log normal distribution) infer the way channel estimation is done? For a dense urban environment with rich scattering and no Line of sight (LoS)
path, the Rayleigh distribution can be used to model the fading.
Rician fading is a stochastic model for radio propagation anomaly caused by partial cancellation of a radio signal by itself.

How can the knowledge of the statistical distribution of the noise impact the way the channel estimation is done ?

I need to model a type of channel with the knowledge of its statiscal distribution: for that i use a Discrete FIR Filter block, and then impaired with the AWGN block.

My personal point of view: what matters much more than the distribution is the second-order statistics of the channel. In particular things like delay spread / coherence bandwidth and Doppler spread / coherence time. These parameters tell you a thing or two that you must take into account when estimating your channels:

The coherence bandwidth (roughly equal to the inverse of the delay spread) tells you for what bandwidth you can consider the channel approximately flat. If the bandwidth you want to transmit is below that, you can consider the channel a scalar quantity, which is simple to estimate. If your bandwidth is higher, you either need to go to frequency-domain techniques like OFDM and then estimate the channel on several frequencies (spaced no more than the coherence bandwidth) or estimate channel taps in time domain, e.g., when using a Rake receiver.
The coherence time (roughly equal to the inverse Doppler spread) tells you how quickly the channel changes in time, i.e., for how long you can consider it approximately constant. After the coherence time has passed, your channel estimate is no longer valid and you need to train again (or employ tracking if you can).

That said, the channel conditions (like LOS/NLOS) do influence these quantities as well. If you have a strong LOS, chances are your coherence time is higher since the LOS will not change rapidly in time.


Channel equalization affect on input signal I have the following problem.

Am trying to understand how the channel block affects the input signal s(n).
I know that x(n), the input signal to the filter is basically:
\begin{equation}
x\left(n\right)\:=\:s\left(n\right)* h\left(n\right)+v\left(n\right)
\end{equation}
However, the end goal is to find the Wiener filter tap weights: w0, w1, and w2.
So, in order to find that, I need to find the correlation matrix R of the filter tap inputs and the cross-correlation vector p relating the input of the filter to the output
My process to solving this problem was first by decomposing the channel block using partial fractions. As in the following:
\begin{equation}
H\left(z\right)\:=\:\frac{1}{1-1.2z^{-1}+0.35z^{-2}}=\frac{1}{\left(1-0.5z^{-1}\right)\left(1-0.7z^{-1}\right)}=-\frac{2.5}{1-0.5z^{-1}}+\frac{3.5}{1-0.7z^{-1}}
\end{equation}
This in turn, using the inverse z transform, would turn into the impulse response:
\begin{equation}
h\left(n\right)\:=\:-2.5\left(0.5\right)^n+3.5\left(0.7\right)^n
\end{equation}
Am stuck at this point because I don't know what to do next given the concept of convolution stated above.
The essence later is to find the correlation matrix R which should be in this form:
\begin{equation}
\left[R\right]\:=E\left[\begin{pmatrix}x\left(n\right)x\left(n\right)&x\left(n\right)x\left(n-1\right)&x\left(n\right)x\left(n-2\right)\\ \:x\left(n-1\right)x\left(n\right)&x\left(n-1\right)x\left(n-1\right)&x\left(n-1\right)x\left(n-2\right)\\ \:x\left(n-2\right)x\left(n\right)&x\left(n-2\right)x\left(n-1\right)&x\left(n-2\right)x\left(n-2\right)\end{pmatrix}\right]\:
\end{equation}
and also, to find the coss-correlation matrix p, which is in this form:
\begin{equation}
\left[p\right]\:=\:\begin{pmatrix}E\left(x\left(n\right)s\left(n\right)\right)\\ E\left(x\left(n-1\right)s\left(n\right)\right)\\ E\left(x\left(n-2\right)s\left(n\right)\right)\end{pmatrix}
\end{equation}
Please help me move on the right path on this problem. Thank you in advance.

This appears to be a homework problem so I will avoid providing the complete solution here, but the following should suggest a path toward a solution:
Typically this problem is done when the channel is NOT known and results in a least squared approach given the effects of the added noise $v(n)$. Solving for the equalizer coefficients when the channel is a known quantity results in a similar process to solving for inverse convolution using a matrix equation approach as is linked here: Compensating Loudspeaker frequency response in an audio signal with a convolution equation of the form:
$$h[n] * w[n] = \delta[0]$$
Where $h[n]$ represents the non-causal channel response and $w[n]$ represents the non-causal equalizer coefficients. I use "non-causal" to mean each is centered about $n=0$ which simplifies the equation to remove any delay in the result, as in $\delta[0]$ vs $\delta[n-m]$. The result in turn for the causal solution simply changes the indexing back such that the first sample starts at $n=0$ rather than the middle sample.
In basic English this is saying knowing $h[n]$, solve for $w[n]$ that would result in an impulse response when the two systems are cascaded, since the output and input of a system whose impulse response is an impulse would be identical.
This would result in $y[n] = s[n]$ when the noise is not added, and since the noise is white, this would also be the solution that any least squares approach using actual signals $s[n]$ and $x[n]$ would converge to. As constructed I do not see this as a least squares problem but simply an inverse convolution problem. This is similar to the use of least-squares versus zero-forcing equalizers and has similar construction: the least-squares solution converges to the zero-forcing solution as SNR increases.
In this case we solve for $w[n]$ using inverse convolution.


Recursive filter with truncated polynomial impulse response How can we recursively implement a causal discrete-time truncated infinite impulse response (TIIR) filter with an arbitrary truncated polynomial impulse response:

$$h[n] = \begin{cases}\displaystyle\sum_{k=0}^K c_k n^k&\text{if }0\le n<N,\\
0&\text{otherwise,}\end{cases}$$

where real $c_k$ are the arbitrary coefficients of the polynomial, nonnegative integer $K$ is the degree of the polynomial, $0^0 = 1$, and positive integer $N$ is the length of the part of the impulse response that may be non-zero?

Perhaps a recursive form allows time complexity to be not dependent on $N$, giving an efficient implementation when $K$ is small, even when $N$ is large. Moving average can be implemented recursively as an integrator with tail cancellation. That's a $K=0$ degree truncated polynomial impulse response. Perhaps higher degree polynomials are also possible.

A cascade of identical moving average filters gives a B-spline impulse response, which is a piece-wise polynomial of any desired degree but lacking arbitrary control of the polynomial coefficients. Perhaps some form of tail cancellation allows to isolate a single piece of a higher degree B-spline impulse response and to linearly combine the piece with the same of lower degree B-splines to get arbitrary polynomials of any degree. Or maybe there is another approach.

Some promising literature I found:
Oscar G. Ibarra-Manzano, Yuriy S. Shmaliy, "Implementation of Digital Unbiased FIR Filters with Polynomial Impulse Responses", Circuits Systems and Signal Processing, April 2012.

I don't think you can, I don't see any special property of the truncated polynomial representation that would help here.

The impulse response coefficients are related to the polynomial coefficients through a set of linear equations, so we can simply write this as a matrix multiplication.

$$h = M \cdot c$$
where $h$ is the impulse response vector, $c$ the polynomial coefficient vector and $M$ a matrix that looks something like
$$M = \begin{bmatrix}
1 & 0 & 0 & ...  \\
1 & 1 &1  & ...  \\
1 & 2 & 4 &  ... \\
1 & ... & ... & ...
\end{bmatrix}$$

That matrix appears to be invertible, so we simply express the polynomial coefficients as

$$c = M^{-1} \cdot h$$

for any given impulse response. That means that every FIR impulse response can be represented in truncated polynomial form and there is really nothing special about it. If there were a way to implement this recursively, it would mean that every FIR filter can be implemented recursively which, to the best of my knowledge, is not the case.


Determine the Z-Transform for the following sequence: $ |n|(\frac{1}{2})^{|n|} $ Determine the Z-Transform for the following sequence:

$$ |n|(\frac{1}{2})^{|n|} $$

I have tried to solve the above problem. However, the answer that I got is the negative of what is given in the solution manual. What I may have done wrong?

SOLUTION FROM SOLUTION MANUAL:

It's always good to do a sanity check on such results. E.g., you could try to see that $X(1)$ equals the sum of all time domain samples:

$$X(1)=\sum_{n=-\infty}^{\infty}x[n]\tag{1}$$

With $x[n]=|n|\left(\frac12\right)^{|n|}$ it is clear that the result of $(1)$ must be positive. However, for the solution from the manual you get $X(1)<0$, whereas for your solution you obtain $X(1)>0$. So I think you can be confident that your solution is the correct one.


When do we divide by $2\pi$ while solving Nyquist problems? I slightly lost the tempo in the course for a couple of days, so I apologies for the noob question if seemed so.
Why do we divide by $2\pi$ in Nyquist problems sometimes and sometimes not? I saw there is a frequency and a radial frequency. But I don't understand when do we divide max frequency by $2\pi$ and when do we not?
[Edit] for instance:
We have $cos(t\pi/3)$ for what T can we sample it with a dirac train, such as that we can recover the original signal afterwards?
On my lecture the max T is 3. I see the 3 but how the rest of $2\pi/3$ is gone?

Well!  all signals in this world are made up of sum of different rotations(sinusoidals) - different in three senses:
a. how big is the amplitude (A)
b. how fast is the rotation ($\omega$)
c. where is the starting point of the rotation (phase $\phi$)
Fourier made this very clear.
How do we measure rapidness of the rotating signals(sinusoidals) : by their angular velocity $\omega$ which is given in radians/seconds. This is the correct and most appropriate unit for measuring rotations - how much of angle is being covered per unit time! If a constituent signal is rotating by $200\pi$ radians in 1 second, we say $\omega = 200\pi$ rad/sec.
But we can also measure the rapidness of these rotating signals as number of rotations per unit time. That is $f$ expressed in $Hertz$ : number of rotations per second.
How are $f$ and $\omega$ related? One complete rotation covers a full circle meaning an angle of $360^o = 2\pi \ radians$. That means if the signal is making $200$ rotations per second then it is covering an angle of $200*2\pi = 400\pi$ per second. So, the relation between $f$ and $\omega$ is basically:
$$\omega = 2\pi f$$
That is why you divide $\omega$ by $2\pi$ when you want to express the frequency in $Hz$.
The example you have given is $\cos{\frac{\pi t}{3}} = \cos{\frac{2\pi t}{6}}$. The period of this sinusoidal is $T = 6$sec, therefore, the frequency in $Hz$ will be $f = \frac{1}{6}$ and in radians/sec will be $\omega = 2\pi f = \frac{2\pi}{6}$ rad/sec.
(Think about why the period of this sinusoidal is $6sec$. Figure out that the sinusoidal will repeat after $t = 6sec$. Figure out that the sinusoidal makes 1 full rotation of angle $2\pi$radians in $6sec$.)
And the sampling frequency is twice the frequency of the sinusoidal giving :
$$f_{sampling} = \frac{2}{6} = \frac{1}{3}sec$$
Meaning Sampling period is $\frac{1}{f_{sampling}} = 3sec$.


Why shouldn't the magnitude of a practically realizable frequency selective filter be constant for a range of frequencies? According to Paley-Wiener theorem, a practically realizable frequency selective filter shouldn't have a constant magnitude for a range of frequencies. Why is this so? How do we prove and understand it?!

  To summarize, causality has very important implications in the design
  of frequency-selective filters. These are: (a) the frequency response
  $H(\omega)$ cannot be zero, except at a finite set of points in
  frequency; (b) the magnitude $|H(\omega)|$ cannot be constant in any
  finite range of frequencies and the transition from passband to
  stopband cannot be infinitely sharp [this is a consequence of the
  Gibbs phenomenon, which results from the truncation of $h(n)$ to
  achieve causality]; and (c) the real and imaginary parts of
  $H(\omega)$ are interdependent and are related by the discrete Hilbert
  transform. As a consequence, the magnitude $|H(\omega)|$ and phase
  $\Theta(\omega)$ of $H(\omega)$ cannot be chosen arbitrarily.

Point b in the quote. Source: Proakis ch.10 4th edition

Note that all-pass filters with a perfectly constant magnitude response over the whole frequency range are realizable. A first-order discrete-time all-pass filter has the following transfer function:

$$H(z)=\frac{1-az^{-1}}{a^*-z^{-1}}\tag{1}$$

If $|a|>1$, the system $(1)$ is causal and stable. It is straightforward to show that

$$\left|H\left(e^{j\omega}\right)\right|=1,\qquad\forall\omega\tag{2}$$

i.e., the system's magnitude response is constant.

What is not possible is to have a constant magnitude response over only a part of the possible frequency range $[0,2\pi]$ (or, in continuous time, $(-\infty,\infty)$).

This can be seen as follows. First, design an ideal frequency-selective filter by computing the IDFT of some ideal frequency response (with a constant gain in the pass band(s), and a gain of zero in the stop band(s)). Let's consider discrete-time filters, but the discussion is basically the same for continuous-time filters. Let the resulting ideal impulse response be $h_{id}[n]$. First, we shift that impulse response to the right to have most of its energy at positive indices $n>0$:

$$\tilde{h}_{id}[n]=h_{id}[n-n_0],\qquad n_0>0\tag{3}$$

In order to obtain a causal filter, we need to multiply that shifted impulse response by a sequence $f[n]$ satisfying

$$f[n]=\begin{cases}0,&n<0\\1,&n>n_1\ge 0\end{cases}\tag{4}$$

with some non-negative $n_1$. The values of $f[n]$ between $n=0$ and $n=n_1$ can be chosen arbitrarily, but in practice we would choose them such that $f[n]$ increases monotonically in the interval $[0,n_1]$. The final impulse response of our causal filter is

$$h[n]=\tilde{h}_{id}[n]f[n]=h_{id}[n-n_0]f[n]\tag{5}$$

Note that $h[n]$ is still infinitely long.

In the frequency domain this corresponds to

$$H(e^{j\omega})=\left(H_{id}(e^{j\omega})e^{-jn_0\omega}\right)\star F(e^{j\omega})\tag{6}$$

where $\star$ denotes convolution. This convolution destroys the piecewise constant property of the original frequency response $H_{id}(e^{j\omega})$, and the magnitude of the resulting causal frequency response $H(e^{j\omega})$ is neither piecewise constant, nor does it have infinitely sharp transitions.


What are the steps for applying a 50% Overlap to an FFT? I'm currently working with audio where I apply an Hann window before running an FFT on a sample window of 4096 samples. I've been told that since I'm applying a Hann window I also need to overlap each FFT by 50%, however I can't seem to find a clear explanation of the steps involved in applying an overlap. Could anyone give me an idea of how to do so?

For a 50% overlap, offset successive 4096 sample windows by 2048.  e.g. the first half of the data in each FFT will be duplicated from the last half of the previous FFT window (before applying the Von Hann windowing function to either).


Solve equations with conjugates multiplications If I have a two variables $x_1$, $x_2$ and two equations. In the first equation the first variable $x_1$ is multiplied with the conjugate of same number which is multiplied with $x_2$, and in the second equation, the exact inverse happens.  For example, we have a number $1-2j$ so:
$(1-2j)x_1 + (1+2j)x_2 = y_1$  ... (1)
$(1+2j)x_1 + (1-2j)x_2 = y_2$  ... (2)
I'm confused, what's the proprieties of two above equations, I mean can we solve $x_2$ and $x_1$ if  $y_2$ and $y_2$ are known? .. If not, can we get the value of $y_2$ based on the equation (1), it means we still have two variables but with one equation, we need to get the second equation which is very similar to equation (1)

can we solve x2 and x1 if y2 and y2 are known?

Yes, provided that you mean $y_1$ and $y_2$ are known. If you have two linear equations your can (in most cases) solve for two unknown variables.

If not, can we get the value of y2 based on the equation (1),

No. $y_2$ doesn't show up in equation (1) at all, so you can't get it from there.

which is very similar to equation (1)

Similarity has very little to do with it. You can always solve for two unknown variables unless the two equations are linear dependent. The two that you have are not.


Is possible to compare to audio samples with cross-correlation Is it possible to compare two audio samples(meaning find word from one in another) realtime with cross-correlation?

No. If the words are pronounced with a different speed or intonation (prosody), or by different speakers, the cross-correlation won't have any significance.

A solution more robust to changes in prosody is to use audio features which are less dependent on pitch (such as MFCCs), and a comparison method that is robust to local time stretching (such as DTW or a distance between HMMs). It will still perform poorly when comparing a recording of two different persons saying the same word.


Answer to the first post

I guess that is pretty much dependent on the problem that you are dealing with. Boundaries are always a problem that needs extra care.

In most applications, it would be an option to set those values to zero (and handle normalization factors of the filter appropriately).

Other options would be to reflect the data, so that the index -1 becomes 1, -2 becomes 2, and so on. Alternatively, there can also be cases where it would make sense to assume that it is periodic, so the index -1 would become N, -2 would bekome N-1 etc.

Hints regarding your update and follow-up doubt

You are right, that does not look right. Honestly, I havn't had a closer look at the formula you provided. However, I wonder if they are correct. I doubt that, honestly. The sums over $m$ and $n$, respectively, are just going in a positive direction. Having said that, it becomes apparent that you won't get a proper result, eg for $C^\prime(1,1)$, since the indices of $I$ are always below 0. That does not look right. Instead of coming up with lots and lots of formula here, I would like to point you to two links:

Quite visual explaination of 2D convolution/filtering

The obligatory Wikipedia link ;)

Nice visualization of the process, yet with bad audio

I like to think about the process in a rather visual way, as described in the first link. In my mind, I place the middle point of the filter ($K$) on the pixel I would like to evaluate. The other filter entries are then over the adjacent image pixels. Everything outside the original filter dimension is zero. Then you multiply this "extended" filter pixelwise with the image, and sum all up. That is the filter value at the given $(x, y)$ coordinate. You then move the filter's mid point by one pixel position and do the same thing again. Repeat this process until you have covered the full input image. See link #3


Quantization SNR of bipolar signal I want to estimate the SNR of a bipolar signal. I know that each bit increases the SNR at 6 dB. Do i have to subtract one bit for sign, so that for bipolar signals the SNR is (n-1) x 6 dB ?

Your question appears to be about the signal to quantization noise ratio for uniform quantization. You have to count each bit, but that generally doesn't mean that the SNR is $6N$ dB, where $N$ is the number of bits. You normally have a constant offset:

$$\text{SNR}_{\text{dB}}\approx c+6N\tag{1}$$

where the constant $c$ depends on the Crest factor of the signal.

To give you an example, assume that the signal is in the range $[-x_p,x_p]$, and assume uniform quantization with $N$ bits. Then the quantization step size is

$$q=\frac{2x_p}{2^N}\tag{2}$$

If we further assume that the quantization error is uniformly distributed in the interval $[-q/2,q/2]$, the quantization noise power is

$$P_n=\frac{1}{q}\int_{-q/2}^{q/2}\xi^2d\xi=\frac{q^2}{12}=\frac13\frac{x_p^2}{2^{2N}}\tag{3}$$

With $x_{\text{RMS}}$ the RMS value of the signal, the signal to quantization noise ratio is given by

$$\text{SNR}=\frac{x^2_{\text{RMS}}}{P_n}=3\frac{x^2_{\text{RMS}}}{x^2_p}\cdot2^{2N}=\frac{3}{C^2_x}\cdot2^{2N}\tag{4}$$

where $C_x$ is the signal's Crest factor. From $(4)$, the SNR in dB is

$$\text{SNR}_{\text{dB}}=10\log_{10}\left(\frac{3}{C^2_x}\right)+20N\log_{10}2=10\log_{10}\left(\frac{3}{C^2_x}\right)+6.02\cdot N\tag{5}$$

which has the form of Eq. $(1)$.


Scalograms in python I am reading this paper to learning basic concepts of dsp and I want to reproduce the following scalogram of a test signal (fig 4.2 of the paper):

It has been produced from the discretization of the formula:
$$
\frac{1}{s} \int_{0}^{\Omega} f(t) \overline{g\left(\frac{t-\tau}{s}\right)} d t
$$
Where $g(t)$ is the wavelet:$$
g(t)=w^{-1} e^{-\pi(t / w)^{2}} e^{i 2 \pi \eta t / w}
$$
And $f(t)$ is the function:
$$
\begin{array}{l}
\sin \left(2 \pi \nu_{1} t\right) e^{-\pi[(t-0.2) / 0.1]^{10}} \\
\quad+\left[\sin \left(2 \pi \nu_{1} t\right)+2 \cos \left(2 \pi \nu_{2} t\right)\right] e^{-\pi[(t-0.5) / 0.1]^{10}} \\
\quad+\left[2 \sin \left(2 \pi \nu_{2} t\right)-\cos \left(2 \pi \nu_{3} t\right)\right] e^{-\pi[(t-0.8) / 0.1]^{10}}
\end{array}
$$
I have  tried to reproduce them using the function scipy.signal.ctw(). I have coded $g(t)$ as a function of t and $\omega $. Nevertheless, I get only non-zero coefficients for the last part of the wave, I guess that I could have something wrong. How can I compute the scalogram using scipy.signal.ctw()?

scipy's cwt is primitive and error prone; below is via ssqueezepy.cwt:

Code:
Note that if you seek to code the wavelet yourself, you'll need to take it to the frequency domain first (preferably analytically via Fourier transform and then sampling it, rather than via FFT), then pass it like wavelet = Wavelet(my_func); cwt(x, wavelet).
import numpy as np
from ssqueezepy import cwt
from ssqueezepy.visuals import plot, imshow

#%%# Helper fn + params #####################################################
def exp_am(t, offset):
    return np.exp(-pi*((t - offset) / .1)**10)

pi = np.pi
v1, v2, v3 = 64, 128, 32

#%%# Make `x` & plot #########################################################
t = np.linspace(0, 1, 2048, 1)
x = (np.sin(2*pi * v1 * t) * exp_am(t, .2) +
     (np.sin(2*pi * v1 * t) + 2*np.cos(2*pi * v2 * t)) * exp_am(t, .5)  +
     (2*np.sin(2*pi * v2 * t) - np.cos(2*pi * v3 * t)) * exp_am(t, .8))
plot(x, title="x(t) | t=[0, ..., 1], %s samples" % len(x), show=1)

#%%# Take CWT & plot #########################################################
Wx, scales = cwt(x, 'morlet')
imshow(Wx, yticks=scales, abs=1,
       title="abs(CWT) | Morlet wavelet",
       ylabel="scales", xlabel="samples")


Would a visible light beamforming "telescope" be feasible? If you had a phased array of high speed photodiodes that were spaced sufficiently far apart, would it be possible to apply beamforming theory to create a steerable visible light "telescope"?

Traditional, earth-bound telescopes are limited by the size of their primary mirror due to the mechanical properties of figured glass as the size increases.  Additionally, the larger the primary mirror the more robust the tracking mechanism needs to be to handle the weight.

It seems that if you could instead "virtually" track the object by beamforming the output of a horizontally mounted photodiode array, you are no longer limited by the size of a primary mirror and/or infrastructure to physically move the telescope.  Angular resolution could be increased simply by making the array larger.  Furthermore, your tracking accuracy is only dependent on your software processing and electronics, not moving parts.

Of course, the problem likely becomes dealing with immense quantity of data collected from the photodiode array, and the sheer physical size of the array.

Light travels at about 299,700 km/s in air, so ~29.97 cm/µs.  That seems within the realm of high-speed electronics given a large enough array (10x10m?)

Difficult, huge amount of data...but it doesn't seem impossible unless I'm misunderstanding how beamforming works?

Not in the way you think, no.

Digital Beamforming works because the digital receiver is able to sense the phase of the incident radio wave; that phase is retained.

Photodiodes don't sense the phase – they only sense the intensity.

Also, this requires the whole incident radio or light wave to be phase-coherent across its front – that's the case for laser light, but not for two random photons emitted from a huge fusion fireball, I'm afraid. So, the second requirement also falls flat. Generally, even non-laser sources of light have a coherency period – but the more different wavelengths contribute the smaller that gets.

In TX, that principle works – if you feed a laser through an array of digital phase modulators (optical crystals with electrically adjustable refractive indices), you can build phase-coherent multiemitters.

In RX, that might also work, if the received coherent light can be made to fall through a microscope objective, and then onto an acousto-optic deflector, you might, via modulating that deflector's state, be able to build a selective/phase-shifting summer, leading to beamforming. But: that's particle physics sized electronics normally, i.e. suitable for situations where you have a very strong laser illuminating very small points (in fact, I learned about AODs just today, and they're usually used the other way around: take one beam of laser, split it into multiple parallel beams, focus these, and trap quantums in the focus point. Cool stuff.).

Anyway, there's also mathematical accuracy limits to digital beamforming, and I'd must admit that I don't know whether they'd even allow higher precision with the kind of timing accuracy we can do today. Then again, phyisicists invent the coolest devices, and maybe you can


Using voice recognition algorithms to reduce noise in financial time-series? Is it possible to use voice recognition algorithms to reduce/remove the noise in financial time-series?
Does such an approach even make sense?

If "yes", how could that be done technically?

First, could you please explain what you mean by "noise" in financial time series? Are you talking about extracting a trend from time series?

In itself, speech recognition consists in converting a signal into a higher-level semantic, symbolic representation of this signal. It does not achieve any denoising, so it would be surprising to use speech recognition techniques to perform any modification of a signal.

Maybe you have something else in mind?

Anyway, a speech recognition system is made of the following layers:

A feature extraction module, converting chunks of audio signal into a vector of features capturing a representation of sound which is somewhat noise-robust and very roughly pitch invariant.
An acoustic model (usually, GMMs) computing the match (a likelihood score) between the audio feature vector and a very large set of phone models.
Optionally, an adaptation module "closing the loop" between the first two steps - learning a better parametrization that makes the extracted features conform more to the known model.
Then a hierarchy of HMM-like models (phones models, word pronunciation models, language model) - all combined into a single giant finite state transducer.

The first step doesn't make sense at all for financial data, since the parametrization used for speech recognition is deeply rooted into the way speech sounds are produced, and perceived by the human ear.

There could be financial applications of the kind of tools (Hidden Markov Models and their variants) used for the following steps... But then, it's a very well known artillery that is not specific to speech recognition and they are indeed widely used in time series analysis - for segmentation, change detection, detection of outliers, or automatic labeling of time series with the kind of "patterns" used by technical analysts.


Why Does ASK Modulation Create Fourier Sidebands? I know why analog amplitude modulation has side bands, it is related to (fc+fd) and (fc-fd). But what about DAM?
ASK(DAM) is a type of digital modulation, and there are only two state: carrier signal for "1" or nothing for "0". The expected behavior on a spectrum is seeing only one peak on the fc, but there are many side bands at fc+3d, fc-3fd, fc+5fd, fc-5fd etc.
Can someone explain that why we see these Fourier series pattern side bands?
fc: carrier frequency
fd: data signal frequency

If you did a continuous on off keying of a 10101010... pattern, then you would see sidebands as described since this is simply an up-conversion of the Fourier Transform of a 50% duty cycle square wave (moved to any carrier frequency). However if the data pattern for this case of a rectangular on-off keying was random, the resulting spectrum would be continuous and under the envelope of a Sinc function (Sinc squared as a power spectrum). The first case itself is also under the envelope of the same Sinc function in frequency, but since the repetition rate of the underlying pattern was not random (in this case the single pulse on off cycle) and repeated at a specific repetition rate, only integer harmonics of that repetition rate can be in the spectrum, with all the even harmonics falling on the nulls of the Sinc envelope (and thus a square wave with 50% duty cycle only has odd order harmonics).
This is demonstrated below, showing a single rectangular pulse and the magnitude (in dB) of its Fourier Transform. If this was a carrier frequency turned on and off, the same spectrum would result just at that particular carrier (DC is just another carrier frequency). When we repeat that pulse in time at a periodic rate, the spectrum can only have non-zero content at integer multiples of that periodic rate (just as implied with the Fourier Series Expansion), but will still be within the envelope of the base pulse (in this case a Sinc function).  So in the second case it is repeated at a 50% duty cycle, such that the even harmonics are in the nulls of the Sinc (the waveform repeats at $1/(2T)$). The third case is the same pulse duration repeating at a 25% duty cycle, so the 4th harmonic lands on the first null (with the waveform repeating at $1/(4T)$).

For example, if we generated a long pseudo-random sequence of such pulses with a pattern that repeats once per second, with each pulse having a duration of 1 millisecond, we would see spectral lines every 1 Hz under the envelope of a Sinc function with it's nulls at 1 KHz spacing. Thus we see in the extreme case of a random pattern that never repeats, that the spectrum will again be continuous such as in the first case of a single pulse (with the primary difference that it will be a power spectrum rather than an energy spectrum since if the pulse continues for all time, it will have infinite energy but finite power).


How to implement a digital sawtooth oscillator? In reference to this topic: How to implement a digital oscillator?, I would like to know how can I implement a solution similar to the one in the topic on a Sawtooth oscillator.

Refer to think link, and use the truncated phase word directly (no look-up-table needed):

Numerically Controlled Oscillator (NCO) for phasor implementation?


Why do we express complex noise as Gaussian distributed while the channel is Rayleigh distributed? I have a confusion regarding the distribution of noise and channel fading. As we know channel fading is usually modeled with Rayleigh fading since it has Gaussian distributed real and imaginary components. On the other hand, we express a complex noise as a Gaussian distributed signal (at least for AWGN channels) while it also has the same real and imaginary components.

Why do we express complex noise as Gaussian distributed while the channel is called Rayleigh distributed? Am I missing something here?

Thank you.

There are too many things named after Rayleigh floating around in this problem.

Consider first additive white Gaussian noise with two-sided power spectral density $\frac{N_0}{2}$.  This is generally taken as the model for channel noise, though in fact the source of the noise actually is thermal noise in the front-end of the receiver. NASA once even experimented with cooling the receiver front end in a liquid helium bath to reduce the noise temperature to about $5^\circ$K. Be that as is may,
the receiver has filters (and amplifiers) so that the noise at the receiver output is definitely band-limited. Because the front-end noise has been amplified greatly, the contribution of the impedances of later receiver stages to the thermal noise is negligible. Indeed, the noise at later stages is a
(stationary) band-pass Gaussian process $\left\{N(t)\colon -\infty < t < \infty\right\}$ where
$$N(t) = X(t)\cos(2\pi f_ct) - Y(t)\sin(2\pi f_ct)\tag{1}$$ with $\left\{X(t)\colon -\infty < t < \infty\right\}$ and $\left\{Y(t)\colon -\infty < t < \infty\right\}$ being independent zero-mean low-pass Gaussian processes with common autocorrelation function $R(\tau)$. And Yes, I do want a $-$ sign in $(1)$.

The autocorrelation function of the $\{N(t)\}$ process is
\begin{align}\require{cancel}R_N(
\tau) &= E[N(t)N(t+\tau)]\\
&= E\big[\big(X(t)\cos(2\pi f_ct) - Y(t)\sin(2\pi f_ct)\big)\\
&= \hspace{0.2in}\cdot
\big(X(t+\tau)\cos(2\pi f_c(t+\tau)) - Y(t+\tau)\sin(2\pi f_c(t+\tau))\big)\big]\\
&= R(\tau)\big[\cos(2\pi f_c t)\cos(2\pi f_c(t+\tau))+\sin(2\pi f_c t)\sin(2\pi f_c(t+\tau))\big]\\
&= \hspace{0.2in}~ - \cancelto{0}{R_{X,Y}(t)}\big[\cos(2\pi f_c t)\sin(2\pi f_c(t+\tau))-\sin(2\pi f_c t)\cos(2\pi f_c(t+\tau))\big]\\
\\&= R(\tau)\cos(2\pi f_c \tau)
\end{align} which is a band-pass function.  Note that we can also express the $\{N(t)\}$ process in terms of its complex baseband equivalent process $\{\mathcal Z(t)\colon -\infty < t < \infty\}$ where
$\mathcal Z(t) = X(t) + jY(t)$ as
$$N(t) = \mathrm{Re}\left(
\mathcal Z(t)e^{j2\pi f_c t}\right)\tag{2}$$ But, if we choose to represent the  complex number $\mathcal Z(t)$ in polar form as $B(t)e^{j\Theta(t)}$ where $B(t) \geq 0$ and $\Theta(t) \in [0,2\pi)$ are real numbers, then
$B(t) = \sqrt{X^2(t)+Y^2(t)}$ is a Rayleigh random variable and $\Theta(t)$ is uniformly distributed on $[0,2\pi)$. (This idea is the basis of the Box-Muller method of simulating Gaussian random variables. The usual probability transform method does not work well for Gaussian random variables because the Gaussian CDF does not have an inverse that is an elementary function whereas a Rayleigh CDF $1-\exp(-t^2/2\sigma^2)$ is readily invertible. So, the Box-Muller method calls a random number generator twice to simulate $B$ and $\Theta$ and then sets $X = B\cos(\Theta)$ and $Y = B\sin(\Theta)$ to get realizations of two independent Gaussian random variables).

The other place where Rayleigh's name comes in is in the transmission of signals over fading channels. Real-life signals are real-valued, not complex-valued.  When a sinusoid (say $\cos(2\pi f_ct), -\infty < t < \infty$) is transmitted over a (slowly fading) channel, the received signal is of form
\begin{align}r(t) &=  b(t)\cos(2\pi f_c t + \theta(t)) \tag{3}\\
&= b(t)\cos(\theta(t))\cos(2\pi f_c t) - b(t)\sin(\theta(t))\sin(2\pi f_c t) \tag{4}\end{align}
where $b(t)$ and $\theta(t)$ are slowly varying signals. For the purposes of digital communications, it is reasonable to treat $b(t)$ and $\theta(t)$ as constants over the duration of transmission of one bit/symbol, so that $(4)$ simplifies to
$$b\cos(\theta)\cos(2\pi f_c t) - b\sin(\theta)\sin(2\pi f_c t).\tag{5}$$
Slow fading occurs in ionospheric channels because the transmitted signal
is reflected back to earth from the ions in a ion cloud and the received signal is a sum of a very large number of very weak reflections with nearly equal delays. As a consequence, with a lot of handwaving and genuflections in the direction of the Central Limit Theorem, it is possible to regard
$b\cos(\theta)$ and $b\sin(\theta)$ in $(5)$ as realizations of two independent Gaussian random variables which allows us to treat $b$ and
$\theta$ as realizations of Rayleigh and uniform random variables respectively.  Hence, the complex baseband equivalent of a Rayleigh fading channel is that complex input $\mathrm{Re}(\mathcal Ae^{j2\pi f_c t})$ is received as output $\mathrm{Re}(\mathcal{AB}e^{j2\pi f_c t})$ where
$\mathcal B = Be^{j\Theta}$ with $B$ and $\Theta$ being independent
Rayleigh and uniform random variables respectively.  Fading is sometimes referred to as multiplicative noise since the fading multiplies  the input amplitude instead of adding to it as additive noise does.

"But, but, but," you say, "In a Rayleigh channel model, the fading parameters $B_{n+1}$ and $\Theta_{n+1}$ for the next bit/symbol interval are independent of the parameters $B_n$ and $\Theta_n$ for the current interval. How does that fit into your model of $b(t)$ and $\theta(t)$ being slowly varying functions? If $b(t)$ and $\theta(t)$ really vary so slowly that you can treat them as constant over the interval, then $B_{n+1}$ and $B_{n}$ should be very highly correlated, as should with $\Theta_{n+1}$ and $\Theta_n$."  Well, the dirty little secret that no one will ever tell you about is that communication systems intended for use over Rayleigh fading channels are very heavily interleaved so that what appears to be the next symbol in the mathematical model is actually a symbol that occurs $N\gg 1$ bit durations later where $N$ is chosen to be large enough that heavily-faded consecutive symbols appear in different frames rather than wiping out four bars of
Bethoven's Fifth Symphony or a whole stripe from an image of Old Glory if interleaving were not used.


Am I supposed to normalize FFT in Python? I'm sorry if this question is posed often, I can't seem to understand the answers already out there. I'm working with several forms of the Fourier transform, including the FFT, PSD, and spectrograms. I'm not sure whether, upon computing the Fourier transform of my signal, I'm supposed to normalize the result by some factor. For context I don't have any interest in comparing signals, I just need my outputs to be mathematically accurate.

I've seen some answers that say the DFT has to be normalized by a 1/N factor - I'm interpreting this as the 1/N factor in the IFFT, is that correct? Or is there some other reason you would have to normalize by 1/N?

What's really confusing me though, is answers where I'm told to divide by the sampling frequency. This answer implies that it is equivalent to dividing by N- does this have something to do with the transform being computed over 1-second segments or something along those lines?

A tutorial I'm working from normalizes this way, and I'm not sure where it comes from. They construct a whitening procedure like so:

def whiten(strain, interp_psd, dt):
    Nt = len(strain)
    freqs = np.fft.rfftfreq(Nt, dt)

    # whitening: transform to freq domain, divide by asd, then transform back,
    # taking care to get normalization right.
    hf = np.fft.rfft(strain)
    norm = 1./np.sqrt(1./(dt*2))
    white_hf = hf / np.sqrt(interp_psd(freqs)) * norm
    white_ht = np.fft.irfft(white_hf, n=Nt)
    return white_ht

where interp_psd is a power spectral density interpolant of the noise, and is computed in the following way:

NFFT = 4*fs
Pxx_H1, freqs = mlab.psd(strain_H1, Fs = fs, NFFT = NFFT)
Pxx_L1, freqs = mlab.psd(strain_L1, Fs = fs, NFFT = NFFT)

# We will use interpolations of the ASDs computed above for whitening:
psd_H1 = interp1d(freqs, Pxx_H1)
psd_L1 = interp1d(freqs, Pxx_L1)

Why do they multiply the whitened strain in the frequency domain by dt*2? I understand that 2 could come from discarding negative frequencies, but wouldn't the RFFT take that into account?

The tutorial later on computes the spectrogram of the filtered data:

# Plot the H1 spectrogram:
plt.figure(figsize=(10,6))
spec_H1, freqs, bins, im = plt.specgram(strain_H1[indxt], NFFT=NFFT, Fs=fs, window=window,
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
plt.xlabel('time (s) since '+str(tevent))
plt.ylabel('Frequency (Hz)')
plt.colorbar()
plt.axis([-deltat, deltat, 0, 2000])
plt.title('aLIGO H1 strain data near '+eventname)
plt.savefig(eventname+'_H1_spectrogram.'+plottype)

Here, they don't normalize, but I recognize that they might not need to since they're only interested in the visual of the spectrogram.

Later they normalize by the sampling frequency when performing a matched-filter exercise, but then reverse it.

# Take the Fourier Transform (FFT) of the data and the template (with dwindow)
    data_fft = np.fft.fft(data*dwindow) / fs

    # -- Interpolate to get the PSD values at the needed frequencies
    power_vec = np.interp(np.abs(datafreq), freqs, data_psd)

    # -- Calculate the matched filter output in the time domain:
    # Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin.
    # Taking the Inverse Fourier Transform (IFFT) of the filter output puts it back in the time domain,
    # so the result will be plotted as a function of time off-set between the template and the data:
    optimal = data_fft * template_fft.conjugate() / power_vec
    optimal_time = 2*np.fft.ifft(optimal)*fs

I apologize if this is too much information. All in all, my questions are these: do I have to normalize the output of a FFT in python (numpy, scipy, matplotlib) in order to be mathematically accurate, and by what factor? And does that normalization differ for transforms like the PSD or the spectrogram?

Your decision to normalize or not does not change the accuracy of your answer, as it is simply a scaling factor. If you use the common scaling of $1/N$, then the output for each DFT bin will represent the average of the portion of the input signal that is at the frequency defined by that bin, scaled to the same units as the input. So that is convenient and gives a certain meaning to what the magnitude of the output represents, but if you didn't divide by N it doesn't make the answer incorrect, as long as your approach is consistent.

For example, it is also common to use $1/\sqrt{N}$ as this would make the DFT and IDFT completely symmetric as the same scaling can be used for both.

To understand what is happening, consider the simplest case of the computation of the first bin of the DFT, which represents DC, and consider a DC signal that is 1 for all $N$ samples given. If this DC signal represents a magnitude in voltage, then this is samples of $f(t) = 1$ volt, given as $f[n]=1$, for all $N$ samples, $n = 0$ to $N-1$.

The first bin of the DFT without any normalization is simply the sum over all N samples:

$$F[k=0] = \sum_{n=0}^{N-1}v[n] = N$$

So if we were to scale this by $1/N$, then the result would represent the average for the frequency component of our signal within that bin in the same units as the input. In this case all values were at bin 0 "DC", so the average value was 1V, consistent with what we get when we scale by $1/N$. (It is clear from the formula above that when scaled by $1/N$ this is simply the average of $f[n]$).

Similarly all the other bins are the result of the average after multiplying by the time domain signal by $e^{-jk\omega_o n}$, which can be viewed as frequency translating any signal at that bin to DC and then taking the average of that result.


What do you mean by noise vairance and mean in an image? I know the sudden variation in the intensity level of an image is termed as noise. But I couldn't understand what do mean and variance of noise in an image refer to? What is the role of noise variance? As far as I have observed till now, most of the noise reduction algorithm's results vary depending on noise variance. Can someone help me in understanding the concept of noise mean and variance?

As I am very new to image processing, it would be better if someone could give a more elementary explanation. Thank you in advance!

Noise is not just a sudden variation in intensity; that could also be an edge of an object, a texture pattern etc. DC bias variations could also be considered as noise.

Noise is mostly an unrelated (independent of the signal) and random signal, most typically, imposed on the true signal of interest. Note that distortion (which depends on the signal) is something else.

For practical reasons, noise is mostly assumed to be white, and uncorrelated for simplicty but there are various different noise types as well, such as bandlimited (colored) noise, correlated noise etc.

For a WSS, zero mean, noise source, the variance gives its power. Stated mathematically: The variance of a WSS random process (noise) $v[n]$ is
$$ \sigma_v^2 = \mathcal{E} \{ (v[n]-\mu_x[n])^2 \} = \mathcal{E} \{ v[n]^2 \} - \mu_x^2 $$

When the process is zero mean; i.e. $\mathcal{E} \{ \mu_x[n]\} = 0$,  then
$$ \sigma_v^2 = \mathcal{E} \{ v[n]^2 \}  $$

Which is the definition of power of the WSS random process (noise) $v[n]$

Mean of ergodic noise is, its time (sample) average in the most practical sense. Or it's exptected value as ensemble average, in a more theroretical sense.

It's difficult to separate noise from true signal. A class of techniques to reduce noise, relies on statistical methods to characterize and differentiate noise and true signal, and adjust their filter characteristics according to noise variance (noise power) versus signal variance (signal power).

Optimal Wiener filters exemplify a class of linear (LTI) filters that compute minimum MSE expectations of true signals in noise, for given signall and noise statistical characteristics.


Is concept of "bit" in computer programming similar to the concept of "bit" in information theory? until today I knew that one bit is a variable, or a space in memory that can hold a value of either One (high) or Zero (low). This is the concept I learned from studying computer programming, microprocessor or DATA bus etc.

But after starting the course on information theory, I found out that bit is expressed as the information content of a symbol in message. This is calculated taking the logarithm (base 2) of the inverse of the probability of occurrence of the symbol.

Are these two concepts same ? On one hand one bit is a variable that can store either zero or one. On the other hand, one bit is the uncertainty associated with one of two symbols with probability of occurrence of 0.5. So, does 1 bit in computer programming or ASCII code mean 1 bit in information content of source or information theory?

A little edit:
here is one thing I am finding trouble understanding this topic. See, in data transfer of English alphabets, if we use ASCII code, we basically represent each symbol with 8 bits. Suppose that's 00000000 for a, 00000001 for b etc. So we are essentially allocating 8 quantization levels for each symbol.

But when the information theory comes into play, we take the probability of each symbol into account. 'E' has the highest frequency, where 'Z' has the lowest. So average information content comes down to 3 or 4 bits, right ?

My book says, 'Entropy or average information content is the minimum average number of bits required to represent each sample without distortion'.So, in this case, for efficient data transfer, are we creating maximum four quantization levels for each symbols? Because, on an average they carry information worth 4 bits. If that's so, isn't bit in information theory the same as the one in computer programming, data transfer or ASCII code etc ?

You probably get that I am clearly a noob here :p

They are not the same, but they're related. In particular, if you look at a computer memory holding $M$ "computer" bits, where each bit can be considered random and independent of all other bits, and there are roughly 50% of zeros, then the memory also holds roughly $M$ "information theory" bits.

Of course, this is often not the case: computer bits are usually correlated, and not uniformly random. This is why they can be compressed. Compressor programs such as LZW ("source coders" in information theory parlance) work, in a sense, by making each computer bit hold one information bit.

Edited to add: This example may make the distinction clearer. Consider a memoryless source with two outputs, $m_1=000$ and $m_2=001$, with probability 0.5 for each. Clearly, the information in each message is one (information) bit, but its length is three (computer) bits. A source coder, such as the Huffman algorithm, will readily code the messages to $c_1=0$ and $c_2=1$, compressing the source output. You can easily extrapolate this example to a source that produces ASCII-encoded text.

Note that, in the case of written languages in general and English in particular, nobody knows what the actual source entropy is, because there is no model for it. That is why there are contests for the best compression of large bodies of text; nobody is really sure what the optimum compression algorithm for English is.


How to add awgn noise to a signal in Fourier domain for given SNR I want to know the equations required to add awgn noise to a signal in Fourier domain for different SNR values. I know the procedure in time domain but not able to do the same in  frequency domain for the signal $x(t) = e^{2 \pi j f t}$ .

Thank you

Since the Fourier transform (and its discrete version) is linear, there is literally no difference in the way you do this in time and frequency domain:

You observe the energy in your signal, set the Gaussian's variance according to the desired SNR, and add the noise. Not a single flop was given about in which domain you were that day!


Wireless communications outage For communication systems, there is usually a target minimum rate $R$ bits per second. An outage occurs if the actual rate ever falls below $R$. Since $\frac{1}{2}\text{log}_2(1+\text{SNR})$, sometimes outage is also talked about in terms of $\text{SNR}$ and received signal power needing to be greater than some threshold $P$.

My question is: Say I have a very small packet of five symbols, and during the second symbol I experience an outage. The other four symbols received power was all greater than $P$. Does this mean that my entire packet was unable to be decoded?

EDIT:

For example, if I use the encoding scheme where: 0 = 00000, 1 = 11111, then I have 1 bit per 5 channel uses. So I can use the formula: $\frac{1}{2}\text{log}_2(1+\text{SNR})=\frac{1}{5}$ bits per channel use, and can solve for the required $\text{SNR}=0.3195$. Does this mean that all five symbols should have $\text{SNR}\geq0.3195$ in order to successfully receive the single bit?

No. It means 1 out of 5 symbols was not able to be decoded. Beyond that it depends on your encoding scheme and criteria for deciding if that packet is usable. For example what if you used an encoding scheme with BPSK symbols such that 1 1 1 1 1 = "1" and -1 -1 -1 -1 -1 = "0"?   (You could replace those with other codes as well 1 1 -1 -1 1 = "1" and -1 -1 1 1 -1 = "0" for example). Typically the packet is longer and error correction encoding is used so that individual symbol errors would not cause the entire packet to be lost.


Counting pixel connectivity groups How many 4-connected and 8-connected components are present in the binary image below?

Using definition of connectivities from the same textbook

"A connected component is a group of pixels that are connected in the sense that there is a connected path between any two pixels in the component. 4-connected means that the path comprises only horizontal or vertical neighbours (i.e. NSEW directions). In 8- connected regions the path may include diagonal neighbours as well (i.e. NW, SW, NE, SE directions)."

I get four 4-connected groups and three 8-connected groups.

The answer though is four 4-connected regions and one 8-connected region and it comes with the following diagram.

I realise I might be double counting some of the clusters, but I don't see how does the author get his/her answer from the above picture. Looking at it I get two 4-connected groups and one 8-connected groups. Unless the correct groups are grouped like this.

N-connectedness means that from any pixel to any other within the component, there should exist a path composed of N-connected steps. In other words, as long as two groups have one connected neighbor, they should be merged into a single group. It does not mean a 2 or 3 pixel group, but rather a component which consists of linked pixels.

See the connectedness garph theory.

With such definition, it is obvious that the author is right.


How to calculate bandwidth of a wire, "air" or other channels? I'm really curious of how engineers identify what is the frequency range that could be sent through "wire" as channel? Do they put wire in front of them and make test that consists of sending sinusoidal signals through it by gradually increasing frequency, till maybe wire suddenly "burns off" and they measure that as highest frequency of a "wire's bandwidth"? How do they test "air" as channel where EM waves travel?

First of all note that bandwidth may have multiple (but similar) definitions. A practically accepted one defines the bandwidth of a baseband lowpass channel (indeed a filter as an LTI system) as the frequency at which the magnitude of its frequency response falls to $1/\sqrt{2}$ of the value at the frequency $f=0$. This point is conventionally called as the $-3$ dB point, or as the cutoff frequency of the channel (filter).

A physical conducting wire can be electrically modeled (depending on the frequency of application) as an LTI system with relevant parameters of $R$, $L$, and $C$, according to its resistance, inductance and capacitance respectively. Such a model therefore will define an LTI filter with a frequency response magnitude which will have theoretical $-3$ dB point(s). In practice the bandwidth of the filter could indeed be measured by a set up similar to what you have defined.

The air as a channel means either the empty space (vacuum) or the atmosphere. The former case ideally has an infinite electromagnetic bandwidth which is practically limited by the ability to create EM waves and the physical space's intrinsinc space-time and energy-mass characterization. The latter has a much narrower bandwidth determined by the electromagnetic characterisation of atmospheric air composition and its different layers, furthermore the sperical shape of Earth and the conducting soil also affects the channel characterisation. It's well known that the channel bandwidth of atmosphere depends on frequency, especially the ionosphere layer.


which filters are these? Location of poles are given in unity circle and 4 zeros are given at origin for every plot,How to check which plot shows which filters?

For the Option (C) i am doing like this

I can write the Tranfer function from the plot like this

$H(z)=\frac{z^4}{Z^4-1}$

Now at z=1,$H(z)=\infty$

at z=j at mid frequecny,its its also coming $H(z)=\infty$

at z=-1,$H(z)=\infty$

How can i conclude that which filters it should be ?

What should be the general procedure to do this?

When only a rough pole-zero plot (without exact locations of them) is given for a filter, then the only procedure that remains to evaluate the frequency response magnitude will rely on the geometrical analysis which I shall show here for the case of a single zero and pole:

Given the Z-transform of a causal LTI filter which has a single zero at $r^{j\phi}$
$$ H(z) = (1 - r^{j\phi} z^{-1} ) $$

It's frequency response magnitude is
$ | H(e^{jw}) | = | 1 - r^{j\phi} e^{-j \omega} | $

In order to evaluate this magnitude geometrically, first treat the quantities as vectors on the complex-z plane, with real and imaginary parts making up their components.

Therefore the term $1 - r^{j\phi} e^{-j \omega}$ is a vector from the location of the zero $z = r^{j\phi} e^{-j \omega}$ to the unit circle for each frequency $\omega$ of analysis denoted by the vector $z = e^{j\omega}$.

The magnitude of this vector when $\omega$ spans a $2\pi$ interval from $0$ to $2\pi$ provides the magnitude of the Frequency response $|H(e^{j \omega})|$.

For this single zero filter, it's obvious from the vector that its magnitude will be smallest when $\omega$ points along the angle of $\phi$ and it will be a maximum when the opposite.

For the transfer function with a pole and zero of the form
$$ H(z) = \frac{ 1 - r_z e^{j\phi_z} z^{-1} } {1 - r_p e^{j\phi_p} z^{-1}  } $$

An similar procedure follows with the care taken to divide the magnitudes of the zero vector and pole vector at each frquency to yield to magnitude of the Frequency response.

Consider the geometric (rough) evaluation of the frequency response magnitude of a filter whose $4$ poles are as given in the plot of case A. There will be $4$ vectors from those $4$ poles to the unit-circle for each frequency of evaluation. I'll show $4$ frequencies $\omega = 0, \phi, \pi/2 , \pi$ respectively where $\omega=\phi$ is a frequency verl close to the pole angle locations. These four angle (frequency) points will let us varify rougly the frequency response magnitude of the BPF.

Each 4 vectors will have the following colors: $v1$ blue, $v2$ red, $v3$ green and $v4$ orange. And then we shall evaluate the frequency response magnitude as:

$$ | H(\omega)| = \frac{ 1}{ |v1| \cdot |v2| \cdot |v3| \cdot |v4| } $$
Based on the following graphics:

for $\omega=0$ (DC response): All four vectors are similar and nonzero in length and hence their reciprocal,$|H(\omega)|$, will be a very small number.

for $\omega=\phi$ (passband-response): Blue and Red vectors gets vanishingly small as $\omega \to \phi$ which means that their reciprocal, $|H(\omega)|$, will be very large at those frequencies justifying the bandpass charactheristics.

for $\omega=\pi/2$ (stob-band) similar to the DC case, none of the vectors gets small and all of them are nonzero large values (orange and green being larger) hence their reciprocal, $|H(\omega)|$, will be very small in line with the stop-band charactheristics.

for $\omega = \pi$ (stop-band): similar to first and third cases this time all four vectors are attaining their largest values hence making their reciprocal, $|H(\omega)|$, as the smallest value justifying the BPF filter characteristics.

From which we can conclude the following graph for $|H(e^{j\omega})|$:
(below is a MATLAB freqz plot for a $4$ pole filter with $r_1 = 0.95$, $\phi_1 = 0.9 \pi/3$, $r_2=0.95$ and $\phi_2 = 1.1 \pi/3$)


Frequency resolution and timestep in DFT I am learning DFT and I am not sure I understand it correctly. I am just going to write out the process of DFT and I am asking you to tell me if I am thinking the right way.

Suppose we have a signal on which we run DFT to get the frequencies incorporated in the signal. We use a window of size N samples and calculate DFT on it. The result are frequencies which occur in this chunk of audio. Now we can display the spectrum (like Winamp does). But we would like to display the spectrum while the song is playing thus calculating DFT for many such windows. Because it would take a long time to do this with all the windows, we say we have a timestep (dt) between windows. This way we loose some frequencies (not all samples are processed) but we gain in computing time and this way it is possible to display spectrum "realtime" while the audio is playing.

Again, I am not sure if what I wrote is correct so please shed a light.

Even if everything what I wrote above is true indeed, where does the frequency resolution (df) step in? My first answer would be the window size (N) because that way we can test more samples.

The "timestep" $dt$ (using your nomenclature) between DFT windows can really be anything you choose. There are a few different cases to consider:

$dt < N$: In this case, there is overlap between successive windows of samples that you input to the DFT. This can be used if you would like to get answers more frequently than every $N$ samples. It's quite common to see something like 50% overlap between windows being used.
$dt = N$: In this case, there is no overlap; you just take successive blocks of $N$ samples and use those as inputs to the DFT. This would give you the minimum processing load while not skipping any input data.
$dt > N$: As you noted, in this case, you're skipping data between windows (sometimes called "underlap"). This can reduce your processing load further, but since you're skipping data, you run the risk of missing the analysis of some feature in that portion of data that was discarded.

The frequency spacing of the DFT is related to the sample rate and the DFT length. Specifically:

$$
df = \frac{f_s}{N}
$$

where $f_s$ is the sample rate of the input data. $df$ in the equation above corresponds to the frequency spacing between each output (i.e. with $f_s$ measured in Hz, the outputs will have a spacing of $df$ Hz).


How to express linear convolution using positive frequencies from channel and symbol DFT vectors? Sorry for the long question, but I'm struggling to implement something related to this question. Your help would be appreciated.

Background

Notation: In the sequel, the uppercase letters represent frequency domain vectors, while lowercase represent time-domain.

Suppose I have two $N$-length hermitian symmetric DFT vectors $\mathbf{X}$ and $\mathbf{H}$, one representing a symbol and the other a channel through which the symbol goes. Thus, the received symbol in the time domain could be expressed as:

\begin{align}
\mathbf{y} = \mathbf{h} * \mathbf{x}
\end{align}
where $\mathbf{h}$ is the channel impulse response, $\mathbf{x}$ is the IDFT of $\mathbf{X}$ and "$*$" denotes the circular convolution.

This result could also be obtained by:

\begin{align}
\mathbf{y} = \tilde{H} \mathbf{x}
\end{align}
where $\tilde{H}$ is the circulant matrix.

Now, considering only the positive frequencies from both $\mathbf{X}$ and $\mathbf{H}$, that is, $\mathbf{X_+} = \{\mathbf{X(1), X(2), \cdots, X(N/2+1)}\}$ and $\mathbf{H_+} = \{\mathbf{H(1), H(2), \cdots, H(N/2+1)}\}$, it is possible to obtain $\mathbf{y_+}$:

\begin{align}
\mathbf{y_+} = \mathbf{h_+} * \mathbf{x_+}
\end{align}
where $\mathbf{x_+}$ is the IDFT of $\mathbf{X_+}$ and $\mathbf{h_+}$ the IDFT of $\mathbf{H_+}$, both having size $(N/2 + 1)$.

Then, $\mathbf{y}$ could be perfectly reconstruced from $\mathbf{y_+}$, since its IDFT:

\begin{align}
\mathbf{Y}_+ = \mathbf{H_+}\cdot\mathbf{X_+}
\end{align}

The $\mathbf{y_+}$ vector can also be obtained by making a circulant matrix with the IDFT of $\mathbf{H_+}$, which is a complex vector (because $\mathbf{H_+}$ does not have hermitian symmetry), and multiplying the new circulant matrix by the IDFT of $\mathbf{X_+}$, which is also complex.

Goal:

The ideal received symbol $\mathbf{y}$ is of my interest. However, I also need to know the intersymbol interference occurring when $\mathbf{x}$ is linearly convolved (not circularly) with the channel. I want to know this using $\mathbf{x_+}$ and $\mathbf{h_+}$, since I want to pre-equalize the symbols and using $(N/2 +1)$ vectors requires less operations.

Question

Using the IDFT of the positive frequencies from the original hermitian symmetric DFT vectors  as the actual transmit symbols, how could the result of a linear convolution be recovered perfectly?
The circular convolution is perfectly recovered.

I currently think it is not possible, because the interference will be different.

Sample code

The circular convolution is perfectly recovered:

N=8;

x = rand(N,1); % Time-domain random symbol
h = rand(N,1); % Time-domain random channel
X = fft(x); % Symbol FFT
H = fft(h); % Channel FFT
Y = H.*X; % Ideal Received Symbol (Ideal because of circular convolution)

% Time-domain complex vectors from positive frequencies
x_plus = ifft(X(1:(N/2 +1)));
h_plus = ifft(H(1:(N/2 +1)));
% Time-domain complex received symbol
y_plus = cconv(x_plus,h_plus,(N/2 +1));
Y_plus = fft(y_plus); % FFT of the complex received yields positive freqs.

% Comparison
[Y_plus; flipud(conj(Y_plus(2:end-1)))]
% Should be equal to:
Y

It doesn't work for linear convolution:

N=8;

x = rand(N,1); % Time-domain random symbol
h = rand(N,1); % Time-domain random channel
X = fft(x); % Symbol FFT
H = fft(h); % Channel FFT

y = conv(h,x); % Linear convolution
y = y(1:N); % Remove convolution tail
Y = fft(y); % Ideal Received Symbol (Ideal because of circular convolution)

% Time-domain complex vectors from positive frequencies
x_plus = ifft(X(1:(N/2 +1)));
h_plus = ifft(H(1:(N/2 +1)));
% Time-domain complex received symbol
y_plus = conv(x_plus,h_plus); % Linear convolution
y_plus = y_plus(1:(N/2+1)); % Remove convolution tail
Y_plus = fft(y_plus); % FFT of the complex received yields positive freqs.

% Comparison
[Y_plus; flipud(conj(Y_plus(2:end-1)))]
% Should be equal to:
Y

I think already your fist equation is potentially wrong. The output of the channel is given by LINEAR convolution signal and channel impulse response, not by circular convolution. Circular convolution can be used but only if signal and impulse response of the channel are sufficiently small and enough zero padding is applied to avoid wrap-around errors.

Dumping the imaginary part of the spectrum is also problematic. This basically an "analytic" signal (at least if the negative frequencies are set to 0). This typically is a non-causal operation and it becomes very difficult to determine when you will be getting wrap around errors and when you don't. You will almost certainly end up with a non-causal impulse response of the channel.

What exactly are you trying to achieve and why?


Why are LTI filters not able to output frequencies not present at the input? I know this a very basic question but I'm just having troubles understanding why.

Thanks in advance!

Because for an LTI system, every input sine has an output which is a sine of the same frequency, with a potential shift and an amplification factor.  In other words, input $x(t) = a\sin(\omega t+c)$ will be turned into some $y(t) = aH_\omega \sin(\omega t+d)$, where $H_\omega$ is the amplification factor for every sine  with frequency $\omega $.

This come from the theory: LTI are diagonalized by complex exponentials (in a pedantic way). As a result, an absent sine at a given frequency is a zero amplitude sine (in other words, $a=0$), and the output will be a shifted sine, with   zero amplitude whatever the amplification $H_\omega$  (since $aH_\omega=0$), still a zero amplitude sine.

More details:

Are complex exponentials the only eigenfunctions of LTI systems
Impulse Response to LTI


What is the difference between monocular and stereo image? I read somewhere in SE about stereo and mono videos. Stereo videos are captured by using two cameras with which we are allowed to measure depth accurately. But what about still images? what is difference between stereo image and mono one?

Another question: if I apply an algorithm on a stereo data set, the same result will be acquired as mono data set? to what extent segmentation of stereo images is different in comparison with other type of images?

Have you ever wondered why people and many animals have two eyes? By using two eyes we can have a perception of depth. The principle is very simple, while you are looking the screen keep your head fix and close one of your eyes, try to watch all the things behind your screen, now open your closed eye and close the other one and try to watch the same things, can you see the difference?

Each eye send information to our brain. This information is almost the same, if you could have two pictures of your both eyes and you try to overlap this pictures you will see that there are big displacements in the nearest objects to your eyes while the displacements in the objects that are far away will be very small. These observed displacements are translated to the depth of the objects.

Now other experiment that could help you to understand this concept is this: Try to keep your view fix close your eyes for a while, now open only one eye but try to don't move your view, try to touch the things around you, but remember keep your view fix, is it simple?. Now with one of your eyes closed, move your head and try to do the same, is it easier?. By moving your head, you are giving information to your brain of how far are the objects because in the nearest objects the observable displacements will be bigger. Of course you should consider that the eyes have also other kind of cells that help us to perceive depth. But in a simple stereo pair of cameras the information that we could recuperate is only the distance of the same point in each image, the disparity map, and there are a lot of ways to recuperate it, optical flow is very common. This kind of information is normally used to know the 3D space and the segmentation is applied comparing the the distance between the cameras and the objects.

While by using only a camera you can compare only the intensities changes and try to detect the main changes by using a lot of metrics. Is very common to use this kind of systems in tracking of objects by keeping the camera fix and trying to omit the things that aren't moving in the scene.


How to "scale" the FFT when using it to calculate discrete convolution? As you probably know, the discrete convolution $ H = F \ast G $ of some $ F \left[ x \right] $ and some $ G \left[ x \right] $ can be calculated using the Fast Fourier Transform (FFT).

To do this, one would ...

Set $ M $ to the length of the sequence $ F $ and $ L $ to the length of the sequence $ G $.
Choose an $ N = 2^n $ with $ n \in \mathbb{N} $ such that $ N \ge \left(M + L \right) - 1 $.
Extend both sequences $ F $ and $ G $ to a length of $ N $ by appending zeroes (on the "right"). (This is necessary to prevent time-domain aliasing from occuring due to the fact that we perform circular convolution instead of linear convolution.)
Calculate $ \hat{F} = FFT \left( F \right), \hat{G} = FFT \left( G \right) $.
Calculate $ \hat{H} = \hat{F} \odot \hat{G} $, where $ \odot $ is the Hadamard operator.
Calculate $ H = FFT^{-1} \left( \hat{H} \right) $.

Now, the thing is that there are different variants of the DFT (and therefore FFT) and they differ from each other regarding their scaling factors.

The "standard" version is as follows.

DFT: $$ X_{k} = \sum_{n = 0}^{N - 1} x_{n} \cdot e^{\frac{-i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

IDFT: $$ x_{n} = \frac{1}{N} \cdot \sum_{k = 0}^{N - 1} X_{k} \cdot e^{\frac{i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

However, there is also an "orthonormal" version of the transforms.

DFT: $$ X_{k} = \frac{1}{\sqrt{N}} \cdot \sum_{n = 0}^{N - 1} x_{n} \cdot e^{\frac{-i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

IDFT: $$ x_{n} = \frac{1}{\sqrt{N}} \cdot \sum_{k = 0}^{N - 1} X_{k} \cdot e^{\frac{i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

And, as far as I know, in principle, one could in fact use any scaling factors to transform back and forth between time and frequency domain, as long as the product of both is $ \frac{1}{N} $.

As long as you only do one transform, then perform linear manipulations in the frequency domain, then do an inverse transform, it does not really matter, which way you scale the transforms, however, when we calculate the convolution of two vectors by taking the DFT/FFT of both, multiplying them, and then taking the inverse transform, the scaling seems to matter.

When using the "standard" transforms, both $ F $ and $ G $ get scaled by $ 1 $ (not scaled at all) on the "forward" transform, then multiplied, then get scaled by $ \frac{1}{N} $ on the inverse transform. Thus, in total, you have one scaling by multiplication with $ \frac{1}{N} $.
When using the "orthonormal" transforms, both $ F $ and $ G $ get scaled by $ \frac{1}{\sqrt{N}} $. Then the transformed functions are multiplied, so that the product of the two effectively got scaled by $ \frac{1}{N} $. Then, on the inverse transform, they get scaled by $ \frac{1}{\sqrt{N}} $ again, so that in total, you got a scaling of $ \frac{1}{N \cdot \sqrt{N}} = N^{\frac{-3}{2}} \neq \frac{1}{N} $ so that the result will be different from that obtained using the "standard" transform.

Also note that, either way, we always scale by $ N $ (or some - possibly fractional - power of $ N $), which is neither (a power of) the length of the sequence $ F $ (which is $ M $), nor (a power of) the length of the sequence $ G $ (which is $ L $) but something rather unrelated to it. In the end, we chose $ N $ so that it fits our algorithm. Most FFT algorithms require the vectors to have a size which is a power of two, so we chose one. I'm not sure if (and how) one would have to compensate for that.

Summary: How should I scale the FFTs (or apply additional scaling before/afterwards to compensate) so that $ FFT^{-1} \left( FFT \left( F \right) \odot FFT \left( G \right) \right) \stackrel{!}{=} F \ast G $?

this formula $FFT^{-1} \left( FFT \left( F \right) \odot FFT \left( G \right) \right) \stackrel{!}{=} F \ast G$ is not always true, it's true when you scale at the inverse operation as you mentioned in your first case, in second case you have to multiply the right side by $\sqrt{N}$ to make it right.

Generally when you apply the Fourier transform over the convolution you see any scaling would be used for the first signal and no scaling remain for the second signal so if you multiply and divide by the scaling and use the scaling in the denominator for standard Fourier transform, the scaling remain in the nominator.


Signal amplification when summing low SNR images I'm trying to add up a few black & white images of the same object in order to amplify the SNR. I'm working on relatively low SNR images (about 1.1) and I wish to amplify it somehow.

My first try was to add up the images grey-levels while removing the noise through a simple threshold but this isn't really working since this threshold isn't quite well defined.

Are there any known algorithms for this kind of task? If so is there any pseudocode or even open source code somewhere?

Are the images properly aligned? Then you can usually add them all up and (for normalization) divide each pixel value by the amount of images that were added together. For $N$ images this will give you an SNR gain of

$$SNR_{+} = 10\cdot \log_{10}(N) \text{dB.}$$

This is under the assumption that the "signal" part of the images is always the same and that the noise in each image has the same power and is perfectly independent from the noise in other images. Hence, it is probably an upper bound on what you are actually going to achieve.


How to analyse the sound of Motor I have a motor, that I want to analyse, using its sound, so I can  its 'character' I was using the the simple fft, and figured out, that it wouldn't help, so I tried the order analysis, at least I incremented the speed of the motor and I can say it in spectrogram like here:

My question isn't about codes but about theory, what are the methods that I can use to analyse the sound of motor?
Thanks for any help !
link : related question

In my opinion, using the STFT (or spectrogram) as you already did is a good start. Try to change the parameters of your STFT, i.e. window size and frame shift until you are able to detect certain frequencies and their repetitions.

This is an example for a decelerating car on the highway. The inertial sensor (not a microphone) was placed at the axle:


Simplifying frequency response equation I am trying to to find a magnitude and phase response of a FIR filter as in problem 5.18 Digital Signal Processing 4th Edition by John Proakis. I got the frequency response as:

$$
H(\omega ) = 1 - e^{-j4\omega}
$$

I checked the manual solution and I got the same answer, but in that manual solution, the answer could be simplify to be:
$$
H(\omega)=2e^{-j2\omega}e^{j\pi/2} sin2\omega
$$

I don't understand what are the steps to simplify my answer into that solution, can someone derive it from my answer into the simplified solution. Thank you so much, hope that my question is clear enough

$$1 - e^{-4j\omega} = e^{-2j\omega}(e^{2j\omega}  - e^{-2j\omega}) \tag{1}$$

Now,

$$ \sin(2\omega) = \frac{e^{2j\omega}  - e^{-2j\omega}}{2j} \tag{2}$$

Equation 2 is a consequence of Euler's formula. Multiply and divide by $2j$ in (1) and use the identity (2) in equation 1 we have:

$$1 - e^{-4j\omega} = 2je^{-2j\omega}\sin(2\omega) \tag{3}$$

Now $j = e^{\frac{j\pi}{2}}$, substituting in 3 we have

$$1 - e^{-4j\omega} = 2e^{-2j\omega}e^{\frac{j\pi}{2}}\sin(2\omega) \tag{4}$$

The term $e^{-2j\omega}$ in (1) is represntative of the delay of the envelope in time domain, the term $(e^{2j\omega}  - e^{-2j\omega}) $ is reprsentative of the magnitude response and actual shape of the filter in frequency domain.

A general trick in simplifying expressions of FIR responses type is to tap expressions around the mid point of the impulse response taps, you could use this strategy to pair terms in responses with more than 2 terms as well


Explanation of LidarBoost Algorithm? I am trying to understand the LidarBoost algorithm as explained in this paper (PDF warning).

I don't understand how they take the original depth-images $Y_k$ and transform them into the up-sampled images $D_k$. I get how optical flow is used to align the $Y_k$ into a chosen reference frame, but is the transform between $Y_k$ and $D_k$ just a standard image up-sampling with nearest neighbors for the "interpolation" step? If that's the case, can someone explain how they get the term $W_k$ in their data term of the energy function? They say this about it, but I don't understand how to construct $W_k$ for each $k$:

  $W_k \in \mathbb{R}^{\beta_n \times \beta_m}$ is a banded matrix that
  encodes the positions of $D_k$ which one samples from during
  resampling on the high-resolution target grid

I'm trying (so far frustratingly ineffectively) to get enough of a sense of the algorithm to try to implement it.

It sounds like $W_k$ is just a matrix of 1's and 0's.  When $W_k$ is a 1, then the corresponding range value in $D_k$ is used. If it's a 0, then that range value is not used.

Similarly, $T_k$ should be mostly 1, and 0 when the corresponding $D_k$ values are unreliable.

Does that tally with your understanding?  If not, can you elaborate on your question with more details of what you suspect?


Pitch salience definition/references What is the pitch salience? I see it a lot when looking for papers about multiple pitch detection but haven't found a proper definition of it.

The perception of pitch, in particular pitch salience, is not a simple idea.  As recommended by Robert, if you can cite the paper where you saw the term pitch salience it would give us a chance to understand the context of where it was used and perhaps allow us to give you more succinct answers.

If you look at most textbooks on the subject of psychoacoustics (in this post I am referencing "Psychoacoustics, Facts and Models" by Fastl and Zwicker, Third Edition, Springer Publishing, 2007), you'll find descriptions of pitch perception for simple pure tones, complex tones, and bands of noise.  To give you an idea of the complex nature of pitch perception, the following paragraphs are brief summaries of the ideas found in Chapter 5 of the Fastl and Zwicker textbook.

As described in Section 5.1.1 of Fastl and Zwicker, the pitch of pure tones can be measured using different procedures.  In particular, measurements for "half pitch" and "double pitch".  In the "half pitch" case, a subject is presented with a pure tone of frequency f1 and has to adjust the frequency f1/2 of a second tone so that the second tone produces half the pitch of the first tone.  For example, if a pure tone of 440 Hz is used as sound 1 and a pure tone of variable frequency as sound 2, then the subject's average setting for the second tone is 220 Hz (halving of the pitch sensation corresponds to a ratio of 2:1 in frequency).  At higher frequencies, this relationship changes.  At an f1 frequency of 8 kHz, the "half pitch" frequency, or f1/2 frequency is 1,300 Hz.  Measurements at other frequencies above 1 kHz confirm this general observation:  for the perception of "half pitch", a ratio of the corresponding frequencies larger than 2:1 is necessary.  If the Fastl and Zwicker textbook is available to you, see Fig. 5.1.

The pitch of pure tones also depends on sound pressure level.  This effect can be examined by comparing the pitches of pure tones at different presentation levels.  The pitch of a pure tone at a presentation level L and frequency fL is measured by matching its pitch to that of a tone with level of 40 dB and a frequency of f40 dB.  For example, if a 200 Hz tone is alternatively presented at levels of 80 dB and 40 dB, the louder tone produces a lower pitch that the softer one.  On the other hand, if the same experiment is performed with pure tones at 6 kHz, the effect is reversed.  See Figure 5.3.

The pitch of complex tones is, well, a bit more complicated.

As described in Fastl and Zwicker, complex tones can be thought of as the sum of several pure tones.  If the frequencies of the pure tones are integer multiples of a fundamental frequency, then the complex tone is typically called an harmonic complex tone.  Like pure tones, the pitch of complex tones can be assessed by pitch matches with pure tones.  Although complex tones contain many pure tones, they do not produce many pitches, but only one single or perhaps one prominent pitch.  Here, the pitch of the complex tone is close to the pitch of the fundamental frequency or the differences between adjacent pure tones that make up the complex tone.  However, the relative frequency difference between the pitch of the complex tone and the pitch of the pure tone does depend on the fundamental frequency.  For example, the difference at 60 Hz is almost -3% (i.e. a pure tone, with a frequency of 58.2 Hz produces the same pitch as the harmonic complex tone with a 60 Hz fundamental frequency.  With a 400 Hz fundamental frequency, the relative frequency difference is about -1% (i.e. a pure tone of 396 Hz produces the same pitch as a complex tone with a 400 Hz fundamental frequency).  Above 1 kHz, the frequency of the pure tone and the fundamental frequency of the complex tone of the same pitch are equal.  As with pure tones, the pitch of harmonic complex tones is also influenced by presentation level.

Unlike the harmonic complex tones for which the component frequencies are integer multiples of the fundamental frequency, complex tones can also be produced by shifting the spectral components by some amount to produce an inharmonic complex tone.  With inharmonic complex tones, from which the lower components have been removed, ambiguous virtual pitches can be sometimes produced.  As HOTPAW2 alluded to, some complex tones have no discernible pitch and sound like noise.

If you are interested in the role that pitch salience has in music, read Trainor et. al (2014),"Explaining the high voice superiority effect in polyphonic music:  Evidence from cortical evoked potentials and peripheral auditory models", Hearing Research, Vol. 308, pages 60-70.  Further, see Degani et al. (2014), "A Pitch Salience Function Derived from Harmonic Frequency Deviations for Polyphonic Music Analysis", Proceedings of the 17th International Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5.

I hope this helps.

Michael.


Coherence for order spectra The scipy signal coherence function computes for two time series the coherence as a function of the frequency. For rotating systems however, it is sometimes meaningful to analyze the two recorded signal measurements as order spectra. For an intuitive difference between lines of constant frequency and constant order, I use Figure 2 from Order analysis in ArtemiS Suite - Head acoustics.

What is the most suitable adaptation of the scipy signal coherence in order to get the coherence as a function of the order, instead of the frequency?

The topic is to the best of my knowledge an area of open research. The research contribution which helped me most in solving my problem was: The spectral analysis of cyclo-non-stationary signals. Specifically, the key seems to be:

"In this paper, it is shown that the optimal way of describing
cyclo-non-stationary signals is jointly in the time and the angular
domains."


Z transform of finite signals I was trying to solve the Z-transform for u[n] - u[n-N], where u[n] means discrete unit step function, and N is some finite integer. I solved this using 2 methods.

Method 1 :

Taking z- transform using time-delay property and keeping in mind that delta[n] has z-transform = 1; I get:

which suggests that ROC is |Z| > 0

Method 2:

I know that z transform of $u[n] = \frac{z}{z-1}$ with ROC |Z| > 1

using this and the time-delay property on both u[n] and u[n-N] I say that:

which on making the denominators of the 2 fraction equal and simplifying becomes :

Which suggests ROC |z| > 1 .

The result in method-1 makes sense since the signal is a finite duration signal and taking z=0 would essentially mean a divide by zero situation while calculating the z transform.

But Method-2 is something that results from simply applying the properties of z-transform on some pre-known result for a special signal.

Why are the results different then?

Note that

$$1+z^{-1}+\ldots + z^{-(N-1)}=\sum_{n=0}^{N-1}z^{-n}=\frac{1-z^{-N}}{1-z^{-1}}\tag{1}$$

where I've used the formula for a finite geometric series.

So both your results are identical and correct.

The ROC is $|z|>0$, which is the case for all causal sequences of finite length. Note that in the expression on the right-hand side of $(1)$ there is a pole-zero cancellation at $z=1$, so in fact there is no pole at $z=1$, hence the ROC $|z|>0$.