short time energy function explanation Will someone be able to explain short time energy function for me in a simple way explaining all the variables in the formula?

The issue that I have is in the formula the summation varies m like,

m=n-(N-1) to n

from my research what I understood is that n is the center sample of the window, so if the center of the window is, lets say 80 and if the window size is 160. Then there is a situation where “m” get a negative value discontinuing a loop in a programming environment. Is my understanding incorrect or is there a reason? Thanks in advance.

If you're trying to identify birds, you should look into the Mel-frequency Cepstral Coefficients (MFCCs). They are commonly used in speech recognition, and there's actually a fair amount of academic papers on bird song recognition using the MFCCs and various other features.

Short-term energy is definitely useful to extract areas of bird sound, but it probably won't be sufficient to just take an FFT and look for coefficients. For example, imagine the same two  bird sounds, but one is slightly higher pitch, or there is more/different noise present. The FFT coefficients will be significantly different. The MFCCs basically capture the shape of the spectrum, on a scale that emphasizes regions of the spectrum where more detail is needed, and discard some of the unimportant details.


Data fusion using 2d discrete wavelet transform (DWT) I am working on a project that employs a linear array sensor that provides data from the same object at two different energies.  Collected in time, I end up with two images (16-bit sensor values, MxM image), call one HIGH the other LOW.  My engineering background is Material Science and metallurgy, so I am coming up to speed on DSP, image analysis, and all of the bits that bring all of that together.

I have been researching "image fusion" and have been working from Gonzalez's fabulous Digital Image Processing (3rd ed.) text.  I have developed an understanding of DWT and have successfully written code to replicate the 3-scale DWT transforms in the text (major accomplishment for me).

Now, I am trying to understand how to "fuse" the DWT coefficients from the LOW and HIGH images.  From the IEEE Xplore database, I pulled several papers on the topic; one of the more useful is by Zhiyu Chen, et.al. "A Combinational Approach to the Fusion, De-noising and Enhancement of Dual-Energy X-Ray Luggage Images".

Succinctly, they recommend "averaging the corresponding approximation coefficients of L and H" (sigma)... and "summing the corresponding detail coefficients of L and H" (psi).  But, that is all of what is said on the matter.

I have several other sources, but they must presume some obvious bit of understanding that I am lacking!

My question is, in order to successfully fuse the HIGH and LOW images: what does this mean?  In the case of a 3-scale DWT, do I average or sum the 3-rd scale approximation/detail coefficients?  Do I somehow average across each scale, despite the fact that the size of the coefficient matrices are different by powers of 2?

If I remember my wavelet transform correctly, this is what happens. Let's imagine that we are working with only the LOW image at the moment. I perform the first-scale decomposition which gives me the LOW_A0 (approximation) and the three detail coefficient images which are LOW_H0, LOW_V0 and LOW_D0. Then, to perform the second scale decomposition, you take the approximation image LOW_A0 and obtain LOW_A0_A1 (approximation) and the corresponding three detail coefficient images LOW_A0_H1, LOW_A0_V1 and LOW_A0_D1. Finally, in the last scale decomposition, you will split LOW_A0_A1 into its corresponding approximation coefficients LOW_A0_A1_A2 and the final detail coefficient images LOW_A0_A1_H2, LOW_A0_A1_V2 and LOW_A0_A1_D2.

So at the end you are left with the following:
1) The lowest approximation coefficients LOW_A0_A1_A2
2) The detail coefficients at each level which are 9 in number (LOW_H0, LOW_V0, LOW_D0, LOW_A0_H1, LOW_A0_V1, LOW_A0_D1, LOW_A0_A1_H2, LOW_A0_A1_V2 and LOW_A0_A1_D2).

Now, also consider that in parallel you have also performed the wavelet decomposition for the HIGH image and you have the same end stage products which are:
1) The lowest approximation coefficients HIGH_A0_A1_A2
2) The detail coefficients at each level which are 9 in number (HIGH_H0, HIGH_V0, HIGH_D0, HIGH_A0_H1, HIGH_A0_V1, HIGH_A0_D1, HIGH_A0_A1_H2, HIGH_A0_A1_V2 and HIGH_A0_A1_D2).

You want to now reconstruct a FUSED image which contains information from both the LOW and the HIGH images. At this stage, you would apply the fusion algorithm. Average the lowest approximation coefficients k*LOW_A0_A1_A2 + (1-k)*HIGH_A0_A1_A2, where k= 0.5 for equal weighting (you may want to weight one more than the other). Then, sum the detail coefficients at each of the corresponding LOW and HIGH detail coefficients. For example,
LOW_H0 + HIGH_H0.

Now take these approximation coefficient and detail coefficients and do the wavelet synthesis to get the final result. Does this make sense? I haven't read the paper/thesis you pointed out to though.


How much can the FFT noise affect the timbre analysis of a harmonic sound? The timbre of an harmonic sound (e.g. the steady part of a violin note) is defined by the amplitude of each harmonic (including the fundamental). In order to compare the timbre of two instruments (e.g. two violins playing the G string), I could compare those amplitudes. If the instrument are similar, I suppose that there will be a little difference in the harmonic distribution. Now, is the FFT a reliable tool for this kind of analysis? Should I take into account the noise that it introduces, or is it negligible?

An FFT is a finite length transform.  Therefore one does need to take into account any artifacts introduced by the window (implied rectangular, or optional Von Hann, et.al.).

There is some controversy, but you might consider sampling using a greater ADC bit depth than 16 bits at a rate higher than 44.1 ksps to feed you FFTs, and use one of the typical double precision FFT libraries to reduce any numerical noise.  (...and how you mic the violin may make a far greater difference than any of the above.)


Emulating guitar pickup position with delays If I have a vibrating string, then there will be no motion at the endpoints. If I move away a bit from an endpoint then I'll get a signal, but particularly the fundamental will still be attenuated. Higher modes may suffer less or even benefit. This explains why the bridge pickup (closer to the edge) of an electric guitar suppresses the fundamental more than the neck pickup and sounds "thinner".

The amplitude of the signal obviously depends on the position of the pickup and the frequency. The relationship is some $sin(\frac{\omega x}{c})$ function.

The question is: if I have a string signal, can I emulate this position dependant behavior with delays?

Suppose the pickup is at $\frac 14$ of the string. I can then add a delay with an amplification of $-1$ which reflects (sic!) the wave bouncing back from the near edge of the string. I will get an attenuated fundamental and all that, but is it physically correct?

Why not do the same with the far end of the string and use a delay which stands for $\frac 34$ of the string. It will also attenuate the fundamental, but the notches will be at different frequencies.

Now I realize, that I am about to model the string itself. If this model was done correctly, I could excite the string and pick it up anywhere I want. But I already have a string signal (from a Karplus-Strong algorithm).

What would I have to do, to add just the effect of the pickup position to this?
What wiring of delays would produce this  $sin(\frac{\omega x}{c})$ frequency response?

i won't divulge some proprietary information i have about how to do this, but i'll tell you this hint:

the vibrations on the guitar strings are standing waves.
standing waves as well as any waves on a string are made of of sinusoidal components with frequencies at harmonic intervals, that is each harmonic sinusoid has a frequency that is an integer multiple of the fundamental which is the frequency we will normally associate with the pitch of the note.
with standing waves, there are nodal points for each harmonic at distance that is a submultiple of the length of the string between the bridge on the right and the fret that is terminating the string on the left.
even-numbered harmonics will have a nodal point at the middle of the string and odd-numbered harmonics will not have a nodal point there.
there is a relationship between ration of the even-harmonic energy to the odd-harmonic energy and the position of the pickup and with the position of the pick.

have phun figgering it out.


Is this system linear? I have this system:

$$ y[n] = -\frac{1}{2}  x[n+2] - y[n+1] $$

I have no idea how to prove if the system is linear or not, because, it depends on future outputs...

Thanks for the help

For linearity it is irrelevant if the output depends on future values. This only affects the system's causality. For linearity you just need to check whether the response to the signal

$$ax_1[n] + bx_2[n]$$

equals

$$ay_1[n] + by_2[n]$$

where $y_1[n]$ and $y_2[n]$ are the individual responses to inputs $x_1[n]$ and $x_2[n]$, respectively.

And of course you can write the input/output relation as

$$y[n+1] = -\frac12 x[n+2] - y[n]\tag{1}$$

I will now characterize all sequences that satisfy the recursion (1) for the specific input signal $x[n]=\delta[n]$. Applying the $\mathcal{Z}$-transform to (1) we get (with $X(z)=1$)

$$Y(z)=-\frac12\frac{z^2}{1+z}\tag{2}$$

The function $Y(z)$ corresponds to one left-sided sequence and to one right-sided sequence. Note that left-sided and right-sided doesn't necessarily mean anti-causal and causal, but it only means that $y(n)=0$ for $n>N$ or for $n<N$, respectively, for some positive or negative integer $N$. These two sequences are given by

$$y_{p1}[n]=\frac12 (-1)^{n}u(n+1)\tag{3}$$
and
$$y_{p2}[n]=-\frac12 (-1)^n u(-n-2)\tag{4}$$

Note that neither (3) nor (4) is causal. The sequences $y_{p1}[n]$ and $y_{p2}[n]$ are particular solutions of the difference equation (1) for $x[n]=\delta[n]$. They are in fact impulse responses of two different linear time-invariant (LTI) systems described by the difference equation (1). The impulse response $y_{p1}[n]$ was mentioned by Dilip Sarwate in a comment below, and something similar to $y_{p2}[n]$ was mentioned by Robert B.-J. I completely agree with them that the difference equation (1) describes an LTI system, if it is used to describe a system with either (3) or (4) as impulse response.

However, the sequences (3) and (4) are not the only sequences that satisfy (1). Looking at the homogeneous difference equation

$$y[n+1]=-y[n]\tag{5}$$

we see that

$$y_h[n]=c(-1)^n,\quad c\in\mathbf{R}\tag{6}$$

is a solution of (5). Consequently, the general solution of (1) (with $x[n]=\delta[n]$) can be written as

$$y[n]=y_{p1}[n]+y_h[n]=\frac12 (-1)^{n}u(n+1) + c(-1)^n\tag{7}$$

Of course we could as well have used $y_{p2}[n]$ to define the general solution. Note that $y_{p2}[n]$ is obtained from (7) by choosing $c=-1/2$.

Now assume that we scale the input signal and use $x[n]=2\delta[n]$ instead. Obviously, the solutions $y_{p1}[n]$ and $y_{p2}[n]$ will change accordingly, but the general solution (7) will not (unless $c=0$ or $c=-1/2$), because the homogeneous solution (6) does not depend on $x[n]$. This is why I claim that the difference equation (1) does not necessarily describe an LTI system. In our case, it describes exactly two LTI systems and an infinite number of systems which are neither linear nor time-invariant.


How to practically correct for group delay imposed by a digital filter (Python) I'm using Python to filter a signal using a Butterworth filter (scipy.signal.butter) and subsequently passing this through a forward-backward filter with cascaded second order sections (scipy.signal.sosfiltfilt). Obviously, this introduces a group delay which I don't like.
I was told to use just the scipy.signal.filtfilt as it has a zero phase shift, though it is less numerically stable than the sosfiltfilt method.
I understand there are ways to calculate the group delay, for example scipy.signal.group_delay. This gives me an array of the phase shift for each frequency in my sample domain.
So this brings me to my general question - how is knowing the group delay useful? How can I practically  use the group delay at each frequency to correct for the phase shift in the filtered signal? Is there a method for this, or is it as tedious as delaying each frequency component separately and then summing them up to get a zero phase shifted signal?
Cheers

A causal filter such as what must be used in real time filters will always have group delay. Because it is unavoidable in such a circumstance it’s beneficial to know what it is.
Group delay cannot be corrected for in real time, as it would require knowledge about signal which hasn’t happened yet.  Variable group delay can be made constant by applying a filter with a phase response such that the sum of the applied filter and the incident filter is linear, but this will result in a constant, almost certainly non-zero and positive group delay.
In rendered applications, a method such as filtfilt will apply the filter, then flip the filter about zero time and apply it again, which will square the magnitude responds and result in no phase shift. There should not be anything about this process that results in stability problems.
Hypothetically, one could construct a zero phase version of a filter in a rendered application by way of some maths. One way would be to take the frequency transform of the resulting signal and subtract the phase response of the filter, and transform it back to the time domain.  This might have some gnarly edge artifacts.  You might try creating an all pass filter with the negative of the filter’s phase response, use that to filter the signal.  That would result in a signal that’s delayed by half the filter length, which you could remove manually. However, it sounds like you are using an IIR filter, so picking an appropriate length for prototyping the all pass might be tricky.


Gaussian Equivalent of Convolving an Image 50 Times with a Box Filter I have a image that I'm convolving 50 times with a box filter, and I like to replace this with a single gaussian filter. Because both box filter and gaussian filter are separable I can just study the 1D case.
Therefore...
Let $l \geq 1$ be fixed and $n = 50 \approx +\infty$ The box filter I'm using is defined as
$$
h_l(x) = \frac{1}{2l+1} \begin{cases}
1 & -l\leq x \leq l \\
0 & \text{otherwise}
\end{cases},
$$
here $x \in \mathbb{Z}$.
After a while I got that the gaussian equivalent for given $l,n$ (large $n$) is given by
$$
f_{n,l}(x) = \frac{1}{\sqrt{2\pi \sigma_{l,n}^2}} e^{-\frac{x^2}{2 \sigma_{l,n}^2}},
$$
where
$$
\sigma_{l,n} = \sqrt{\frac{n(l+1)l}{3}},
$$
the proof is in the answer of Fourier Transform of a Sequence and Inverse Fourier Transform.
Therefore I approximated the operator
$$
A_{h_l,n}={\underbrace{h_l*h_l*\ldots *h_l}_{\text{n times}} * (\cdot)}
$$
with
$$
G_{l,n} = f_{l,n}*(\cdot)
$$
However for practical purposes I had to consider the $f_{l,n}$ defined on a compact support (like $\left[-k\sigma,k\sigma \right], k \in \left\{1,2,3\right\}$.
The questions are:

Does this approach make sense?
If yes, and I've already implemented it for some reason using a single convolution with a gaussian (large kernel though) with 50 times the one with the box filter (you can assume $l = 5$ for my experiments) yields worse results, timewise speaking, assuming the approach is correct what can the reason be?

So either my approach is wrong or I'm probably missing something.
Just to clarify I'm trying to replace
for(int i = 0; i < n; ++i)
 blur(image,image,Size(2*l+1,2*l+1);

With
GaussianBlur(image,image,Size(2*k*sigma,2*k*sigma),sigma,sigma);

in opencv, and hopefully gain a speed-up.

What you should do is work with the Variance.
The Variance of Discrete Random Variables with support of $ \left[ -l, l \right] $ is given by (Notice that since the Expected Value is $ 0 $, the Variance is given by the 2nd moment)):
$$ \sum_{k = -l}^{l} \mathbb{P} \left( k \right) {k}^{2} = \sum_{k = -l}^{l} \frac{1}{2l + 1} {k}^{2} = \frac{1}{2l + 1} \sum_{k = -l}^{l} {k}^{2} $$
Convolution of a function which is the PDF of a random variable is equivalent of addition of the Random Variables (Assuming they are independent).
The Variance is linear with addition of Independent Variables hence the Variance of $ 50 $ Convolutions (Which is addition of 51 variables) is given by: $$ \frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2} $$
Now generate a Gaussian Blur with $$ \sigma = \sqrt{ \frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2} } $$ apply it on the image and it will be the best approximation for applying 50 times the Box Filter with Radius of $ l $.
By the Central Limit Theorem we can get better and better approximation as the number of Kernels gets higher.
Remarks

Usually Gaussian Blur with large STD is approximated by few (4-5) iterations of Box Blur. Yet in this case since there is a large number of convolutions the proper way is to calculate the STD of equivalent Gaussian Blur and then apply it either directly or by accelerated method (For example, approximation by few box blur iterations).
Box Blur can be applied in O(n) manner. Namely complexity is independent of the box kernel radius. It can be done by Integral Image or Running Sum.
Pay attention that if we do 50 Convolution of Box Blur with itself it is equivalent of 51 additions of Random Variables. Yet above the first convolution is done between the Box Blur Kernel and the images hence 50 Convolution means addition of 50 Random Variables.

MATALB Code
boxRadius   = 4;
boxGrid     = [-boxRadius:boxRadius].';
boxSupport  = (2 * boxRadius) + 1;
vB          = ones([boxSupport, 1]) / boxSupport;

boxStd = sqrt(sum((boxGrid .^ 2) .* vB));

vT = vB;

numConv = 9;

for ii = 1:numConv
    vT = conv2(vB, vT, 'full');
end

kernelRadius    = floor(size(vT, 1) / 2);
vKernelGrid     = [-kernelRadius:kernelRadius].';
kernelSupport   = (2 * kernelRadius) + 1;

kernelStd = sqrt(sum((vKernelGrid .^ 2) .* vT));

kernelStd / boxStd %<! Should be 'numConv + 1

vGaussianKernel = (1 / sqrt(2 * pi * kernelStd * kernelStd)) * exp(-(vKernelGrid .^ 2) / (2 * kernelStd * kernelStd));

figure();
plot(vKernelGrid, [vT, vGaussianKernel]);
xlabel('Kernel Support');
ylabel('Kernel Value');
legend({['Convolved Box Blur'], ['Gaussian Blur']});

Closed Form Formula for the Variance of Centered Uniform Discrete Random Variable
As seen above if $ x $ is a centered uniform Random Variable on the support $ \left\{ -l, -l + 1, \ldots, 0, 1, \ldots, l - 1, l \right\} $ then its variance is given by the second moment (As it is centered which means the mean value is $ 0 $):
$$ \operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2} $$
The above sum of squares can be solved either by integral or nice trick form the Binomial Identity:
$$ {\left( k - 1 \right)}^{3} = {k}^{3} - 3 {k}^{2} + 3 k - 1 \Rightarrow {k}^{3} - {\left( k - 1 \right)}^{3} = 3 {k}^{2} - 3 k + 1 $$
Summing each side from $ k = 1 $ until $ n $:
$$ \sum_{k = 1}^{n} \left( {k}^{3} - {\left( k - 1 \right)}^{3} \right) = {n}^{3} = 3 \sum_{k = 1}^{n} {k}^{2} - 3 \sum_{k = 1}^{n} k + \sum_{k = 1}^{n} 1 = \sum_{k = 1}^{n} {k}^{2} - 3 \frac{n \left( n + 1 \right)}{2} - n $$
Rearranging yields:
$$ \sum_{k = 1}^{n} = \frac{n \left( n + 1 \right) \left( 2 n + 1 \right)}{6} $$
Hence, for $ x $ as above:
$$ \operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2} = \frac{2 l \left( l + 1 \right) \left( 2 l + 1 \right)}{6 \left( 2 l + 1 \right)} = \frac{l \left( l + 1 \right)}{3} $$


Why are vector quantizers not used in video codecs? Using a vector quantizer (VQ), groups of image samples (or motion compensated and/or transformed data) can be quantized as vectors. I'm curious about why VQs are not preferred in the modern, highly popular video codecs such as H.264/AVC and H.265/HEVC. A quick Wikipedia search showed that there are some codecs (such as Daala) using VQs. However, the best codecs, in terms of compression efficiency, do not use them. Is that because of the difficulty of designing a good codebook and efficient searching of the codebook?

Theoretically, VQ is always better than scalar quantization (SQ). However, as you already mentioned, the complexity of VQ is higher than the complexity of SQ. It increases exponentially together with the dimensionality of vectors. Therefore, the question is how to design a practical coding system which can provide a best rate-distortion performance under a given complexity constraint.

A designing of a codebook is another very important issue in VQ. The main practical concerns about using a codebook in VQ are:

How to generate a universal codebook for a large number of images.
How to scale the codebook for various rate-distortion requirements.
How much storage space is reasonable for storing the codebook in memory.

All this practical constraints make it difficult for VQ to force through standardization process.


SINAD calculation in MATLAB I want to understand how SINAD is calculated. In MATLAB there is a funciton called SINAD and I would like to see if I can obtain "the same" answer as the inbuilt function does.
From the defintiion, I understand that it is the ratio of RMS of the signal to the noise and it takes all the values including the harmonics of the signal but the DC value.
So given the following code:
Fi = 2500;
Fs = 48e3;
N = 1024;
X = sin(2*pi*Fi/Fs*(1:N)); %Signal
Y= 0.001*randn(1,N); % Noise
x=X+Y;              %Signal + Noise
SINAD = sinad(x,Fs)

Matlab gives  me the value of SINAD for this particular signal. I believe that DC value is the first input to the signal, so if I wanted to calculate it "by hand" I would calculate the RMS of all the samples but the first one like this:
A=X([2:end]);
B=Y([2:end]);
z=sqrt(mean(A.^2));
w=sqrt(mean(B.^2));
sinad=20*(log10(z/w))

The answers are close, usually within 0.5 but are not the same. Is what I am doing correct? If not, how would I reproduce the SINAD function?

Is what I am doing correct?

Yes

If not, how would I reproduce the SINAD function?

By replicating the exact same algorithm that the Matlab function uses.
Your answer is more correct, but you also "cheated". You have access to the original sine and noise signals. sinad() ONLY has access to the sum of both. Hence it needs to do an estimation using windowing and Periodogramm analysis. That's less precise than calculating it directly from the original signals. If you really want to replicate what sinad() is doing, try an algorithm that only uses the combined noise+sine signal.


Why is X(0) the DC component Why exactly is  X(0) the DC component of a signal?

How is it equal to N times x(n)'s average value and why it is at X(0)?

Follows from the DFT definition. It's defined as

\begin{equation}
X(k) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi \frac{kn}{N}}
\end{equation}

So $X(0)$ is

\begin{equation}
X(0) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi \frac{0 \cdot n}{N}}
\end{equation}

Having $k=0$ gives $e^0=1$ all the time so that

\begin{equation}
X(0) = \sum_{n=0}^{N-1} x(n) 1
\end{equation}

Comparing this to the average

\begin{equation}
\overline{x} = \frac{1}{N} \sum_{n=0}^{N-1} x(n)
\end{equation}

shows that $X(0) = N \overline{x}$


FFT for a non coherent wave First of all, I am very new to signal processing.
I use FFT to get various characteristics (SNR, THD, etc.) of a sinusoidal wave, represented as a series of double numbers. However, this works only if the wave is "coherent", i.e. the last sample completes the wave. In other words, the number of periods is an integer number. If I just remove a few samples from the wave, I get a complete different FFT result, and the characteristics also look totally different.
Is this a known issue in FFT? If yes, how should such waves be analysed?
Thank you in advance.

Oh yes.  This is a known issue.  What you have just discovered is called "spectral leakage". Entire books could be written on this topic. What are you going to want to do is to use a "window" to reduce the effects of spectral leakage. I have some tutorial material on my website that may help you.

http://mechanicalvibration.com/windowed_DFT.html


Linear Convolution in DSP and Hann window My understanding is that when dealing with discrete finite sequences, convolution done thru FFT is circular. To obtain linear convolution, one would pad the input with zeros up some appropriate length. Then take FFT of two padded sequences, element-by-element multiplication, then IFFT.

For very long signal, to do Linear convolution, one might use the Overlap-Add (OLA) method, which breaks the long signal into segments of equal length, then do the linear convolution as above and sum all the overlapping section together to get the final output.

My question is that isn't that when one breaks a long signal into parts and to use FFT on it later, one should apply a Hann window to minimize Gibbs phenomena due to discontinuities that might occur at two end points ?? So far all the documents about OLA that I come across has not shown use of Hann window in conjuction with OLA. I wonder if it's somehow not needed or all the OLA document just want to illustrate one idea clearly, so as not to confuse reader with addition of Hann window. I understand that if we use both OLA and Hann, it will be more complicated. Please enlighten. Thank you.

This is an insightful question, one that I remember pondering in the 1980s when I first learned of the Overlap-Add (OLA) and Overlap-Save (OLS, sometimes called Overlap-Scrap) methods of "fast convolution".

It turns out that for strict LTI filtering, that you need not window with a Hann, but you could if you wanted to, and for a general frequency-domain modification algorithm, you may want to window with Hann anyway.

Let $w[n]$ be your window and have non-zero length $L-1$ which is less than the length of the FFT, $N$.  If it's Hann, the symmetrical definition (with $L$ even) is:

$$ w[n] = \begin{cases}
\tfrac12 \cos\left(\tfrac{2\pi}{L}n \right) + \tfrac12 \qquad & \text{if } |n| <  \tfrac{L}{2} \\
\\
0 \qquad & \text{if } |n| \ge \tfrac{L}{2} \\
\end{cases} $$

You can see that

$$ \sum\limits_{m=-\infty}^{\infty} w[n-m\tfrac{L}{2}] = 1 \qquad \forall n \in \mathbb{Z} \tag{1}$$

then

$$\begin{align}
 x[n] &= x[n] \sum\limits_{m=-\infty}^{\infty} w[n-m\tfrac{L}{2}] \\
\\
&= \sum\limits_{m=-\infty}^{\infty} x[n] w[n-m\tfrac{L}{2}] \\
\\
&= \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{L}{2}]\\
\end{align}$$

where $$ x_m[n] \triangleq x[n+m\tfrac{L}{2}] w[n] $$ is a "grain" of $x[n]$ slid over to be centered at $n=0$.

You can see you can take each segment or grain, $x_m[n-m\tfrac{L}{2}]$, slide it over to 0 to be $x_m[n]$, pass that through your filter (however it's done, with an FFT or not) to get $y_m[n]$ (which is not normally as short as $L-1$, but must be no longer than the FFT size, $N$), which you slide back to $y_m[n-m\tfrac{L}{2}]$ and add up each term to get $y[n]$:

$$ y[n] = \sum\limits_{m=-\infty}^{\infty} y_m[n-m\tfrac{L}{2}] $$

But this works whether the segments overlap or not as long as Eq. $(1)$ is valid.  You can make this work with a rectangular window.

But, for a rectangular window, the hop length (which was $\tfrac{L}{2}$ for the overlapping Hann) will be the same as the non-zero window length for a rectangular window (which was $L-1$ for the Hann).  So overlap-add will be more efficient with a rectangular window.  But you might want to do other frequency-domain things (like a phase vocoder) so I might recommend using the Hann window anyway for OLA.  There is something elegant about it.


Frequency modulation of $m(t)$ in MATLAB I have a message signal defined in the below interval
$$
m(t)=\begin{cases}1&, 0\leq t < \frac{t_0}{3}\\[2ex]
-2&, \frac{t_0}{3}\leq t <\frac{2t_0}{3}\\[2ex]
0&, \frac{2t_0}{3}\leq t < t_0\end{cases}
$$

where $t_0 = 0.15\ \rm seconds$, and sampling frequency $F_s=2\ \rm kHz$.

My MATLAB code in order to construct $m(t)$

t0= 0.15;
Fs= 2000;
t1= 0:1/Fs:t0/3- 1/Fs;
t2= t0/3:1/Fs:2*t0/3 - 1/Fs;
t3= 2*t0/3:1/Fs:3*t0/3-1/Fs;
m1= ones([1, length(t1)]);
m2= -2*ones([1, length(t2)]);
m3= zeros([1, length(t3)]);
m= [m1, m2, m3];
t= [t1, t2, t3];

I'm trying to find frequency modulation and the following formula is wanted to be used
$$
x_{FM}(t) = \cos\left(2 \pi  f_c t + 2\pi  k_f \int_{-\infty}^t m(x) dx\right)
$$
with carrier frequency $f_c=200$Hz and deviation constant $k_f= 50$.

Only for or while loops are allowed while we find integral. In order to integrate, I have made this for loop

result=0;
j=1;

for i= 0:1/Fs:t0- 1/Fs
   result= result + m(j);

    %% result2 will be used as we integrate it over 0 to 0.15
   result2(j) = result;
   j= j+1;
end

Please tell me, is something wrong with this integration?
I plotted this and it seems correct to me.

If the integration is true, then why frequency didn't change at all when plotting this signal.

xFM= cos(2*pi*fc*t + 2*pi*kf*result)

Full code here

Integration was incorrect. I was integrating it with $dx= 1$. Doing m(j)*1/Fsoperation in the for loop, I have made $dx= 1/Fs$ now.

As can be seen, frequency changes for different amplitudes now.

Moreover, resultis a constant, it should be changed toresult2 so that it gives result with respect to time $t$.


Why can't we have IoT from 2005? but begin now,or future? In 2005,ITU announce the "IoT",but why can't we have IoT at that time?

Is it because the speed of internet is not fast enough?wireless communication is not good at that time?Does anyone have some good reason which can explain that we just have the IoT concept in 2005,but we can not have the IoT in the real world at that time?

IoT requires inexpensive low-power communication. There were no solution with those 3 characteristics in 2005. Also, it was harder to monetize data so the perception was that IoT would never payback.


Correlation (or something similar) analysis between data with different sampling rates I tried to do correlation analysis on the data with different sampling rates.

For example, I have two sensors, A and B. Each sensor collects data with different sampling rates, e.g., A: 10Hz, B: 1Hz.

I turned on and off these two sensors at the same time, and collect data for 2 minutes.

Then, sensor A will give 1200 samples of data, and sensor B will give 120 samples.

I want to know the relation of these two data samples, like the values or trends from sensor A and sensor B are very similar.

Correlation may be one method to show such a relation, but how can I do correlation analysis in case of different sampling rates between two data?

You can adjust the sampling rate of a sequence. For example, the data from sensor B can be interpolated to the same sampling frequency as sensor A.

Sample rate increase is called interpolation and, if working in matlab, there is a function for this.

Lowering the sample rate is called decimation, also supported in matlab.

Once at the same sampling rate, correlation is easier to calculate.


Amplitude Modulation with large carrier and suppressed carrier From what I have read so far Amplitude modulation with large carrier is favaroble because envelope detector can be used at the receiver to derive the message signal and therefore another mixer is not needed at the receiver.

According to my text book one way to find the envelope of a signal is to lowpass filter the absolute value of the modulated signal.

Using low pass filter on the frequency spectrum $\mathcal F\left\{\lvert (m(t)+1)\cos(\omega t)\rvert\right\}$ will give $m(t)+1$ which is the envelope offset by constant unity.

My question is why can't we use the same technique for AM with suppressed carrier modulation?
Wouldn't low pass filtering of the spectrum of $\mathcal F\left\{\lvert m(t)\cos(\omega t)\rvert\right\}$ give us $m(t)$?

In short: you can use envelope detection on a suppressed-carrier signal, but you'll recover $|m(t)|$, not $m(t)$.

One easy way to see it is by drawing the signals. Draw an arbitrary signal $m(t)$, which takes positive and negative values. Now draw the signal $s(t)=\cos(2\pi f_ct)\cdot m(t)$. To do this, draw $-m(t)$ on top of $m(t)$, and draw a cosine wave that is "enveloped" by $m(t)$ and $-m(t)$.

To find the envelope of $s(t)$, take its positive part only and draw lines connecting the peaks of the carrier. You'll see that what you obtain is $|m(t)|$.

By the way, the same is true for large carrier modulation; it just turns out that $1+m(t)$ is always positive, so $|1+m(t)|=1+m(t)$.


How do the magnitude and phase spectrum of an imaginary function look like? Say I have the function
$$x(t)=j \operatorname{rect}(t)$$
Is the phase spectrum even or odd?
I am confused whether the phase spectrum is an odd/even function of $\omega$ (angular frequency, Fourier transform variable).

This might be helpful: https://users.aber.ac.uk/ruw/teach/340/ft_symmetry.php
Your
example is straight forward. You are just multiplying with a constant factor $j$ and that factor carries through the Fourier Transform, i.e. it's simply
$$\mathcal F\big\{j \cdot x(t)\big\} = j \cdot \mathcal F\big\{x(t)\big\}$$


What is the significance of Z-transform? As we have in Laplace transform that the roots decide the stability of the system i.e. if the roots are complex and lie in the left side of the plane you get a sinusoidal response with decreasing amplitude

similarly is there any significance of the roots , zeros and  ROC of the z-transform and the stability criteria . All i read in books is how to find the ROC and the properties of z transform like linearity ,time reversal,time shifting . It has nowhere mentioned why we are even using z-transform

my  apologies if this question is too basic or if it shouldn't belong here.

First of all, I think you're reading the wrong books. Almost any basic text on DSP has a chapter on the $\mathcal{Z}$-transform and its significance to describe linear time-invariant (LTI) discrete-time systems. If you're looking for good (and free) books, take a look at this answer.

I will not repeat all the details you can find in those books (and in many other places), but let me just point out a few very basic things to get you started. Each (single) pole $p$ of the transfer function $H(z)$ of a causal LTI discrete-time system contributes a term

$$c\cdot p^nu[n]\tag{1}$$

to the system's impulse response, where $c$ is some constant, $p$ is the (possibly complex) pole, and $u[n]$ is the discrete-time unit step function. From (1) it is clear that this contribution only decays with time if $|p|<1$. So for a causal system to be stable we require that all the poles of the transfer function are inside the unit circle of the complex plane, i.e. they have magnitudes smaller than $1$. So if you're looking for analogies with the Laplace transform, the inside of the unit circle corresponds to the left half plane of the complex variable $s$. Furthermore, the unit circle of the $z$-plane corresponds to the $j\omega$-axis. Knowing these two things, it becomes very easy to carry over everything you know about transfer functions of continuous-time systems (Laplace transform) to the discrete-time domain ($\mathcal{Z}$-transform).


Band limited discrete signals division If I take two analog, band limited signals $ x(t) $ and $ y(t) $ with non-zero $y(t)$ I can define the division:
$$
 z(t) = \frac{x(t)}{y(t)}
$$
In general $z(t)$ may have a very wide bandwidth, also if the two signals had a small bandwidth.
I ideally may lowpass filter out a limited band part of $z(t)$ and make a discrete version $ z[n] $ without aliasing.
If I instead make division from the discrete versions $ x[n] $ and $ y[n] $ (they are band limited so I can also have a sampled discrete version that completely describe them) the result can be very different from the other discrete version $ z[n] $ after the lowpass antialias filter, in fact I have aliasing.
Is there a way to estimate in digital domain the $ z[n] $ that I would obtain making the division in analog domain and filtering out the lowpass part?

I've found an estimator that can work for a lot of similar situations, the idea is that you can rewrite the division:
$$ x(t) = z(t) \cdot y(t) \tag1$$
The problem now is to estimate $ z(t) $ that solves the equation.
One can start to find a constant coefficient that minimize the 2-norm. This coefficient will be the linear regression slope without intercept term ($\beta$), hence you can write an approximation of $ z(t) $ that minimize the 2-norm:
$$ x(t) \approx \beta \cdot y(t)\tag2 $$
This estimator is constant, hence it is band limited, and you can estimate the same $ \beta $ from the sampled data, in fact be valid that
$$ x[n] \approx \beta \cdot y[n]\tag3 $$
Obviously this estimation is quite useless but you can extend this approach in a kind of FIR way , first.
You could make a rectangular window and estimate the linear regression slope on $N$ elements.
$$ x[n] \approx \beta(x[(n-N),..,N], y[(n-N),..,n]) \cdot y[n] \tag4$$
I don't like a lot this estimator but should limit aliasing if $ N $ is  large enough.
I instead prefer an IIR like approach, using an exponential window. You can do this using a 0-order RLS algorithm. You can choose a bandwidth choosing the forgetting factor $ \lambda = e^{-\frac{1}{W_s}}$,
where $ W_s $ is the width in samples of the exponential window, that can be expressed in function of cut off pulsation ($\omega_c$) using the relation
$$ W_s = \frac{f_s}{\omega_c} \tag 5$$
where $f_s$ is the sampling frequency.
With RLS you can estimate the slope $\beta$ in a time-variant way, with a limited band in discrete domain, where the bandwidth is limited by $ \omega_c $. Note that in this case (0-order) RLS becomes a scalar problem so it is very easy to implement.


Get all local maxima of a DFT I'm doing a Machine Learning project that incorporates many DSP elements in feature extraction.  I am looking at the STFT of a piece of audio and trying to search for specific spectral features within it - for example, the exponential decay of the fundamental frequency of a kick drum - to determine what kind of sound the piece of audio is.

One thing I've looked into is finding the frequencies of all the local maxima of a DFT, unquantized.  This includes being able to identify when the peak is "between" bins based on contextual information.  A possible path to go down is to find all of the changes in direction within the spectrum, and then for N changes fit an $(N+1)$-degree polynomial using a regression, then finding all the zeroes of the derivative of this polynomial.

But does an algorithm for this already exist?  It seems like something that might have been done before.

That sounds a lot like you shouldn't be using the DFT, but rather a parametric spectrum estimator.

For example, there's spectrum estimators that actually give you a function to evaluate at real-valued points rather than a vector of a sampled spectrum (your sentence "DFT, unquantized" doesn't make all that sense to me – the DFT is discrete, inherently, so what you describe as search algorithm only works if you do a massive DFT and then start somewhere in the and look for peaks from these points on).

My go-to example here is either

ESPRIT, which, given knowledge about the number of significant signals will just give you a vector of real-valued frequencies of peaks or
ROOT-MUSIC (dunno if that's the official term, it's what we used, based on [1]) which uses the pseudospectrum function generated by the MUSIC algorithm and finds the signals by looking for roots of a specifically crafted polynomial [1].

[1]: Barabell, Arthur J.: Improving the resolution performance of eigenstructure-based direction-finding algorithms. In: Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP ’83., Band 8, April 1983, available online.


How do you interpret FFT of an FFT of a discrete signal? I am new to DSP & Radar domain. I have one question.
Say I have a discrete time signal $x[n]$ and I take an FFT of it say the output is $X[k]$. The result what i get is frequencies present in $x[n]$.

Now what will happen if I take fft on $X[k]$ ? How would the spectrum look and what would you interpret from it?

thanks in advance !

Essentially, what you do is performing the DFT twice. Performing DFT twice amounts to time-reversing the signal:

N = 16; sN= np.sqrt(N)
x = np.random.randn(N)

X = np.fft.fft(x) / sN

x2 = np.fft.fft(X) / sN

plt.plot(x.real)
plt.plot(x2.real)

blue is original signal, green is two-times DFT. As you can see, the signals are time-reversed in the sense of $y[n]=x[-n \mod N]=x[N-n], n=0...N-1$. (y is green, x is blue). Here, I used the unitary version of the DFT, to keep the signals in the same scale. Otherwise, with the "normal" definition of DFT, the green signal would be scaled by 16.


Design and a notch filter using two Butterworth filters in parallel (MATLAB) I am doing a project to filter out a 1kHz noise inserted to a 30sec .wav audio file.
To filter out the 1kHz noise, I am required to design using two Butterworth filters connected in parallel where one is a low pass and the other is a high pass Butterworth Filter.

I have tried to notch out the 1kHz noise by setting the bandwidth to (+/-10Hz) of the 1kHz noise.
Therefore, the range of the noise signal to be attenuated will be from 990Hz to 1010Hz.

I am planning to obtain the transfer function in s-domain first and then transform it to z-domain.

The designed transfer function will be inputted into MATLAB and the filtered audio file will be obtained.

May I know the steps in designing this filter? Cost-effectiveness and ease of implementation of the filter is emphasised.

Thanks.

Cost-effectiveness and ease of implementation of the filter is emphasised.

It that's a your goal then you really shouldn't be using two Butterworth filters. Use a notch filter instead. If there is a "requirement" to use other filters, please explain why this is required. I think your best bet here is to change the requirement.


How to extract pixel-specific time-stream from mp4 video I am interested in extracting specific timechange information from an mp4 1900x1200 20fps video.The only pixels I want to process can be found in its borders. Therefore I only need to process approximately 6200 pixels/frame. I am using Matlab and VideoReader function, which is terribly slow. Is there a way to extract only one pixel from each frame of the video throughout its duration, in order to speed up the process? Or am I obliged to load the whole frame? To rephrase the question, if the video is considered to be a mxnxt complex array and the frames are mxn size, can I get a 1x1xt stream instead of mxnx1? I tried to bruteforce by finding information considering the way mp4 is encoded, with no success. If I changed the video encoding or used any other program would that be possible?

In video, pixels are not usually encoded independently; you will need to decode full frames. If you are going to run several experiments on a single video, I would suggest that you decode the video once frame-by-frame and write the border pixels into a file in your own raw format for easy access. This should give a manageable file size perhaps of a few gigabytes for a movie-length video.


How to get general projection of rect function? To be more clear I want to get a general expression for

PƟ(t)= ∫∫rect(x)rect(y)δ(xcosƟ+ysinƟ-t)dxdy

I also want to find particular projection for Ɵ=0 and Ɵ=45 and also Fourier Transform of the projections at 0 and 45. I kinda new at image processing and projections, thus I could not figure out a way.

Thanks in advance.

Since your expression
$$\begin{align*}P(\theta,t) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos\theta + y\sin\theta -t\right)dxdy\\
\end{align*}$$

uses a delta function with only 1 dimensional argument, in the 2 dimensional integration I assume you mean an infinite impulse line (an impulse "wall") and not a point impulse (an impulse "spike").

Let's just look at the general solution first, and then handle the degenrate case of $\theta = 0$.

I'll be using the following properties of the delta function:

$$\begin{align*}\delta(ax+b) &= \dfrac{1}{|a|}\delta\left(x + \dfrac{b}{a}\right)\\
\\
\int_{-\infty}^{\infty} f(x)\delta(x-a)dx &= f(a)\\
\end{align*}$$

So, I believe the solution is this:

$$\begin{align*}P(\theta,t) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos\theta + y\sin\theta -t\right)dxdy\\
\\
&= \int_{-\infty}^{\infty} \Pi(y) \left[\int_{-\infty}^{\infty} \Pi(x)\dfrac{1}{|\cos\theta|}\delta\left(x - \left[-\dfrac{\sin\theta}{\cos\theta}y +\dfrac{t}{\cos\theta}\right]\right)dx\right]dy\\
\\
&= \dfrac{1}{|\cos\theta|}\int_{-\infty}^{\infty} \Pi(y) \Pi\left(-\dfrac{\sin\theta}{\cos\theta}y +\dfrac{t}{\cos\theta}\right)dy\\
\\
&\mbox{and if the rectangles have any overlap for the given $t$ and $\theta$}\\
\\
&= \dfrac{1}{|\cos\theta|}\int_{\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}^{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)} dy\\
\\
&= \dfrac{y}{|\cos\theta|}\biggr\rvert_{\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}^{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}\\
\\
P(\theta,t) &= \dfrac{\min\left(\frac{1}{2},\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)-\max\left(-\frac{1}{2},-\frac{1}{2}\frac{\cos\theta}{\sin\theta}+\frac{t}{\sin\theta}\right)}{|\cos\theta|}
\end{align*}$$

So for the degenerate case of $\theta = 0$, the above expression runs into problems, but the integral for this special case is a bit easier:

$$\begin{align*}P(0,t) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Pi(y)\Pi(x)\delta\left(x\cos{0} + y\sin{0} -t\right)dxdy\\
\\
&= \int_{-\infty}^{\infty} \Pi(y) \left[\int_{-\infty}^{\infty} \Pi(x)\delta\left(x-t\right)dx \right]dy\\
\\
&= \Pi(t)\int_{-\infty}^{\infty} \Pi(y) dy\\
\\
&= \Pi(t)
\end{align*}$$

The one dimensional Fourier Transform of $P(0,t)$ is a simple table look-up to find that it is a $\mathrm{sinc}()$ function.

The one dimensional Fourier Transform of the more general $P(\theta, t)$ is obviously more complicated.  There are 4 cases there.


What data should be collected to determine the autocorrelation sequence? For a current project I need to implement an algorithm that needs the autocorrelation sequence of the channel, as part of a key step.

We're measuring the frequency utilisation in our local area - simply be scanning the spectrum (from 440 MHz to 780MHz) with an SDR and collecting the power at each frequency. The goal of the algorithm is to determine how the spectrum may be split into bands to transmit on.

I want to ask, is this data sufficient for what we need to do or is there anything else we would need to collect?

All you need to do an autocorrelation is a number of evenly spaced samples at a given frequency, so you will need some dwell time at each frequency (which I assume that you are doing already).  That is all you need.


Estimation of Amplitude, Frequency and Phase of Linear Combination of Harmonic Signal Beyond the Leakage Resolution of DFT How can I find a rough ( as accurate as possible) Amplitude of each frequency when there is spectral leakage.
Currently, I am dealing with a system that contains special leakage which seems unavoidable as I am measuring a real signal with any possible frequency. currently, I can identify the component frequencies to a decent level of accuracy even with the leakage, but I can't find the amplitude of this frequency as the energy is spread over multiple frequency bins.
Is it possible to get a rough prediction of a frequency amplitude either directly from the FFT, or afterwords using the identified frequencies?
I have been experimenting with Windowing and found this useful.
https://www.youtube.com/watch?v=VxTx9QW8Zx8&ab_channel=Adash
using the Hann window and using a frequency correction equation. it gives good approximations for frequency but no amplitude or phase.
A sample of the code currently used to test and generate the spectra.
import numpy as np

def gendata(FreqR,AmplR ,Ssize):
    PhaseR = np.pi # phase range
    FreqR = (FreqR*10) +1
    AmplR = (AmplR*100)+1
    PhaseR = (PhaseR*100)+1

    freq_OR = []
    ampl_OR= []
    phase_OR =[]

    for i in range(Ssize):
        freq_OR. append((((FreqR-1)/Ssize)*(i+1)/10)) # currently just equaly spread but can be any value between 0-9 rad/s
        ampl_OR.append(1) # unit amplitude to check for variation and error.
        phase_OR.append(np.random.randint((-PhaseR),(PhaseR))/100) # random phase

    return(freq_OR,ampl_OR,phase_OR)

def pltsignal(freq_OR,ampl_OR,phase_OR,t_OR):

    z = 0 # original signal flat sea's
    if len(freq_OR) == len(ampl_OR):
        for i in range(len(freq_OR)): # for each frequency generate a regular wave.
            wave_OR = ampl_OR[i] * np.cos(2*np.pi*freq_OR[i]*t_OR + phase_OR[i]) # generated wave signal using the frequncy, amplitude and phase of each wave.

            z = z + wave_OR # superpostion each wave ontop of each other.
    else:
        print("amplutide and frequency are diffrent lengths")

    return z  # returns superimposed irregular wave.

FreqR = 1.5 # htz or (9 rad/s ish)
AmplR = 1 # currently unit amplitude to measure the error
Size = 12 # sample size, also used to move around frequencies

E_OR = 20 # time over wich signal is measured (max values is about 120 ish)
Fs = 6 # Sample rate in Hz

t_OR = np.arange(0,E_OR,1/Fs) # starts at 0, ends at E, steps by 1/Fs

data_OR = gendata(FreqR,AmplR,Size)
z = pltsignal(data_OR[0],data_OR[1],data_OR[2],t_OR) # original signal ( simulation of mesurment )

# anaysisng the signal using NumPy fft, then scaling and finding peaks.
Ramp = np.fft.fft(z) #real aplitudes used to find phase
Rfeq = np.fft.fftfreq(z.shape[-1]) #real frequency domain to find phase

Signal Model:
$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \cos \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$
Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).
One could assume:

The SNR is very very high.
The observation time is ~ 120 [Sec].
The number of signals is 1-20.
The frequencies are up to 2 [Hz].
The gap between frequencies can be as small as 0.005 [Hz].

Solving the Linear Combination of Real Harmonic Signal
The data model is given by:
$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \sin \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$
Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).

Estimating the parameters above requires 4 steps solution:

Estimate the model order ($ M $) / number of signals.
Estimate the frequencies accurately using direct method (Given the number of frequencies).
Estimate the amplitude and phase of the model given the frequencies (Usually using Non Linear Least Squares solver).
Estimate all amplitudes, frequencies and phase (Usually using Non Linear Least Squares solver).

The above works only in very very high SNR cases. For low SNR this problem is probably unsolvable (Certainly for large $ M $).
If the frequencies are well spread (Far from each other) then steps (1) and (2) can be achieved using DFT (fft()).
Yet in case the frequencies are close in the best case the estimation of the frequency won't be accurate (Leakage from near by frequencies will shift the peak) and in the worse the case frequencies will be blended and neither the number of signals nor the frequencies is achievable.

In the images above I used observation time of 12 [Sec] with 4 signals with frequencies of [0.95; 1.00; 1.05; 1.10] [Hz]. Since I took 1 / 10 of observation time in your question, it is equivalent of frequencies being 0.005 far apart in the 120 [Sec] observation time. I just wanted to deal with less samples.
As can be seen in the image above, in the DFT there is no way estimating the number of signals and certainly extracting the frequencies.
But using Super Resolution methods one could solve this case!
See some information about Frequency Domain Super Resolution in the section below.
Yet for this case I used different method for Super Resolution which I can't reveal (Used commercially by me as an advisor).
I also developed a method to estimate the number of signals.
In the case above indeed my method estimated 4 signals and the estimated frequencies were [0.9503, 1.0016, 1.0483, 1.0996] [Hz].
Then applying the next 2 steps I could even farther improve results:

As can be seen from above, the estimated signal is almost perfectly aligned with the ground truth (The model with known parameters with no noise).
Indeed the SNR here is very high yet still, without the Super Resolution method, it wouldn't work.
The full code is available on my StackExchange Signal Processing Q74024 GitHub Repository.
Unfortunately, I can't disclose the code of the estimation of the number of parameters and the Super Resolution.
Frequency Domain Super Resolution
Old non modern methods will use some Windowing tricks which falls short in most cases.
What you really need is Super Resolution. Super Resolution basically is equivalent to a "Window Method" with very narrow main lobe and almost no side lobes.
You can employ Compressed Sensing / Sparse Representation for Super Resolution in Frequency Domain.
One way to do so is solving the problem:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| F \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$
Super Resolution means, in that context, being able to resolve frequencies which are closer than what the observation time suggests (overcoming Leakage issues):

In the above you can see the DFT of a sum of 2 sines with the given frequencies. The Gaussian model is using $ {L}_{2} $ for regularization (Which is basically damped zero padding).
You may see that the $ {L}_{1} $ could resolve the 2 sines even when they are only 0.5 [Hz] apart with an observation windows of 1 [Sec].
The answer was taken from my answer for Super Resolution in Frequency Domain Using Compressed Sensing.


How to characterize the regularity of rhythm? Again, my snoring detector:

When I look at the rhythm of snoring vs the rhythm of talking in terms of peaks/events (which would represent snores) per hour, there's not much difference -- they range from roughly 500 to 900 per hour.  But it occurred to me that perhaps there's a significant difference in the regularity of the rhythm.  In particular, I suspect that talking would have a more irregular rhythm than snoring, and I was wondering what might be some ways to characterize that difference.

Note that this is a bit different from tempo/beat detection in music in that there is no master clock, and individual "beats" are not "skipped" but elongated, so, eg, an autocorrelation would likely not reveal much.

The simplest approach would be to measure standard deviation of the event intervals, but I'm unsure to what extent that would discriminate between apneic breathing (snore, snore, snore, pause, snore) vs simply irregular intervals.  And, as you know, my statistics background is weak.

Any other suggestions?  (Take as a given that I have already discriminated the "events" and can produce a train of intervals, vs deriving rhythm from raw data.)

You have already found train of intervals between the events. What you are interested in is a respresentation of that in time. I would suggest, from my experience, to run a sliding window differential filter, with some low pass averaging filter at the output, that will give you the rate of changes in time. In this way you will get the rate of changes as they progress in time, which is most likely what you actually need. Yes, standard deviation may be done, but that is just one number, as you said, and that is not sufficient.
One simpler method will be to use just a sliding average window that averages train of intervals and gives that as a result.
With this, youy may, at the end, postprocess obtained result, by finding max. min. and average values, as the informative, yet simple measurements of the signal that has been analyzed.
If needed, you may get more complex solutions, involving cyclostationary data analysis. In a nutshell, finding some statistics functions, (sorry it cannot be avoided if you decide to go this way): for example mean, variance, cyclic correlation function, cyclic power spectral density etc. This will lead to a rate detector, with various outputs that characterize the signal rate changes (and much much more). I hope that this helps.


Baseband vs. Passband I am a bit confused if the following signal is baseband or passband:

$f(t) = x(t) \cos(2 \pi f_c t) - j \: y(t) \sin(2 \pi f_c t),$
where $j = \sqrt{-1}$.

Since it is clearly not real, I cannot say it is passband. On the other hand, since it contains the frequency $f_c$, it should not be a baseband. Am I missing something?

Thanks,

Both real and complex signals can be passband.

Assuming $x(t)$ and $y(t)$ are baseband signals with bandwidth much less than $f_c$, then $f(t)$ is a passband signal, even if it is complex.


How does filter bank work in image processing? I read a paper about directional filter bank in image processing and I don't understand clearly about the process. It says they use '2-band structure' like this image:

could someone here please explain it to me how this process works? I would like to implement it on python or MATLAB.

Sorry, I couldn't access the paper behing the paywall, so here is an educated guess about the figure that you linked in the question.

$H_0$ is a low-pass filter, while $H_1$ is a high-pass filter. You can note that the union of both domains correctly yield the full bandwidth of the input signal.
These filters are directional, i.e., they operate on different subbands that correspond to a set of orientations. This is shown in the figure by the 1,2,3,4 triangles: each triangle corresponds to a set of low/high-frequency of a given set of orientations. For example, 1 corrresponds to lower frequencies, and the slope of the lines in 1 is positive.
Finally, the directional filter bank consists in two operations: separatig low/high frequencies (the filter bank part) and applying a transfer function that also depends on the orientation of the current frequency (the directional part).

I'm note sure that you can implement the figure from this explanation though. A good starting point now would be for you to look for an example code that does similar things.


Histogram of signal with a time shift added One of the problems in Proakis's signal processing book using Matlab asks to plot the histogram of

x3(n) = x1(n) + x1(n-1)

where x1(n) is a random sequence.

To create x1(n), I do the following:

n1 = 1:100000;
a = 0; b = 2; % interval over which to pick random numbers.
x1 = a + (b-a)*rand(100000,1);
N1 = hist(x1,100);
figure, bar(N1);

To generate x3(n), I use a couple of helper-functions defined in Proakis's book.

[x31,n31] = sigshift(x1,n1,1);
[x3,n3] = sigadd(x1,n1,x31,n31);
N3 = hist(x3,nbins);
figure, bar(N3);

Unfortunately, I do not have enough points to post images. However, the histogram of the random sequence x1(n) just shows that x1 is uniformly distributed. The x3(n) however, looks more like a Gaussian distribution.

I don't understand why this happens, and would appreciate if someone could explain this a bit.

Thanks.

In your example, $x_3[n]$ is the sum of two iid random variables, $x_1[n]$ and $x_1[n-1]$. Since each is taken from the uniform distribution, the probability density function (pdf) of each is:

$$
f_{x_1}(x) = \begin{cases}
\frac{1}{b-a}, & x \in [a,b] \\
0, & \text{otherwise}
\end{cases}
$$

This is just a rectangle function that is scaled and shifted to lie on the interval $[a,b]$. One property of the sum of independent random variables is that the pdf of the sum is equal to the convolution of the individual random variables' pdfs. Therefore:

$$
f_{x_3}(x) = f_{x_1}(x) * f_{x_1}(x)
$$

If you convolve two rectangles together, you get a triangle function. Therefore, if you histogram enough samples of $x_3[n]$, you should see the plot converge to a triangular shape.


Definition of DSP Terms can someone define/explain briefly (or with some detail-> up to you) the following terms:

Block Processing? and what are the differences between it and Sample Processing?
(removed)
What is a 3D Audio System and what is relation between it & Head Related Transfer Function (HRTF)? and how can HRTF be applied in Headphones?
(removed)
The relationship between bit rate, sample rate and bit depth? Is there a formula that combines all terms?

Block Processing (sometimes called Vector Processing) is a programming technique similar to, but not exactly the same as First-In-First-Out (FIFO) buffering.  It is double buffering of the samples coming into the process (be they from an A/D converter or a stream) and going out (to a D/A converter or another stream).  Block processing divides the processing time into slots of the entire block.  I.e. in audio, 32-sample blocks might be common in 48 kHz samplerate data.  So in one block time (which is 32/48000 second), the DSP algorithm is presented with 32 valid input samples at the beginning of the time period, and has the entire block time (of 32 sample periods) to compute 32 valid output samples.  This allows the cost of some of the work, like loading/saving states and loading coefficients and initializing accumulators and such, to be amortized over all 32 samples.  This makes the processing more efficient.  But the cost is 2 blocks of delay, due to the double-buffering.  And if there is feedback from one processing algorithm to another, that block of samples being fed back is from the previous block of samples, not the current block.  So feedback has a delay of 1 block of samples.
—
HRTF, I'll have to come back to that.  Maybe someone else can pick it up.
—
Sample rate is a well-known definition.  It is the number of samples per unit time expressed as frequency.  Usually it's samples-per-second expressed as Hz.  Bit depth is the same as the word-width of the samples.  CDs have 16-bit depth.  Decent audio processing DAWs have 24 or 32-bit depth.  Bit rate for uncompressed audio is simply the sample rate times the bit depth (times the number of channels).  It's the data rate of the audio stream.  For compressed audio, the bit rate is much reduces and continues to be the data rate of the audio stream.


Fourier transform with periodicity at the harmonic frequency Let's suppose I have a signal $F(t)$ that is periodic, with two periodicities $P1$ and $P2$, with $P1>P2$

Suppose that I know the values of the two periodicities. Using the Fast Fourier transform I can show the two values as peaks in a power spectrum.

Now, let's suppose the second periodicity $P2$ (the faster one), has exactly the same value as the first harmonic of the fundamental value, or $P2=2×P1$. This means that I will be not able to distinguish it by using the power spectrum, at least not by looking at the frequency of the peak.

My question is: is there a way to separate the contributions in such a case?

For example, is it possible to predict the power of the first harmonic, so that the difference between the predicted power and the observed power of the harmonic peak gives a result significant enough (i.e., greater than $3\sigma$) to say that the first harmonic also "contains" the contribution from a periodicity?

Please, be plain I am not experienced in this (nonetheless some equations/numbers are ok).

If you have two periodic functions $f(t)$ and $g(t)$, where $f(t)$ has fundamental frequency $f_0$ and $g(t)$ has fundamental frequency $2f_0$, then their sum $h(t)=f(t)+g(t)$ is periodic with fundamental frequency $f_0$. Without any further knowledge about $f(t)$ or $g(t)$ there is no way to separate the two, because if the only knowledge you have about $f(t)$ is that it is periodic with fundamental frequency $f_0$, then there is no way to distinguish it from the sum function $h(t)$, which has the same periodicity.


Speech SNR estimation I'm trying to estimate SNR level in db to use it as a switch for an speech enhancement method that i'm working on. The noise level is not known and can be estimated using a posterior SNR algorithm. Considering having noisy speech signal power spectrum and estimated noise power spectrum (using a posterior SNR estimation method) how can one estimate signal SNR level in db? The reason for estimating the SNR in db is that if the SNR level rises above a known value e.g. 10db the enhancement algorithm should be switched to another one. How can i estimate SNR in db accurately?

Here is the C++ code of the answer.

#include <algorithm>
#include <queue>
#include <numeric>
#include <vector>
#include <cmath>
const float N = 4;
float noise_th= 0.35;    //noise threshold
float calculate_SNR(queue<float> &signal_fifo)
{
    vector<float> noise_vect;
    vector<float> speech_vect;
    float SNR;
    float noise_std;
    float speech_std;
    auto fast_std = [](vector<float> input) //https://stackoverflow.com/questions/7616511/calculate-mean-and-standard-deviation-from-a-vector-of-samples-in-c-using-boos/12405793#12405793
    {
        float mean, acc;
        mean = accumulate(begin(input), end(input), 0.0)/input.size();
        acc = 0.0;
        for_each(begin(input), end(input), [&](const float d)
        {
            acc += (d - mean)*(d - mean);
        });
        return sqrt(acc/(input.size()-1));
    };

    while(!signal_fifo.empty()) //Seprate speech and noise samples(fifo has the ABS of input signal)
    {
        if(signal_fifo.front() < noise_th)
        {
            noise_vect.push_back(signal_fifo.front());
            signal_fifo.pop();
        }
        else
        {
            speech_vect.push_back(signal_fifo.front());
            signal_fifo.pop();
        }
    }

    noise_std = fast_std(noise_vect);
    speech_std = fast_std(speech_vect);
    SNR = 20*log10(speech_std/noise_std);
    noise_th = N*noise_std;     //Update threshold value
    return SNR;
}


Variance of deterministic sine wave I have read that variance of a deterministic signal is 0 because it is deterministic and does not vary over time.

In this explanation, I am considering a $\sin(t)$ wave where value at each t is a variable(it will be a random variable when we are doing measurements but I am not sure what to call it here, so just variables). We can have sampled values of $\sin(t)$ for limited number of variables but it is not related to my question.

I understand very well for a single variable that it's variance is zero. For eg: value of sine wave at $t=1$ which is given by $x=\sin(1)$ . Since $\sin(1)$ is constant every time you try to measure it, variance given by the following formula is zero as well:
$$ var(x)= E[(x-\bar x)(x-\bar x)^T] $$
I understand this is zero because $x=\bar x$ every time you calculate x.
Now, suppose you have a deterministic signal($\sin(x)$ from $0$ to $2\pi$). Now, mean of the signal is 0 and variance is just the squared sum of signal values from $0$ to $2\pi$ which is not zero. Even more, this is providing the energy of the signal.
Does that mean variance of a deterministic signal is defined and does not have to be zero. Is variance even defined for a whole signal or is this concept valid only for a single random variable.

If $x(t)$ is deterministic its variance is zero. Note that the mean of $x(t)$ (defined as an expectation) just simply equals $x(t)$:

$$\mu_x(t)=E\{x(t)\}=x(t)\tag{1}$$

where the last equality follows from the fact that $x(t)$ is deterministic. From $(1)$ we get the variance

$$E\{[x(t)-\mu_x(t)]^2\}=E\{[x(t)-x(t)]^2\}=0\tag{2}$$


Audio signal dither and noise shaping I'm trying to get a handle on the importance of the error feedback term in noise shaping operation in typical audio dither algorithms.

I'm thinking in terms of four signals.  The original signal, the output (quantised) signal, the dither signal, and the error signal (the difference between output and original).

A simple approach would be to generate white noise and apply a filter to it to generate a coloured dither signal, add this to the original signal, quantise to the output signal, and disregard the error signal.  It appears to be a fairly common approach to use white (unfiltered) or violet (differentiated) noise functions in this configuration.

Another approach I see implemented is to take the error signal and filter that and add white noise to define the dither signal, then add it to the original, quantise, and calculate the next error term for the filter.  This implies a tight loop for the data flow -- we cannot compute a new error term without resolving the previous one -- making it hard to use SIMD efficiently.

From those I can infer a third option -- forgo the noise altogether and simply filter and reintegrate the error signal directly.  I suppose this comes with some risk of modulation artefacts related to but distinct from the artefacts introduced by quantisation without any dither or shaping at all; however, I suspect this might be how it's done in some high-frequency 1-bit cases.

The point I'm unclear on is that complex filters which attempt to model the threshold of human hearing appear to come under the title of 'noise shaping', and this implies (supported by implementations I have seen) use of the error feedback term.

Is it a flawed or misguided approach to forgo the error feedback but still use the complex filter on white noise for the dither function?  Is the error signal genuinely offering something important here?

i dunno what a "Lipshitz filter" is but i do know who Stanley Lipshitz is.  now, i don't know for sure what Stanley would say to you, but the best advice i got from him is that it does little good to color the additive dither if you are also using error feedback.  if there is no noise shaping feedback, probably the simplest way to get triangular PDF dither is to generate good, independent, rectangular PDF random numbers $v[n]$ and run that through a differentiator.

$$ d[n] = v[n] - v[n-1] $$

the result is high-pass in the spectrum and triangular PDF.  if $\Delta$ is the step size of the quantization, $\frac{\Delta^2}{12}$ is the variance (or power) of raw quantizer, but you have two more of them due to your dither and $\frac{\Delta^2}{4}$ is the total power due to the entire dithered quantizer.  a 4.77 dB increase in quantization error, but the quantization is fully decoupled from both the mean and variance of the input signal.  although this is the same quantization power that you would get from adding two independent random numbers, i think you would prefer that much of the power is steered to higher frequencies (but not all of it, the spectrum of the raw quantizer is still flat).

if you're planning on error feedback, i think the best advice is either use white, triangular PDF dither (requires adding two independently generated rectangular PDF random numbers) and error feedback (designing a good feedback filter takes a bit).

or no dither (if the quantization step size is very small) and error feedback.  the simplest no-trivial example of this is "fraction saving" where you simply hack off the low-order bits (the floor() function) and add those bits back in (zero-extended on the left) to the next sample before that sample is quantized (with the floor() function).  it has infinite S/N ratio at DC and totally kills any limit cycle in feedback filters (like biquads) that get the filter stuck on a DC value even when the input goes to zero.  i've done that a lot.


how does DSD silence looks like? DSD is 1bit PWM Sigma Delta audio format.In PCM that dont have dithering,silence looks like samples with 0 value,I wonder,what does DSD silence looks like.Is it just series of 0 value samples,or does it switch between 0 and 1 constantly,third option,it changes between 0 and 1,but not constantly,more in noisy random way.

Many DACs these days use are multibit sigma delta,how does multibit sigma delta silence looks like?

DSD bits are 0 and 1 indeed, but they actually represent normalized outputs of -1 and 1, respectively, since one-bit samples can only be mid-raiser-quantized. This means DSD samples cannot represent zero. One can generate DSD silence by alternating 0-s and 1-s (your second option), or by a 1-bit delta-sigma modulation of a zero PCM input, with or without dithering (your third option).
Multi-bit modulators can (and normally use) mid-tread quantization, so they can represent zero. Thus, multibit silence is either a sequence of zeroes or the output of a multi-bit delta-sigma modulator presented with a zero input.


How to deal with a negative pole (unstable) in the pre-filter of a control system? So while answering how to design a PI controller for a first order time delayed system   (Question Here )

Here is the closed loop equation to a control system:

$$
G_C(s) = \frac{\frac{K}{T}(1-sT)(s)} { s^3 + (\frac{1}{T} + a - KK_p)s^2 + (\frac{a}{T} + \frac{KK_P}{T} +K_I)s+\frac{KK_I}{T}}
$$

Question: How do you deal with normalizing the numerator in your closed loop transfer function when the filter is unstable? (Pole on RH of plane)

Typically you introduce a filter before your controller that does:

$$ \frac{1} {\frac{K}{T} (1-sT)(s)} $$

to normalize the numerator

But the filter itself is unstable because of the term:

$$ \frac{1}{(1-sT)}$$ is unstable for a step response which would create an issue realizing the system at all.

One way I've thought about dealing with this is multiplying it by its complex conjugate
$$
\frac{(1+sT)} {(1+sT)}$$

but im not really sure about the merits of it.

In general, in order to stabilise a system, however complex, if you have the transfer function $G(s)$, you introduce a feedback loop with a new function $F(s)$.

Write the closed loop transfer function for the new system with the added $F(s)$, and then find $F(s)$ in order for the new system to be stable. This is like the first exercise in any control book in order to offer an example of stabilising a system via negative Feedback.

Check the book from Ogata on control engineering for reference.


Efficient computation of inverse FFT of real sequence using half sized FFT or inverse FFT The inverse FFT of a folded real sequence yields a real symmetric sequence. Are there any ways to exploit this symmetry to perform a smaller complex FFT or complex inverse FFT?

Background: it's possible to do this for a forward FFT of a real sequence (see: Real valued FFT in Matlab):

x = [1, 2, 3, 4, 5, 6, 7, 8];   % some real-valued signal
n = length(x);                  % must be even!
n2 = n/2;                       % assume n is even
z = x(1:2:n)+j*x(2:2:n);        % complex signal of length n/2
Z = fft(z);
Ze = .5*( Z + conj([Z(1),Z(n2:-1:2)]) );        % even part
Zo = -.5*j*( Z - conj([Z(1),Z(n2:-1:2)]) );     % odd part
X = [Ze,Ze(1)] + exp(-j*2*pi/n*(0:n2)).*[Zo,Zo(1)];  % combine
X2 = fft(x);

Then X and X2 yield the same values. With X we have avoided the calculation of the symmetric values as we can recover those easily.

For the ifft:

x = [1, 2, 3, 4, 3, 2];   % some folded/symmetric real-valued signal (of even length)

ifft(x)

This yields a real symmetric sequence.

Would it be possible to take [1, 2, 3, 4] and pack this into a complex FFT or inverse FFT to yield the same result?

As far as I understood your problem right, then the following code (modified from your code) provides an efficient computation of the $N$ point inverse-DFT of a real valued and even symmetric frequency sequence $X[k]$ in terms of an $M = N/2$ point inverse-DFT of an associated frequency sequence $Z[k]$, as the following MATLAB code shows:

It exploits the duality of operations between time to frequency and frequency to time mapping of DFT.

N = 10;                  % length N of frequency sequence X[k]
M = N/2;                 % M = N/2 , half-size
X = [4 3 2 1 0 -1 0 1 2 3];  % some real-valued even sequence X[k]
x = ifft(X,N)            % N-point I-DFT x[n] of N-point sequence X[k]

Z = X(1:2:N)+1j*X(2:2:N);% Half-Length complex frequency sequence Z[k]

z = ifft(Z,M);             % M-point DFT Z[k]

% Construct first M+1 points of real-even x[n] from Z[k]
ze = +0.25*( z + conj( [z(1) , z(M:-1:2)] ) );       % even part
zo = -0.25*1j*( z - conj([z(1),z(M:-1:2)]) );        % odd part
xz = [ze,ze(1)] + exp(1j*2*pi*(0:M)/N).*[zo,zo(1)]; % COMBINE

Is this really an efficient way? Eventhough it exploits the symmetry and avoids computing the half of the sequence, it makes quite large number of pre-processing and post-processing which should be taken into account...


Frequency Response Output I'll start off by saying, yes this is a homework problem, and no I do not want the answer. I simply need to know where to start.

The input to a system with a frequency response of $H(\omega)$ is $x(n)=A\sin(\omega_0 n+\theta)$

What is the expression for the systems output, $y(n)$.

As I said, I just am unsure of where to start.

Hint:

I guess we may assume that the system is real-valued, i.e. its impulse response $h(n)$ is real-valued. This means that the response to the imaginary part of some complex input signal $\hat{x}(n)$ equals the imaginary part of the response to the complex signal $\hat{x}(n)$:

$$\mathcal{T}\{\Im\{\hat{x}(n)\}\}=\Im\{\mathcal{T}\{\hat{x}(n)\}\}$$

where $\mathcal{T}$ is the system operator.

Note that your input signal $x(n)$ is the imaginary part of

$$\hat{x}(n)=Ae^{j(\omega_0n+\theta)}\tag{1}$$

The frequency response is related to $h(n)$ by

$$H(\omega)=\sum_{n=-\infty}^{\infty}h(n)e^{-jn\omega}\tag{2}$$

The output $y(n)$ is given by the following convolution sum

$$y(n)=\sum_{k=-\infty}^{\infty}h(k)x(n-k)\tag{3}$$

Now use (3) with the complex input signal $\hat{x}(n)$ given by (1), make use of (2), and then take the imaginary part to obtain $y(n)$. The response of a system with frequency response $H(\omega)$ to a complex exponential $e^{j(\omega_0n+\theta)}$ is a fundamental result which you should remember.


Why are there only a few high frequency components in the Discrete Cosine Transform? I am learning about the Discrete Cosine Transform (DCT) and I read somewhere that, in general, images have very few high frequency components. Can someone give me an intuitive explanation for why this is so? Is this assumption violated in any scenario?

There are two parts to the answer.  One has to do with the nature of sharp edges, the second has to do with the how imaging systems (lenses and the human eye) work.

As @PeterK's example demonstrates artificial (human created) images tend to have large regions of one color with sharp transitions between the regions.  This is similar to a square wave.  But if you look at the fourier transform of a square wave it will be a decreasing sum of the odd frequencies:
$$
\sum_{n=1, 3, 5, ...}\frac{1}{n} \sin{ n \omega t}
$$

Images of the natural world also often have regions of slowly changing color with sharp transitions where one object overlaps another.  There can actually be a lot of high frequency "texture" in an image (grass, leaves, a pebbly beach, waves on water that is relatively far away, whispy clouds) but the sharp edges of the occlusions (overlaps) obey the same square-wave fall-off in frequencies beyond the fundamental frequency (the average distance between tree trunks, blades of grass, etc.)

The second part of the answer is that imaging systems (camera lenses/eyes) can not make out features smaller than a certain minimum angular resolution.  Images outside of the imaging systems depth of field will be blurred because lenses have focal points.  But even objects inside the depth of field will be blurred.  This is because the aperture through which the light must pass must be reasonably large (several mm) in order to let in enough light.  But this means that point-sources of light actually create a circle of confusion that essentially acts as a low-pass filter.

Because of this the human eye, for example, is more sensitive at some frequencies than others.  In the below image the contrast differences across horizontal scans in the image are the same, but your eye is more sensitive to mid-frequency contrast differences than high-frequency contrast differences.  (You may have to adjust the distance from the image to your eye to see the effect line up with the black line drawn on the image.)


How is the frequency vector being derived? I am having a hard time understanding how the frequency vector is being derived in this code. I would appreciate any guidance from anyone.

 Frv = linspace(0, 1, fix(numel(FTSignal)/2)+1)*Fn;

linspace(0, 1, fix(numel(FTSignal)/2)+1) $\tag{1}$

"$y =$ linspace(x1,x2,n) generates n points. The spacing between the points is (x2-x1)/(n-1)."

So in the above statement it is generating fix(numel(FTSignal)/2) points between 0 and 1 that are equally spaced

  Frv = linspace(0, 1, fix(numel(FTSignal)/2)+1)*Fn;

This generates Frv with just equation (1) scaled/multiplied by Fn

Edited: the FFT signal can usually be represented by half of the total FFT coefficients, since FFT of a real signal is symmetric around mid point, it is uniquely represented by half the coefficients. Adding 1 is for the DC frequency. Take an example an FFT of length 16, is represented by 8+1 coefficients. The "fix" function is present to handle FFT of odd length signal, so that they can be rounded to an integer value. Since these FFT values have been generated under the normalized scale of 0 to 1, it is then normalized by Fn, which can be $\pi$ or $F_s/2$, where $F_s$ is sampling frequency


Why is the DFT used for spectrograms rather than the DCT? Related to Could a DCT be used for an audio magnitude spectrum rather than DFT?, but more on the (audio) spectrogram side of things

I've been reading about the DCT (such as from here) and it seemingly provides quite a few advantages over the DFT such as having more frequency bins for the same number of samples and better power compaction, which sounds useful for signal analysis.

Yet all the programs I've used that have spectrograms seems to use the (magnitude of) DFT rather than DCT. Why is that? Is there something bad or "missing" from the DCT that is important for spectrograms?

(Note: I already understand the benefits of the DCT for compression, and this question is about audio (or other time-displacement) signals, so please avoid focusing too much on compression or images)

Let's define the N-point IDFT $y[n]$ of a signal $Y[f]$ as

$$\begin{align*}
y[n] &= \sum\limits_{f=0}^{N-1} Y[f] e^{j2\pi n\,\frac{f}{N}},& n\in \{0,\dots,N-1 \}\tag{1}
\end{align*}$$
The DCT-II (which is most probably what we're looking at; the others are mathematically useful, but less nice to implement) $\mathbf{Y}[f]$:

$$\begin{align*}
\mathbf{Y}[f] &= \sum\limits_{n=0}^{N-1} y[n] \cos\left(\pi (k+\frac12) \frac nN \right) ,& f\in \{0,\dots,N-1\}\\
&\text{which, Euler says, is}\\
&= \frac12 \sum\limits_{n=0}^{N-1} y[n]\left(e^{j\pi (k+\frac12) \frac nN } + e^{-j\pi (k+\frac12) \frac nN }\right),& f\in \{0,\dots,N-1\}\text.\tag{2}
\end{align*}$$

Now, calculating the inverse discrete Fourier of the discrete cosine transform becomes but a matter of diligence; calculating by inserting $(2)$ into $(1)$, you'll notice that the complex sinusoids are orthogonal w.r.t. the scalar product, ie. the sum over the element-wise products becomes zero for everything but the same exponents.


Inter-symbol energy coding I have scoured my entire library and I can't find where I read about a method of using symbol energy intentionally smeared across multiple bit times in such a way that the symbols can be recovered correctly. Can someone point to any of the classic digital modulation books I am likely to have in my library where I would have read about this technique? I seem to recall it was a part of minimum shift keying techniques. It is frustrating that I can remember one or two of the figures but can't place where I read it.

Thanks.

Ray

The most likely possibility is that you're thinking of partial response signaling, of which GMSK is a special case. This is covered in many textbooks, e.g. in Proakis.


How do I resample an image to a rotated grid? I have an image, I, sampled on a uniform grid:

$\ x_i = i*\Delta x, y_j = j*\Delta y, $

I need to resample this image to a grid rotated counterclockwise by an angle $\ \theta$ around $\ (x_0,y_0)$:

$\ u_i = i*\Delta u, v_j = j*\Delta v $

$\ x(u_i,v_j) = x_0 + cos(\theta)*u_i - v_j*sin(\theta)$

$\ y(u_i,v_j) = y_0 + sin(\theta)*u_i + v_j*cos(\theta)$

How do I best do this?

I understand that I can e.g. do bilinear interpolation to find I from the 4 (x,y) corners around each $\ (u_i,v_j)$.

However I am not free to choose $\ \Delta u $ and $\ \Delta v $ as I wish.

What is the Nyquist criterion for this rotated grid?

Am I right in assuming that it will be?:

$\ \Delta u = cos(\theta)*\Delta x + sin(\theta)*\Delta y$

$\ \Delta v = -sin(\theta)*\Delta x + cos(\theta)*\Delta y$

You should read about image transformations. Anyhow, the idea is quite simple:

Compute the inverse transform $T(x,y)$.
For each point $(x,y)$ on the re-sampled grid:

Calculate $T(x,y)$, and find the 4 closest points on the original grid.
Apply interpolation technique (For example, bi-linear) on the neighbors, and put the value in the re-sampled grid.

Some references:

http://www.cis.rit.edu/class/simg782/lectures/lecture_02/lec782_05_02.pdf
http://www.imageprocessingbasics.com/geometric-transforms/


Discrete algorithm for low pass filter I am working on a position controller for a marine vessel. I have a measurement signal containing the y-position of the vessel that consists of both low frequency (<.1 rad/s) and high frequency (>.6 rad/s) motions. I am currently using a state observer that separates that combined motion in two separate states (a LF and WF states). The gains are calculated using the discrete time Kalman filter algorithm.
At every timestep I receive a new incoming measurement and my state observer corrects the predicted state with the measurement.
I want to compare the performance of this approach with using a low pass filter on the measurement signal instead of using an Kalman filter. The theory seems to be that a Kalman filter has a much lower phase lag, but I am seeing some amplifications in the magnitude bode plot of my outgoing LF signal compared to the incoming measurements and the filter thus amplifies motions that are between the LF and the WF frequencies...
I need some advice on how to implement a nth-order low pass filter in such a discrete time system. Is there an algorithm to filter the incoming measurements at every time step, or do i need to take the entire time trace of the measurement and filter that signal? I am using Matlab.
Any tips are more than welcome,
Scott

You can probably use Matlab's filter command. If you call it like this:

[Y,Zf] = filter(B,A,X,Zi);

the variable Zf contains the state of the filter after filtering, and Zi contains the filter's state before filtering. You could proceed roughly like this (in pseudocode):

Acquire several samples and store them in vector Samples
[Y,Z] = filter(B,1,Samples); % assuming FIR filter stored in B
while 1
    acquire new sample and store it in Sample
    [Y,Z] = filter(B,1,Sample,Z);  % Y contains the updated filter output
end

To design the filter, I'd start by looking at the fir1 command.


Why is second-order Gaussian called Laplacian Gaussian? We usually use Laplacian of Gaussian as the filter for edge detection or blob detection. But the filter itself is essentially a second order Gaussian. So why do we call it Laplacian of Gaussian? Is there some other meanings behind it?

Applying a Gaussian filter is a linear operation, also taking any type of derivative (finite difference) like Laplacian are also linear operations. Now considering the order of applying these operations doesn't matter (a property of convolution), you could combine them in many different ways, e.g. take Laplacian at first then take Gaussian or Gaussian first and Laplacian next, also you could apply Gausian on Laplacian or Laplacian on Gaussian and find a LoG Kernel, then apply the obtained Kernel on image once, which is more efficient in number of computations.


What is meant by "Group delay"?in simple words? Although likes of this question have been asked many times on DSP SE

But i am unable to understand those and i wish to have a crystal clear explanation in simple words with example

The only thing that i am able to understand from those questions is that group delay is negative derivative of phase with respect to frequency. But what does negative derivative means here and why not positive derivative here?

Here is my simplest explanation:
The group delay, as the negative derivative of phase, predicts the time delay of the amplitude envelope of a pulse, as shown in the hand-drawn graphic below. The upper part of the sketch shows a sinusoidal waveform varied in amplitude by its envelope. The lower one is showing this same envelope before and after a system that has group delay.
This applies when the phase of the frequency response can be approximated as linear for the "group" of frequencies within the pulse envelope. Thus for non-linear phase systems, this applies to generally narrower band signals such as the pulse I show where the amplitude transition is gradual.
Consider a single sine-wave with the amplitude envelope such as I show. The time delay of the sine wave itself would be predicted directly from the phase of the frequency response (by dividing by the frequency of the sine-wave: with $\phi = angle(H(j\omega))$,  the time delay is $-\phi(\omega)/\omega$), while the time delay of the pulse envelope is predicted from the negative derivative of the phase with respect to frequency ($-d\phi(\omega)/d\omega$)).

These posts and answers are helpful:
Calculate the time delay introduced by group delay for IIR-Filters
https://electronics.stackexchange.com/questions/135475/physical-significance-of-group-delay
And most helpful to what otherwise seems like a paradox of causality for positive group delay is this paper https://www.researchgate.net/publication/253463703_Causality_and_Negative_Group_Delays_in_a_Simple_Bandpass_Amplifier referenced by Max in this post Physical Meaning of Negative Group Delay for causal LTI systems which I bottom line as causality is not violated but due to the bandwidth restrictions above we create a condition that causes the pulse envelope at the output to precede the input: The output pulse does not appear until the input amplitude varies (and if we have gain in the system the output envelope grows faster) and due to the bandwidth constraint and the result of destructive summation of the input pulse the output will start to decrease before the input does. Very cool DSP magic trick. See this post that illustrates this with a specific example.


Build discrete IIR from a transfer function which is not in terms of $s$ Firstly I should state that my maths knowledge is limited.

I'm currently designing an acoustics modelling application which uses a rectilinear FDTD grid for modelling pressure variations - but I've become stuck upon trying to implement a fibrous layer boundary.  From Kowalcyzk 2008, the fibrous boundary model is defined as a frequency domain continuous impedance filter $\xi_w(\omega)$:

\begin{align*}
 \rho(\omega) = 1.2 + [-0.0364({{\rho f} \over {\sigma}})^{-2} - j0.01144({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}
\end{align*}

\begin{align*}
 K(\omega) = 101320{{j29.64 + [2.82({{\rho f} \over {\sigma}})^{-2} + j24.9({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}} \over {j21.17 + [2.82({{\rho f} \over {\sigma}})^{-2} + j24.9({{\rho f } \over {\sigma}})^{-1}]^{{1} \over {2}}}}
\end{align*}

\begin{align*}
 Z _w(\omega) = [\rho(\omega)K(\omega)]^{{1} \over {2}}
\end{align*}

\begin{align*}
 \Gamma(\omega) = j2\pi f[{{\rho(\omega)} \over {K(\omega)}}]^{{1} \over{2}}
\end{align*}

\begin{align*}
 \xi _w(\omega) = {{Z_w(\omega)} \over {\rho c}} \coth [\Gamma(\omega)d]
\end{align*}

Boundaries are implemented in the grid edge nodes as discrete time domain IIR filters of arbitrary order.  The few resources I've read on the subject have transfer functions that consist of a ratio of polynomials in terms of $s$ or $z$ - so I have no idea where to start on converting this filter, any clues?

I read this post, but I don't know if there's simpler solution for my circumstance.

There exist genetic and evolutionary algorithms that attempt to guess and optimize a set of pole/zero placements that approximate an arbitrary frequency response.

There is also a algorithm called Frequency Domain Least Squares : http://www2.units.it/ramponi/teaching/DSP/materiale/ES_5_2.m

The poles and zeros found can then be used to create an IIR recurance function.


Unable to understand how the paper simplifies the covariance matrix - Kalman filter The paper Convergence Analysis of the Unscented Kalman Filter for Filtering Noisy Chaotic Signals

presents the convergence analysis of Unscented Kalman Filter

download
http://www.eie.polyu.edu.hk/~cktse/pdf-paper/ISCAS07-Feng1.pdf . The Measurement update equations from Eq16 -- Eq19 are clear to me but I don't understand how Eq(20) is derived which is $P_n = K_nR_n$ where $P_n$ is the covariance of the filtering error and $R_n$ is the covariance matrix of the measurement noise. Can somebody please show the steps how Eq(20) is reached ? I tried something like this starting from Eq(18) but I am stuck --

$$ P_n = P_n^- - K_n P_{y_n y_n}K_n^T,$$
$$ = P_{y_n y_n} - R_n - K_n P_{y_n y_n}K_n^T  $$ (substituting Eq 19)
$$ = P_{x_ny_n} - K_n P_{y_n y_n}K_n^T   $$

In order to derive equation (20) you can use the following steps:

From substitution of  equation (16) into equation (18)

$$  P_n = P_{n|n-1} - K_nP^t_{x_ny_n}$$

Now plugging right side equation of (19) into the last equation

$$  P_n = P_{n|n-1} - K_n(P_{y_ny_n} - R_n)=$$
$$ =P_{n|n-1} - K_nP_{y_ny_n} + K_nR_n $$

Now Pluging (16) again we know that
$$ K_nP_{y_ny_n} = P_{x_ny_n}$$
putting this back we get

$$ P_n=P_{n|n-1} - P_{x_ny_n} + K_nR_n $$

Using left side equation of (19) we get the desired result
$$P_n = K_nR_n$$


Time derivative of signal - effect on noise distribution I have an angular velocity measurement that has a certain amount of ripple that yields an histogram shaped like a gaussian curve. I quantify that ripple using a standard deviation. I am interested in knowing more about the angle ripple.

What would be the relationship between the two ripples?

My guess would be that the derivative (from angle to angular velocity) push noises at higher frequency and that the gaussian curve will get a little wider (bigger standard deviation). How can I validate that and what mathematics are involved here?

Let's call $a[n]$ the samples of the angular measurement. Then, an estimate of the angular velocity (in angle/sec) is given by discrete differences, i.e., $v[n] = (a[n]-a[n-1])/t_0$, where $t_0$ is the sampling time.

Let's assume a signal pluse noise process so that $a[n] = a_0[n] + w[n]$. Inserting, we obtain $v[n] = v_0[n] + \bar{w}[n]$, where $\bar{w}[n] = (w[n]-w[n-1])/t_0$. Certainly, if $w[n]$ is zero mean, so is $\bar{w}[n]$ You can now compute the variance of your effective noise samples $\bar{w}[n]$, namely $$\mathbb{E}\{\bar{w}[n]^2\} = \frac{1}{t_0^2}\left(\mathbb{E}\{\bar{w}[n]^2\} + \mathbb{E}\{\bar{w}[n-1]^2\} -2 \mathbb{E}\{\bar{w}[n]\bar{w}[n-1]\}\right). $$ If your noise samples $w[n]$ are i.i.d. and have  constant variance $\sigma_w^2$, this leads to $\mathbb{E}\{\bar{w}[n]^2\} = \frac{2\sigma_w^2}{t_0^2}$. So there is a doubling of variance, but also some kind of rescaling. This depends of course on the unit of your velocity (whether it's angle/second or angle/sample, in the last case your variance is just $2\sigma_w^2$). NB: If they are correlated, this changes to $2\sigma_w^2(1-\rho)$, where $\rho$ is the correlation coefficient.

All this is true if the velocity is simply estimated from differences of adjacent samples. There are more elaborate schemes available, which would have an effect on the noise variance.


Fourier transform of discrete time unit step function To obtain fourier transform of u[n],
u[n] - u[n-1] = delta[n] , taking fourier transform of both sides of the equation results in :
U(w) - exp(-jw) U(w) = 1 , hence :
U(w) = 1/(1-exp(-jw)) which is wrong and the right answer has an extra term.
Which step is wrong in this possible solution?
I know the right proof of fourier transform of u[n], my question is regarding the wrong part of this solution.

The DFT of a unit step response is
$$U(\omega) = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)$$
Applying the shift property as you did will give:
$$\mathcal{F}(u[n] - u[n-1]) = U(\omega) - U(\omega)e^{-j \omega} = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega) - [\frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)]e^{-j \omega}$$
that is
$$\mathcal{F}(u[n] - u[n-1]) =  \frac{1- e^{-j \omega}}{1 - e^{-j \omega}} +\pi\delta(\omega)( 1 - e^{-j\omega}) = 1+\pi\delta(\omega)( 1 - e^{-j\omega})$$
The second term is always zero because for $\omega = 0$, $1 - e^{-j\omega} = 0 $ and it's zero on any other point, so you get,
$$\mathcal{F}(u[n] - u[n-1]) = 1 = \mathcal{F}(\delta[n])$$