Question about a diagram of PWM detection

From the diagrams, the integrator is reset to zero (sometimes this will be called "dumped", as in "integrate and dump") right after its value is sent to the hold.

The reason that I can see for doing this is that it appears that the leading edge of the PWM happens on a steady period, with the trailing edge varying to establish the width.  By sampling on the "A" pulse you insure that the timing of the PAM pulses is consistent; sampling on the "B" pulses would mean that the timing of the output would depend on the    amplitude; this would induce a nonlinearity.


The audible range is about 10 octaves, and usually the center frequencies of a graphic equalizer would be distributed equally spaced on a log scale to cover that range. Common equalizers have either $30$ bands (with $1/3$ octave filters) or $10$ bands (with $1$ octave filters).

If you want $5$ bands, you could choose filters that cover approximately $2$ octaves. The five center frequencies would then be something like

\begin{matrix}
30\,\text{Hz} & 125\,\text{Hz} & 500\,\text{Hz} &2\,\text{kHz} &8\,\text{kHz}
\end{matrix}

If you have an equalizer like that I'm not sure if an additional bass boost button makes much sense. If you want $6$ frequencies, just divide the $10$ octaves into $6$ equally spaced bands (on a log scale).


2D Convolution in the Spatial Domain vs Frequency Domain

To answer your question:

Can easily be done.
One must remember that the short signal (The Kernel) must be padded (With zeros) to have the same size as the image before the DFT conversion.
Once they have the same size all needed is to convert into the Frequency Domain have element by element multiplication and transform back.
A side note would be that this way you assume periodic boundary condition.
I'm not a C# / C Coder.
I can provide you a MATLAB reference if needed.
Just make sure you do the following steps:

Pad the kernel to the image size (Pad it with zeros to the right and left).
Convert into the Frequency Domain.
Multiply Element by element.
Convert back to the spatial domain.

2 Remarks:

In cases of large image relative to the size of the kernel you better (Efficiency wise) apply it in the spatial domain.
The method above describes to do Circular Convolution (See Applying Image Filtering (Circular Convolution) in Frequency Domain). In order to apply Linear Convolution See How Much Zero Padding Do We Need to Perform Filtering in the Fourier Domain?

Enjoy...


MATLAB: $\tt fft$ and $\tt ifft$ scaling

Whether or not to scale the forward FFT by 1/N depends on which result you want for further analysis: energy (preserving Parseval's identity), or amplitude (measuring height or volts, etc.).

If you want to measure or analyze energy, then don't scale by 1/N, and a longer sinusoid of the same amplitude will produce a larger FFT result, proportional to the greater energy of a longer signal.

Slightly more commonly, if you want to measure or analyze amplitudes, then to get a longer sinusoid (thus with more total energy at the exact same amplitude) to produce about the same FFT result as a shorter signal, you will need to scale down the FFT summation by an ratio proportional to the length.  The ratio could be reference_length/N, which is sometimes 1/N if the system input gain is 1.0 for whatever dimensions or units, including the time interval dimensions, that you choose to use in your further analysis.  You need to scale down proportionally because a DFT is a summation: the more you sum up similar items, the bigger the result.

So.  Energy or amplitude.  Which do you want?

Now if you scale down the forward FFT, then you should not scale the inverse so that IFFT(FFT(x)) == x.  Or do vice versa.

The 1/sqrt(N) for scaling seems to me to be for either when one needs a formal symmetry for some proof, or when building some sort of hardware pipeline where the latency and/or number of arithmetic units/gates for the DFT and for the IDFT needs to be identical.  But you get neither a good direct measurement of either energy or amplitude for any typical type of engineering analysis.


Why do we use normalized angular frequency?

We use normalized frequency units in DSP mainly for convenience-- most DSP textbooks are written this way so that the theory can be presented in terms of essentially unitless quantities.

For instance, if you are talking about a sampled signal $x[n] = x(nT_s)$, where $T_s$ is the sample period, the discrete-time signal $x[n]$ now has no notion of time associated with it, because $n$ is an integer. Consequently, there can be no notion of "physical" frequency, in Hz or rad/s.

Practically speaking, what this means is that you end up having to "translate" frequencies back and forth, i.e. find out what cutoff frequency you want your filter to be in Hz, make sure your sample rate is greater than twice that, then find out what that frequency is in normalized units.

For your original question, you just need to make sure that your sample rate is high enough to capture all the frequency content. Your sample rate would need to be greater than $2\times 660 \textrm{ Hz}$ at a minimum, but practically you may want to make it even higher for a number of other reasons.

In MATLAB your script would look something like this:

clear
clc
close all

Fs = 2*700; % set sample rate higher if needed
Ts = 1/Fs;
t_duration = 1; % for now, one second of audio. can change if needed
t = 0:Ts:t_duration;
A = 0.1;
signal = A*sin(2*pi*440*t) + A*sin(2*pi*550*t) + A*sin(2*pi*660*t);

Note that you might want to change the phases of the individual sine terms to get the effect you are looking for.

EDIT 1:

Note that in the example code $t = nT_s = n/F_s$, therefore the argument to the sin function is:
$$
2\pi f t = \left( 2\pi\frac{ f}{F_s} \right) n
$$

EDIT 2:

Technically, in your example if all you are doing is generating a signal, you don't need to explicitly use normalized frequency; you just need to make sure that your sample rate is high enough. This answer was touching on the broader question of, in general, why someone would want to.


Discrete Fourier Transform: Prove that only N distinct values of {X(k)} can be computed

It's not quite clear what you are supposed to consider as starting point for your proof. You have to make some assumptions of what is known before you can prove something else. So unless you give a more detailed description of what you want to prove from what, I will give you a proof that uses a sensible set of assumptions.

Let's characterise the discrete Fourier transform as a map between two complex vector spaces $$\mathcal{F}:\mathbb{C}^N\to\mathbb{C}^M$$ of dimension $N$ and $M$ respectively. The map is supposed to preserve the structure of the vector space, so it must be linear.

A linear map between finite dimensional vector spaces can be completely described using a matrix, which in this case is $M\times N$. The rank of this matrix $F_+$ is less or equal to the smaller dimension, i.e. $\mathrm{rank}(F_+)\le\min(N,M)$.

We also require the map to be a transform, meaning that it can be inverted. So $$\mathcal{F}^{-1}:\mathbb{C}^M\to\mathbb{C}^N$$ must exist and $$\mathcal{F}^{-1}\circ\mathcal{F}=\mathrm{id}_N$$ as well as $$\mathcal{F}\circ\mathcal{F}^{-1}=\mathrm{id}_M$$
and the matrix $F_-$ of the inverse transform is also less or equal to the smaller dimension, i.e. $\mathrm{rank}(F_-)\le\min(N,M)$.

We can write the inversion conditions as a matrix equation
$$F^{-1}\cdot F = E_N$$
$$F\cdot F^{-1} = E_M$$
where $E$ is the identity matrix of the size specified in the subscript. The matrix multiplication can be translated to a rank inequality:

$$\min\left(\mathrm{rank}(F_-),\mathrm{rank}(F_+)\right)\geq\mathrm{rank}(E_N)=N$$
and
$$\min\left(\mathrm{rank}(F_+),\mathrm{rank}(F_-)\right)\geq\mathrm{rank}(E_M)=M$$
because the product of two matrices has at most the rank of the matrix with the lesser rank.
From above we also had
$$\mathrm{rank}(F_-)\le\min(N,M)$$
and
$$\mathrm{rank}(F_+)\le\min(N,M)$$

These four inequalities only allow for one solution, which is $M=N$. To see this, assume that for example $N<M$, then it follows that $\mathrm{rank}(F_-)=\mathrm{rank}(F_+)\leq N$ but also $\mathrm{rank}(F_+)\geq M$, from which $N\geq M$ would follow, making the assumption inconsistent. The same goes for assuming $N>M$.

That means the discrete Fourier transform must be between complex vector spaces of equal dimensions. Or in other words, the number of frequency bins must be equal to the number of time samples.

The proof above can of course be a single line if you already know that vector space isomorphisms require matching dimensions. Of course, some linear algebra is still necessary, but the rank formulas should be rather intuitive.


Perhaps standard deviation is what you want. Standard deviation
$$\sigma= \sqrt{\frac{1}{n}\sum_{k=1}^n (x_k-x_{ave})^2}$$
(where $x_{ave}$ is the average of your sequence) quantifies the fluctuations of your sequence above and below the sequence's average value. The square of the standard deviation, $\sigma^2$, is the variance of your sequence which is a measure of the power of the fluctuations.


One easy way to see the validity of the theorem is to consider the corresponding analytic signals. Let a signal $s(t)$ be defined by
$$s(t)=m(t)c(t)\tag{1}$$
where $m(t)$ is a low pass (message) signal, and $c(t)$ is a bandpass (carrier) signal. Note that we don't need to assume that $c(t)$ is sinusoidal. Let $\hat{x}(t)$ denote the Hilbert transform of $x(t)$, and $\hat{c}(t)$ is the Hilbert transform of $c(t)$.
The equality claimed by Bedrosian's theorem
$$\hat{s}(t)=m(t)\hat{c}(t)\tag{2}$$
is equivalent to
$$s_a(t)=m(t)c_a(t)\tag{3}$$
where $s_a(t)$ and $c_a(t)$ denote the analytic signals corresponding to $s(t)$ and $c(t)$, respectively. Analytic signals are complex-valued signals having only positive frequency components. The analytic signal is given by the original signal plus its Hilbert transform as the imaginary part:
$$s_a(t)=s(t)+j\hat{s}(t)\tag{4}$$
In the frequency domain, the signal $m(t)c_a(t)$ is centered at the (positive) center frequency of the bandpass signal $c(t)$. Its bandwidth equals the bandwidth of $c(t)$ plus twice the bandwidth of the lowpass signal $m(t)$. It has only positive frequency components if the spectra of $m(t)$ and $c(t)$ do not overlap. In that case it equals the analytic signal $s_a(t)$, as claimed in Eq. $(3)$. Otherwise, $m(t)c_a(t)$ has also negative frequency components and $(3)$ is not satisfied.
Note that $(2)$ and $(3)$ are also valid if $c(t)$ is a highpass signal, as long as the spectra of $m(t)$ and $c(t)$ don't overlap.
Also take a look at this related answer and the figure therein.


The small red peak in the CIE standard observer

Probably you've noticed that primarities are $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$, not $\mathbf{R}$, $\mathbf{G}$, $\mathbf{B}$ (which are corresponding to the color values $R$,$G$,$B$). This is the aftermath of original work conducted by Wright and Guild yielded the Color Matching functions: $r(\lambda)$, $g(\lambda)$, and $b(\lambda)$ having a negative value of the red primarity around $522 \;\texttt{nm}$ (figure below). Everyone will agree that it is non-intuitive and somewhat confusing.

Why values of $r(\lambda)$ are negative? Well after tristimulus color matching experiments, it turned out that not every color can be created with the preset primarities ($700\; \mathtt{nm}$, $564.1\; \mathtt{nm}$, and $435.8\; \mathtt{nm}$). So for example, if you wanted to create the blue-green color of wavelength $\lambda = 490 \; \mathtt{nm}$ it cannot be done only by adding blue and green primarity alone. You also need to subtract the red one. Mathematically speaking, the color equation:

$$\mathbf{C}=r\mathbf{R} + g\mathbf{G} + b\mathbf{B} $$

Must by modified by adding the outer mixture $-r_{out}\mathbf{R}$ to both side of the above equation:

$$\mathbf{C} + r_{out}\mathbf{R} = g\mathbf{G} + b\mathbf{B} $$

That is why CIE introduced new set of so-called virtual primarities: $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$. In order to obtain only positive or null values of the Color Matching Functions. That's how $2^{\circ}$ observer Standard Color Matching Functions (also called CIE 1931) were created. By doing so they ensured that (amongst other assumptions):

New SCMF's  $x(\lambda)$, $y(\lambda)$, $z(\lambda)$ have values $\ge 0$ (not like $r(\lambda)$).
Are related to previous CMF's ($x$,$y$,$z$) by a appropriate transform.
$\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$ must be independent (kind of orthogonal).
Only color value $Y$ is proportional the the lightness of a colour.

Newly obtained color quantities $X$,$Y$,$Z$ are called Standard Color Values. The relation between CMF's and SCMF's is given by a following linear transform:
$$
\left( \begin{array}{c}
x(\lambda) \\
y(\lambda) \\
z(\lambda) \\
\end{array} \right)
=
\left( \begin{array}{ccc}
2.7689 & 1.7517 & 1.1302 \\
1 & 4.5907 & 0.0601 \\
0 & 0.0565 & 5.5943
\end{array} \right)
\cdot
\left( \begin{array}{c}
r(\lambda) \\
g(\lambda) \\
b(\lambda) \\
\end{array} \right)
$$

After over a three decades people from CIE had realised that $2^{\circ}$ observation angle might not be enough for adequate color assessment. That's why in 1946 they introduced the $10^{\circ}$ observer, also known as CIE 1964 observer. Latter is recommended for visual angles above $4^{\circ}$, whereas CIE 1931 should be used for any angle below that. Of course there is no sudden change in perception of colour at $4^{\circ}$, this number is just arbitrary. Comparison of these two standards is depicted of a following plot, on which you can see the 'bump' mentioned by you for the $x(\lambda)$ standard color value:


PIsarenko Harmonic Decomposition method for PSD estimation

Since you have $v_{min}$, the eigenvector corresponding to the minimum eigenvalue, compute the function $\hat{P}_{PHD}(e^{j\omega})$.

The vector $e$ is given by $e = [1 \quad e^{j\omega} \quad e^{2j\omega} \ldots e^{pj\omega}]$

Make a plot of $\hat{P}_{PHD}(e^{j\omega})$ and measure the $p=3$ peaks to find the frequencies.


Consider a 1D wave of 10sec duration.

I correct in assuming that the fundamental frequency is the image itself?
Does the 10 sec sample equal the fundamental freq? No. same for images.

And one copy of the image is equivalent to 1 wave cycle?
Repeat above. Answer: No

The sampling frequency is the number of pixels in the image?
One pixel can be thought of as the sampler. But, this could be changed. You could go for sub-pixel sampling also.

With this I think the last question has no relevance. Comparing two images, to the best of my knowledge is an open ended problem.


To insert a JPEG image into another without recompression, you can use a jpegtran version with crop 'n' drop. It is safest to insert into coordinates that are multiples of the 16x16 pixel JPEG block size. On command line:

jpegtran.exe -drop +16+16 pixelated.jpg piano.jpg out.jpg

pixelated.jpg

piano.jpg

out.jpg

Difference images made in Gimp show only black, so there was no lossy recompression.

Difference image: piano.jpg vs. out.jpg

Difference image: pixelated.jpg (aligned) vs. out.jpg

Before all this, you can use

jpegtran -crop

to copy the target region in the big image to a file for modification.

I wonder if jpegtran can handle your huge files. My test files were quite small.


Increasing the number of points in the frequency spectrum

You can also interpolate as many points as you need using a Sinc or windowed Sinc interpolation kernel.


What does "invariant to bias and gain changes" mean?

You need to be careful when reproducing formulas! Your $i_r$, $i_w$ are actually $\bar{i_r}$, $\bar{i_w}$ in the paper and two sentences above $(4)$, it says these are

  their zero-mean versions, which are
  obtained by subtracting from each vector its corresponding
  arithmetic mean.

So, the symbols in your formula are always the same, regardless of whether there is a constant value added to them (that's the definition of bias), because that would automatically shift the arithmetic mean by the same amount.

The division by the individual vector's norm also "undos" any multiplication with a scalar (that's the definition of gain).

So, the overall formula is unchanged (that's what invariant means) when you add a bias to the vectors, and unchanged when you apply gain to the vectors.


What is the difference between a one dimensional and a two dimensional signal?

Basically the dimensional number refers to the number of independent variables (input).

In one dimensional signal $f(t)$, the amplitude $y=f(t)$ is the dependent variable (output), and there is only one independent variable $t$.

In two dimensional image $f(x,y)$, the two independent variables are $x$ and $y$, and the dependent variable $f$ is also called intensity. If you want to make it similar with amplitude vs time plot as you do with one dimensional signal, you can plot a intensity vs x vs y scattered dots with three axis as well.


FFT spectrogram in log frequency space. How?

You can use either interpolation between FFT bins, or something like an MFCC triangular filter bank, except using equal delta_log spaced triangles instead of Mel-spaced triangles, or both.


Is there any transformation to exploit the sparsity of a Gaussian Wave?

Use a dictionary of Gaussian wave atoms, at all different locations and scales (or frequencies). A specific Gaussian should matching maximally with one of the atoms.

This is a generic answer. More details to come if you tell more about your wave: definition, dimension, and purpose of the sparsifying transform


What does it mean in image processing when a filter is called non-linear?

A filter F is called "linear", iff for any scalars $c_1$, $c_2$ and any images $I_1$ and $I_2$:

$F\left(c_1\cdot I_1+c_2\cdot I_2\right)=c_1\cdot F\left(I_1\right)+c_2\cdot F\left(I_2\right)$

This includes:

Derivatives
Integrals
Fourier transform
Z-Transform
Geometric transformations (rotate, translate, scale, warp)
Convolution and Correlation
the composition of any tuple of linear filters (i.e. applying some linear filter to the output of another linear filter $F(G(I))$)
the sum of the result of any two linear filters (i.e. the output of one filter, added pixel by pixel to the output of another filter $F(I) + G(I)$)

and many others.

Examples of non-linear filters are:

the square, absolute, square root, exp or logarithm of the result of any linear filter
the product of the result of any two linear filters (i.e. the output of one filter, multiplied pixel by pixel with the output of another filter $F(I)\cdot G(I)$)
morphological filters
median filter


This serves as potential starting points only, on how to transpose classical signal processing tools to specific kinds of binary data.

First, remember that standard spectrograms, fed as inputs to classifiers, can be seen as means to extract attributes from non-stationary data. There is apparently little need for inversion, so we are interested in analysis only (no synthesis).

Here, we have sequences of 8-bit words $w_0w_1w_2\ldots w_{N-1}$. When dealing with binary patterns, the frequency domain might not be the most appropriate. In that case, as with communication systems, the related notion of sequency could be useful. It relates to blocky oscillating patterns of $\pm 1$ (or $0/1$):

The number of zero-crossing increases in a manner similar to augmenting frequencies in Fourier transforms. I would thus prefer to use those vectors in sliding windows to analyze your binaries. You can find literature in such functions as Walsh, Paley or Hadamard (or WaLeyMard) sequences, depending on their ordering. The proper choice of the transform length is not evident.  This seems to exist already, under the name: "Walsh spectrogram": On Generating Walsh Spectrograms, 1976.

I would think it might depend a bit (pun intended) on the compiling or storage system. There is a possibility to look at dyadic size factors with Haar wavelets or Haar wavelet packets.

Note that there exist transformations for data taking values in finite fields, from The Fast Fourier Transform in a Finite Field, 1971, J. M. Pollard.

Other references:

Harmut, H. F., 1970: Transmission of information by orthogonal functions
Agaian, S. S., 1984: Hadamard Matrices and Their Applications


Measure similarity between quantized grayscale images

This question is very broad, yet interesting. Your way to go are co-occurence matrices.

Co-occurence matrices can be calculated for each image and provide further indicators which specifically allow identification of different textures.

Note that there does not exist a single co-occurence matrix for an image, but several ones, since the co-occurence matrix is defined for specific directions of occurence. Further, you can vary the granularity of the co-occurence matrix. For each co-occurence matrix you can derive several features which enable you to distinguish the textures.

This example is very simple, yet should enable you to understand:

Image you have black-white picture of only a single line of pixels such as:

[0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0]

Now there are there are four possible co-occurences of pixel neighbours:

1) [0,0]
2) [0,1]
3) [1,0]
4) [1,1]

For each these we will now count the times that you can them within the strings:

1) [0,0] ---> 4
2) [0,1] ---> 3
3) [1,0] ---> 3
4) [1,1] ---> 0

This will give us our co-occurence matrix:

[[4,3],[3,0]] which is specific for each pattern.

Based on this co-occurence matrix, you are able to calculate features, such as Contrast, Correlation, Energy, Homogenity and many more.

Matlab provides you with tools to get the co-occurence matrix. Further, wikipedia provides an introduction.

Once you implemented this and encountered your first problems, I am looking forward to your next question. Good luck.


Is match filter a solution of some differential equation?

As the name implies, a matched filter is found to be an exact replica of the signal of interest to be detected; just a reversed version as a result of an optimization.

Hence your question reduces to whether the signal of interest is produced as a solution of some diferential equation.

Even though it's true that there are signals as solutions to differential equations of all sorts (linear, nonlinear, time varying, time invariant etc), whether an arbitrary signal can be imposed to be a solution to some solvable differential equation is not very clear to me.

Since the signal of interest is an arbitrary signal, you must show that every (arbitrary) signal can be shown to be a solution of a corresponing differential equation. I would say no.

Note that as a degenerate case $y(t) = x(t)$ can also be considered to be a (useless) differential equation (of order zero?) but I think it's not what you mean by a differential equation.


why there is no 0 on the right side? (matlab, low pass filter)

Answer: Your derivation is absolutely correct and You will see 0 for both $\omega = -\pi$ and $\omega = \pi$, if you use fvtool() function of MATLAB instead of fft(). The reason is simple : FFT does not calculate values of $H(e^{j\omega})$ at continuous $\omega \in [-\pi, \pi]$, but it calculates DFT and for that the digital frequency resolution is $\omega = \frac{2\pi k}{N}$. So, for $N=6$, it will compute values of $H(e^{j\omega})$ at $\omega = 0, \frac{2\pi}{6}, \frac{4\pi}{6}, \frac{6\pi}{6}, \frac{8\pi}{6} \ and \ \frac{10\pi}{6}$. See that $\omega = \pi$ is occurring only at $k=3$, and hence you are seeing $0$ only once.

For checking this, you can take any symmetrical h, like $h = {1,2,4,4,2,1}$.
Then use the following MATLAB code to plot its DTFT for $\omega \in [-\pi, \pi]$

  fvtool(h,1);

This will plot a Single sided DTFT from $\omega \in [0, \pi]$. You can then change from "Analysis" option to plot for frequency range $[-\pi, \pi]$.

You will get following plot for the example I have taken:

Check that the values are $0$ at both $\omega = -\pi$ and $\omega = \pi$.

And this will happen for any symmetric values, not any particular window function.


Understanding the $\mathcal Z$-transform

Consider a liner discrete-time system. Assume we can define it in terms of an input-output relation as follows (you can assume a more general model but it is enough for our purpose):

$$a_0y[n]+a_{1}y[n-1]+\cdots+a_{N}y[n-N]=b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]\tag{1}$$

When the coefficients $\{a_i\}$ and $\{b_i\}$ are constant, we call it a finite-order constant-coefficient ordinary difference equation.
It explains the current output $y[n]$ in terms of a weighted sum of the current and past inputs and the past outputs:

$$y[n]=\frac{-a_{1}y[n-1]-\cdots-a_{N}y[n-N]+b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]}{a_0}$$

It is very similar to a differential equation in continuous time.

Dealing with such equation to find the output of the system can become complicated. It is a recurrence relation. It defines the values recursively and for arbitrary inputs it is not straightforward to express the output in a closed-form representation.

How to deal with it in an easier way?

Consider the following transform (assume the sum exists):
$$\mathcal{Z}(x[n])=\sum_{n=-\infty}^{+\infty}x[n]z^{-n}$$
You are right, it is just a transform that accept $x[n]$ and gives you $X(z)$.
But it has a useful property. Lets calculate the $z$-transform of $x[n-\alpha]$:
$$\begin{align}\mathcal{Z}(x[n-\alpha])&=\sum_{n=-\infty}^{+\infty}x[n-\alpha]z^{-n}\\[10pt]
&=\sum_{n'=-\infty}^{+\infty}x[n']z^{-(n'+\alpha)}\tag{2}\\[10pt]
&=\left(\sum_{n'=-\infty}^{+\infty}x[n']z^{-n'}\right)z^{-\alpha}\\
&=\left(\mathcal{Z}(x[n])\right)z^{-\alpha}
\end{align}$$
where in $(2)$ I changed the variable $n'=n-\alpha \Rightarrow n=n'+\alpha$.
Assume $\mathcal{Z}(x[n])=X(z)$. We have seen that a property of $z$-transform is: $$\boxed{\mathcal{Z}(x[n-\alpha])=z^{-\alpha}X(z)}$$

Applying this to the difference equation $(1)$ makes it (note that it is a liner transform and can be applied term by term):

$$a_0Y(z)+a_1z^{-1}Y(z)+\cdots+a_Nz^{-N}Y(z)=b_0X(z)+b_1z^{-1}X(z)+\cdots+b_Mz^{-M}X(z)$$

Note that all terms are expressed in terms of $Y(z)$ and $X(z)$. So we can say that

  $z$-transform reduces usually difficult-to-handle recurrence relations to much easier to manipulate algebraic relations.

With an algebraic relation, we can factor the terms:

$$Y(z)\left(a_0+a_1z^{-1}+\cdots+a_Nz^{-N}\right)=X(z)\left(b_0+b_1z^{-1}+\cdots+b_Mz^{-M}\right)$$
and consequently,

$$Y(z)=X(z)\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$

where $$H(z)\triangleq\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$
is called the system function. So we distilled all the system into an algebraic expression as $H(z)$. Hence, the zeros and poles of $H(z)$ have a direct impact on the output.

There are also other reasons that make $z$-transform and $H(z)$ important, and we would prefer working with the $z$-domain rather than time-domain.

Assume the input to the system is $x[n]=z_{k}^n$, where $z_k$ is just a complex number (don't mix it with the $z$ variable in the $z$-transform at this moment! Also there are some restrictions on $z_k$ as we will see next).
The output of the LTI system can be calculated by convolution of the input and the impulse response $h[n]$:

$$\begin{align}
y[n]&=h[n]*x[n]\\
&=\sum_{m=-\infty}^{+\infty} h[m] x[n-m]\\
&=\sum_{m=-\infty}^{+\infty} h[m] z_k^{n-m}\\
&=z_k^n \sum_{m=-\infty}^{+\infty} h[m] z_k^{-m}\\
&=z_k^n H(z_k)
\end{align}$$

What does it really mean?

The output of the system is an scaled version of the input. So the signal $z_k^n$ passes through the system and the only impact on it is a scaling. Note that the scaling factor is the system function when evaluated at $z_k$. Why is it important?

If we could decompose an arbitrary signal (not of the form $z_k^n$) into several components of the form $z_k^{n}$, then we can easily calculate the output of the system to each of the components, and then add them together to find the overall output (benefiting from the superposition property of LTI systems).

How can we decompose the signal $x[n]$ into components of the form $z_k^n$?
It is easy as the following (assume the sum exists for now)

$$\left(\sum_{n=-\infty}^{+\infty}x[n]z_k^{-n}\right)z_k^n$$

Let's refer to this component as $X(z_k)z_k^n$. So the output to this component is $y_k[n]=H(z_k)X(z_k)z_k^n$. We can have as many components as possible. Potentially, we have thousands of $z_k$. So we refer to all of them by a variable $z$ and the component-wise output becomes $H(z)X(z)z^n$. The overall output $y[n]$ is in form of an integral (since we need to add them up). It can be shown that $y[n]$ is given by the following contour integral:

$$y[n]=\frac{1}{2\pi j}\oint_C\frac{H(z)X(z)z^n}{z}dz$$

where $C$ is a counterclockwise path inside the ROC (where $X(z)$ and $H(z)$ exist) that encircles the origin, which is the same as the inverse $z$-transform of $X(z)H(z)=Y(z)$.

Regarding your second question, in the definition of $z$-transform (when it exists on the complex unit circle) if you choose $z=e^{j\omega}$ it becomes the definition of the Discrete-Time Fourier transform Hence is the frequency-domain interpretation.


Can somebody explain SDF and MACE to me?

Your question is allusive at that point: explanation of acronyms (SDF, MACE), need for composite filters... A starting point is  matched spatial filter (MSF). To overcome its limitations, several versions have been proposed. David Casasent is a common denominatoir of the ones you consider.

For the SDF filer, my reference is:

C. F. Hester and David Casasent, “Multivariant technique for multiclass pattern recognition,” Applied Optics, 1980

  A technique for multiclass optical pattern recognition of different
  perspective views of an object is described. Each multiclass
  representation of an object is described as an orthonormal basis
  function expansion, and a single averaged matched spatial filter is
  then produced from a weighted linear combination of these functions.
  The technique is demonstrated for a terminal missile guidance
  application using IR tank imagery.

with comment in:

J. B. Nelson, "Multivariant technique for multiclass pattern recognition: comment"

  In a recent article, Hester and Casasent discussed an interesting
  technique in multiclass pattern recognition. They mentioned that other
  approaches to the construction of an orthonormal basis function set
  were possible and specifically suggested the possibility of using a
  Karhunen-Loéve expansion as an alternative to their own Gram-Schmidt
  procedure. It is the purpose of this Letter to recommend yet another
  alternative, the singular value decomposition (SVD), which has very
  convenient conceptual and numerical advantages, and to show how it may
  serve as a foundation for future extensions of the theory. The scope
  of the present discussion is restricted solely to the digital
  processing aspects of Ref. 1. It will begin with a brief review of the
  problem in the slightly modified matrix notation that is more
  compatible with the SVD technique that follows.

Abhijit Mahalanobis, B. V. K. Vijaya Kumar, and David Casasent, Minimum average correlation energy filters, Applied Optics, 1987

  The synthesis of a new category of spatial filters that produces sharp
  output correlation peaks with controlled peak values is considered.
  The sharp nature of the correlation peak is the major feature
  emphasized, since it facilitates target detection. Since these filters
  minimize the average correlation plane energy as the first step in
  filter synthesis, we refer to them as minimum average correlation
  energy filters. Experimental laboratory results from optical
  implementation of the filters are also presented and discussed


Large spike at beginning of signal after applying Bandpass filter

What you see as a spike at the beginning of the filter output is the impulse response of the bandpass filter itself. This would happen as a transient effect and it will be more pronounced if the filtered signal is sinusoidal + dc in nature. If it were a pure zero mean noise than the spike would be somewhat obscured.

You can either use a group delay shift on the output to extract the relevant portion without that spike, or try to set the initial conditions of the filter to disable the transient reponse but almost always the first approach is preferred.


Finding remaining zeros of linear phase FIR filter

For linear phase FIR filters, each zero at z = z0 will have a matching reciprocal zero at z = 1/zo. And for real-valued coefficients each zero at z = zo will have a matching conjugate zero at z = *zo. Thus for linear phase real-valued coefficients, when you place one zero on the z-plane you determine the location of the other three zeros.


This depends a lot on whether your constraints on the number of taps is computation or memory.
If you have a enough memory, you can simply use an FIR with all taps zero except for the last one (Which is $1$). That's equivalent to using a ring buffer and no filter at all.
If you don't have enough memeory, things are more tricky. YOu can get delays that are constant with frequency with parallel/cascaded alpass filters and you can use "schroeder allpass" filters to crank up the group delay effificently. The best way to do this depends a lot on your specific constraints: how much memory do you have, how wide is your frequency range and how "constant" does that delay needs to be.


Help me understand the phase of an FIR filter

Your answer does not take into account the possible sign of the magnitude term in the polar description of the frequency response $H(\theta) = A(\theta)e^{j\phi(\theta)}$. I guess you mix this definition with the more usual definiton of  $H(\theta) = |A(\theta)|e^{j\phi(\theta)}$ where the magnitude will always be positive and therefore phase will be what you predicted, however for your example $A(\theta)$ may assume a negative value and you shall reflect this as an addition of $\pi$ into the phase where $e^{j(\theta + \pi)} = -e^{j\theta} $


Linear Convolution in DSP and Hann window

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.


How do bode plots work with unstable systems work?

The frequency response can exist, even if there are poles in the right half-plane, as long as there are no poles on the $j\omega$-axis. You cannot use the unilateral Laplace transform, instead you must use the bilateral Laplace transform (with the lower integration limit at $-\infty$), or, simply, the Fourier transform.

The corresponding time domain function is left-sided (or two-sided, if there are also poles in the left half-plane), i.e., non-causal. Note that poles in the right half-plane only indicate instability if you restrict yourself to causal systems. Non-causal systems can be stable if there are poles in the right half-plane.

Take as an example two simple first order systems, one with a pole in the left half-plane, and one with a pole in the right half-plane, such that it is a mirror image of the left half-plane pole:

$$H_1(s)=\frac{1}{s+1}\tag{1}$$
$$H_2(s)=\frac{1}{s-1}\tag{2}$$

If $H_1(s)$ is the frequency response of a causal system (note that we can decide that in this case), then it is stable, and its frequency response is simply

$$H_1(j\omega)=\frac{1}{j\omega+1}\tag{3}$$

In the same way, if we decide that $H_2(s)$ describes a non-causal but stable system, then its frequency response exists and is given by

$$H_2(j\omega)=\frac{1}{j\omega-1}\tag{4}$$

Note that their Bode magnitude plots are identical, because

$$\big|H_1(j\omega)\big|=\big|H_2(j\omega)\big|=\frac{1}{\sqrt{1+\omega^2}}\tag{5}$$

Their phase plots are of course different.


The Dirichlet conditions provide sufficient conditions under which a real periodic function can be assimilated to its Fourier series at each point of continuity.

One of these conditions is that the functions have a finite number of extrema  (maxima or minima) in every bounded interval. This means that, in some sense, the function does not oscillate infinitely on some time segment, even with very tiny oscillations. An example of such an oscillating continuous function can be given by:

$$ f(t)=t^2 \cos (1/t)$$

with $f(0)=0$ by continuity. The ripples around $0$ get narrower and more numerous, and are infinitely many, in any non empty interval containing $0$.

The function is evidently continuous, even differentiable. On each interval defined by $$\left[\frac{1}{(k+1)\pi+\pi/2},\frac{1}{k\pi+\pi/2}\right]$$
the function $f$ vanishes at both ends, and since not identically $0$, it attains at least one extremum on this interval. Hence the function has an infinity of extrema on each interval $[0,a]$, $a>0$. The exact locations of the extremas can be computed with the derivative of $f$.


Finding the impulse response, given a differential function involving the input and output

If you do not select
$$
y'(0) = 0\\
y(0) = 0
$$
then the system will not be a linear, time-invariant filter, so one usually makes that selection and then deals with the special cases for when the initial conditions are not zero.


If low frequency travels longer distance, why our speech is not travelling longer distance?

Well, first of all the Sound Level Pressure decreases by $6 \; \mathtt{dB}$ when doubling the distance - this plays a big role. We do also have sound attenuation coming from our medium - air. Let's take a closer look onto sound absorption coefficient for different frequencies:

Knowing that human speech is mostly concentrated at the range of $300\;\mathtt{Hz}$ - $3.3\;\mathtt{kHz}$ (yes I know - let's forget about sibilants) you can judge that we are really doing our best. Larynx tones and first formants of voiced speech are clearly transmitted with minimum lost of energy. Therefore main problem is the decrease of SPL with increasing of distance from source.

Another issue is the hearing curve of human ear which makes the low attenuation at the low frequencies worse. For example even though air absorption for $2 \;\mathtt{kHz}$ - $3 \;\mathtt{kHz}$ is rather high, we can still hear the ambulance siren as it is in frequency range where human ear has it's resonance. Obviously that depends on many factors, i.e. in city you have lots of barriers for sound.

For example let's say that you are a male speaking with normal level. Exemplary average $dBA$ values in frequency ranges are following:

$250\; \mathtt{Hz}$: $51 \; \mathtt{dBA}$
$1\;\mathtt{kHz}$: $49\;\mathtt{dBA}$

Assume that these are values measured at the distance of $1\mathtt{m}$ We want to know how loud you are at $100\mathtt{m}$. Applying atmospheric attenuation factor for 10% of relative humidity gives:

$L_{250}^{a}=51-0.1=50.9\;[\mathtt{dBA}]$

$L_{1000}^{a}=49-1.4=47.6\;[\mathtt{dBA}]$

Not a big change, isn't it? Now let's apply the following formula for calculating the sound level at a given distance $d_{ref}$ (assuming cylindrical wave):

$$L_d=L_{ref}-20\log{\left( \dfrac{d_{ref}}{d} \right) } $$

Applying it to previously calculated values and substituting $d=1\;\mathtt{m}$, $d_{ref}=100\;\mathtt{m}$ we get:

$L_{250}^{a,d}=L_{250}^{a}-40=10.9\;[\mathtt{dBA}]$

$L_{1000}^{a,d}=L_{1000}^{a}-40=7.6\;[\mathtt{dBA}]$

Both values we've obtained are way below the level of human breathing at a distance of $1\;\mathtt{m}$  ($25\; \mathtt{dBA}$) and environmental noise. You can also try to shout but it still won't give you  big increase in a distance. Please refer to this article: Average Speech Levels and Spectra in Various Speaking/Listening Conditions for detailed info about levels of speech.

One last word, bit out of scope of this question. There are some cases when you can be heard in far distance or you could be heard very poorly in a very close distance. Such phenomena occurs a lot in water (but also in the air), where sound waves can travel tens of kilometres. This is due to presence of wind, temperature inversions, terrain profile and refraction of sound. If you wish to know more, then I am sure there are plenty of guys in Physics SE willing to answer that question. As for start, here is some interesting article: Outdoor Sound Propagation in the U.S. Civil War, as well as bunch of papers regarding Outdoor Sound Propagation.


This appears to be just a matter of projection. Try "imagesc" in matlab.


You need to determine an expression for the samplers. Normally i would think that those are A/D converters, but  they give $A(s)$ and in the variable $s$ it is implied (sort of) that is the laplace transform, and laplace transforms are used for continuos time systems, so the samplers T could also be just converting the signal into a piecewise constant signal, and that is a non linear operation.

Please try to clarify this, and then we will be able to help you more easily.


Strange Response of a fir1 Filter in Matlab

Since you are seeing your filter coefficients and much larger than your signal, this is indicative that you have a very large "impulse" at the start of your signal. The coefficients of your filter is the impulse response for the filter, so that is exactly what you are seeing: the response to an impulse.

Review the start of your time domain data for a very large offset in the first sample or near the start of the file.


Performing image equalisation and then Jpeg compression

I did some experiments. Equalizing, then compressing, resulted in a nearly 50% larger file.  Normalizing increased the filesize but not so much.  Normalizing, then equalizing, also increased the filesize by nearly 50%. All are JPEG quality 50, created with ImageMagick:

 - left: original (4540 bytes); right normalized (5467 bytes):

- left: equalized (7444 bytes); right: normalized, then equalized (7251 bytes)


How do I perform the sifting process in empirical mode decomposition?

It's not shifting, it's sifting. After you find the $m_1$ curve, your subtract it from your signal $X$ to get $h_1= X(t) - m_1$. You're essentially taking away that signal component and are left with $h_1$, which is your new starting point of your next iteration. I assume you were confused by the terminology.

Sifting is a general term in signal processing related to separating out components of a signal one at a time. It is frequently used in context of wavelet decomposition since this process is very similar in this respect.

The reason the next function is $h_{11}$ and not $h_2$ is because you have to keep going with this process until the resulting $h_{1k}$ is an IMF function (criteria given on the top of the page). Once you've reached that point, you get your first IMF basis function $c_1$. Only at this point have you reached your first empirical mode.

Then you start over by finding $m_2$ and so on.


Recently I have seen the paper Hui Yin, Yuanhao Gong, Guoping Qiu - Side Window Filtering.
They suggest a really simple filtering framework for Edge Preserving Filter:

Basically, what they suggest is filtering the image with a set of filter based on the Box Filter.
This filter set is basically composed of 8 filters with different orientations and sub sets of the Box Filter (As seen in the figure above).
Once you apply all filters (Which each of them can be implemented very efficiently by all the efficient implementations available for Box Filtering) you chose, on a per pixel basis, the one most similar to the original pixel.
Applying it by iterations, yields very efficient and very good filter.
I took the Lena Image:

I applied 40 iterations of the filter (Link to 75 Iterations):

The full code is available on my StackExchange Signal Processing Q74674 GitHub Repository.


The in-phase and quadrature components of a bandpass signal are real-valued baseband signals. So if the spectrum of a bandpass signal with bandwidth $B$ extends from $f_c-B/2$ to $f_c+B/2$ (where $f_c>B/2$ is the carrier frequency), the spectra of the in-phase and quadrature components extend from $-B/2$ to $B/2$. Since bandwidth is measured at positive frequencies, the baseband bandwidth is $B/2$.


the best pre emphasis for MMSE equalizer

welcome.
As MMSE estimators are unbiased, try reducing Cramér–Rao bound and maybe your MMSE equalizer can achieve a better performance. Note that in general the MMSE estimator is known to be asymptotically efficient.
Note also that you must specify criteria before say something is "optimal". Even so, it is usually not easy to prove the optimality.


"Is there any exact, i.e. numerical, definition for sparsity?" And by numerical, I understand both computable, and practically "usable". My take is that: not yet, as least, there is no consensus, yet there are some worthy  contenders. The first option "count only non zero terms" is precise, but inefficient (sensitive to numerical approximation and noise, and very complex to optimize). The second option "most elements of a signal are zero or close to zero" is rather imprecise, either on "most" and "close  to".
So  "an exact measure of sparsity" remains elusive, without more formal aspects.
One recent attempt to define sparsity performed in Hurley and Rickard, 2009 Comparing Measures of Sparsity, IEEE Transactions on Information Theory.
Their idea is to provide a set of axioms that a good sparsity measure ought to fulfill; for instance, a signal $x$ multiplied by a non zero constant,  $\alpha x$, should have the same sparsity. In other terms, a sparsity measure should be $0$-homogeneous. Funnily, the $\ell_1$ proxy in compressive sensing, or in lasso regression is  $1$-homogeneous. This is indeed the case for every norm or quasi-norm $\ell_p$, even if they tend to the (non-robust) count measure  $\ell_0$ as $p\to 0$.
So they detail their six axioms, performed computations, borrowed from wealth analysis:

Robin Hood (take from the rich, give to the poor reduces sparsity),
Scaling (constant multiplication preserves sparsity),
Rising Tide (adding the same non zero account reduces sparsity),
Cloning (duplicating data preserves sparsity),
Bill Gates (One man getting richer increases sparsity),
Babies (adding zero values increases sparsity)

and probe known measures against them, revealing that the Gini index and some norm or quasi-norm ratios could be good candidates (for the latter, some details are provided in Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed $\ell_1/\ell_2$ Regularization, 2005, IEEE Signal Processing Letters). I sense that this initial work ought to be further developed (stay tune at SPOQ, Smoothed $p$ over $q$ $\ell_p/\ell_q$ quasi-norms/norms ratios). Because for a signal $x$, $0< p\le q$, the norm ratio inequality yields:
$$ 1\le \frac{\ell_p(x)}{\ell_q(x)}\le \ell_0(x)^{1/p-1/q}$$
and tends to $1$ (left-hand side, LHS) when $x$ is sparse, and to the right-hand side (RHS) when not. This work is published in 2020 as SPOQ: smooth ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry (arxiv preprint, SPOQ published paper, IEEE Transactions on Signal Processing, 2020 and Matlab toolbox).
However, a sound measure of sparsity does not tell you whether the transformed data is sufficiently sparse, or not, for your purpose.
Finally, another concept used in compressive sensing is that of the compressibility of signals, where the re-ordered (descending order) coefficient magnitudes $c_{(k)}$ follow a power law $C_\alpha .(k)^{-\alpha}$, and the bigger the $\alpha$, the sharper the decay.


Wavelet transform wavdec returns matrix of the wrong size

this is due to the fact that the default boundary handling is symmetric extension, which requires several coefficients more in order to achieve the perfect reconstruction. Try

dwtmode('per')

to use periodization to set the representation to be non-redundant. This will however be non-redundant only for input signal lengths being integer multiples of 2^J, where J is number of wavedec levels. Otherwise, very strange padding is used in every step of the wavedec.

P.S.: Sorry for the late answer :)


How to put specific frequency in a specific bin in frequency domain?

I ran your code, there were some errors. f_spac variable is not used after it is filled inside the for-loop. What you can do is to define a variable f_delta=1000. This will make sure the 'f' variable is filled properly with spacing of 1000 and its multiples.

N=8;
r=1;
f_d=1e3;
M=4;
m=0:M-1;
fc=0;
%% freq spacing
for j=1:M
   f_spac(j)=(1+m(j)*r)*f_d;      %freq spacing
end
f_delta = 1000;
%% Generate the frequency vector
for i=1:M
    f(:,i)=[fc+((1+m(i)*r)/2)*f_delta; fc-((1+m(i)*r)/2)*f_delta];
end
B = reshape(f,[1,N]);
fs=20*f(N/2); % Sampling frequency
Ts=1/fs; % Sampling time
t=0:Ts:1e3*Ts; % Time vector

At the bare minimum, you need FFT size to be 8 with Fs = 4kHz. But before that you need to down-convert the signal to baseband.


how important is the pre-emphasis stage in MFCC?

Pre-emphasis is a way of compensating for the rapid decaying spectrum of speech. The experiment is worth trying on real data - you will find that the DCT basis is better at extracting a set of decorrelated coefficients when the spectrum has been whitened by the pre-emphasis filter.

This justification does not really apply in the case of music where the spectrum is more slowly decaying - it is not uncommon in MIR papers to see MFCC computed without the pre-emphasis filter. In the case of music, an argument for pre-emphasis is that it provides a very crude approximation of the outer ear filtering. Some feature extraction modules are indeed using MFCC-like features with a proper outer ear filter (response given by Terhard's formula) in place of pre-emphasis - For example, it looks like this is what is used by the Echo Nest analyzer.


How can we calculate the frequency response of each filter of a Discrete Wavelet Decomposition filter bank?

The decision under the link you provided seems good enough. The main idea is:

we have impulse response of each filter
take the impulse response of the filter with the slowest sample rate as the input signal
route it back through the chain
perform spectrum estimation of the resulting signal

Division by $\sqrt2$ is for normalization only. For the given set of db coefficients their sum is $\sqrt2$, look here for explanation. It has no impact on the sense of the result I think.


Overlapping in real time fourier transform?

Overlapping FFTs, and averaging them, provide an excellent way to more accurately render the noise floor of the signal being observed. That said, you cannot actually start the FFT until you have enough samples, in part because the "last" sample will be used even on the very first butterfly computation.

That said, the time to compute the FFT can be quite short, and shortened further by putting part of it in an FPGA, depending on the signal that can achieve a very short time between the computation and the ability to act on that computation. The latency of that action would then be the number of bins in your FFT multiplied by the sample rate (time to get one FFT) plus the time to execute your recognition code and response. You can compare that time to the rate at which new signals appear or disappear in the sampled signal to see how close you are to "real time."

For audio data recording at 192K samples per second, each sample takes about 5 microseconds to come in, and even a small FPGA can compute a 1024 bin FFT in under a microsecond, so well before the sample after the last one has arrived.


Ideal sampling of audio signals-Periodic spectrum

Discritizing a signal in the ideal sense is equivalent to multiplying by a Dirac delta impulse train.  The Fourier transform of an impulse train is also an impulse train.  Lastly, a multiplication in the time domain is equivalent to a convolution in the frequency domain.  So when we discritize a signal, we are convolving that signal’s frequency transform with an impulse train, resulting in a frequency spectrum that repeats at intervals equal to the sample rate, making it periodic.


Measures that can be used to truncate linear series of functions

We are in the realm of non-linear data approximation. Among the standard litterature, there are Nonlinear approximation, 1998 (R. de Vore), or the more visual On Fourier and Wavelets: Representation, Approximation and Compression, 2006, Vetterli.

Best $N$-term approximation, with respect to some norm $\|f-f_N \|_p$, is quite involved (generically, NP-hard). A basic (suboptimal) approach keeps the $N$ biggest terms (esp. with an orthogonal basis)), which (since $|c_k|\ge|c_{k+1}|$) amounts to finding a proper threshold $\gamma$, and perfom hard-thresholding, up to $|c_N|\ge \gamma >|c_{N+1}|$. There are many strategies to pick an optimal threshold, with additional assumptions on the data, on the noise.

Since you are interested in the Akaike Information Criterion (AIC), you can check A Crossvalidatory AIC for Hard Wavelet Thresholding in Spatially Adaptive Function Estimation, 1998, that compares with other standard wavelet shrinkages.

In general, you have to balance a number of terms (on even bits) with the accuracy. The $\ell_0$ constraint can be converted into other norm penalty. If you want to compress, denoise, you can use other sparsity measures like entropy-based method. Further, you have all sort of matching techniques. You can have a quick overview, for images, in Chapter 4: Redundancy and adaptivity.


U-nets : how exaclty is upsampling performed

The original U-Net paper used nearest neighbor interpolation as far as I know. This is also the default upsampling method in TensorFlow. My own anecdotal advice is to not use transposed convolutions in U-Net. It will only make your CNN slower and won't really increase your F1 score (or whatever metric you're using). I would recommend nearest neighbor interpolation.
Transposed convolution is basically regular convolution with zeros inserted between each input. In U-net after each upsampling, there will be a regular convolution layer. So you don't lose anything if you choose fixed weights like nearest neighbor.
Upsampling outside the neural network domain applies two steps: 1. upsample, 2. lowpass filter. The lowpass filter is here the learnable convolution. In CNN there only seems to be an issue with aliasing (see Making Convolutional Networks Shift-Invariant Again). But this happens also with max pooling / downsampling.


the dominant frequency is overwritten in every iteration using matlab code

You are searching for the maximum on the entire epoched/windowed data-set:

[v,k] = max(datax);

getting the value and index to this maximum. However you should restrict your search on the frequencies of interest:

valid_idx=(freq>=0.1 & freq<=2);
[v,k] = max(datax(valid_idx));

I would leave the find out as you can perform the search on a logical-index variable and thus you do not need to transform the found-index into the original index.

Additionally, I would first demean your data:

data_x = currentwindow_sigx-mean(currentwindow_sigx);

Which would be equal to detrend(currentwindow_sigx,'constant'), without the parameter constant you are detrending and not demeaning. After this I would apply the hamming window. And then calculate the fft.
You should inspect your data afterwards and check if the data looks "physiological", especially in the frequency band of interest. Here you need some experience to judge this. If the data looks okay the maximum search should give you the correct result.

I assume that the fft will be noise, you could apply some smoothing to it before you search for the dominat frequency to avoid outliers.


Why is it called continuous-time frequency?

$\Omega$ is the usual angular frequency in radians per second, and is equal to $2 \pi f$. It is the way to measure frequency for continuous-time signals.

In discrete-time, frequency is measured in radians per sample, and is denoted as $\omega$. Here, a frequency component with $\omega = 2\pi$ is indistinguishable from $\omega = 0$.

$\Omega$ and $\omega$ is the notation in texts such as Oppenheim-Schafer. Other authors use $\theta$ for discrete-time, and the classic $\omega$ for continuous-time signals.


I don't actually have an answer, just ideas, but since this is a riddle, I might of contribute these as (hidden as spoiler text) hints for others, which might have ideas that I miss:

So, MSK is often introduced as

    represented by a half-symbol time-offset BPSK on the I and Q components. That comes in very handy here (avoiding multiplicators) – BPSK can be implemened as $\pm 1$ on the respective axis. The idea is that at the boundary between symbols periods, the phase is continuous, which w.l.o.g cannot be realized in any other way than by having zeros exactly at the point where the other component carries the BPSK constellation point. Thus, the minimal viable MSK modulation has

   $$\begin{align} I&= [- {s_1} && 0 && -{s_2} && 0 && -{s_3}&\dots ]\\Q&= [ 0 && s_1 && 0 && s_2 && 0 & \dots ] \end{align}$$

   as samples, with $s_n$ being the bits to be transmitted represented as $\pm 1$ differentially. That doesn't look very much like CPFSK at all, but one has to realize that the advance between two consecutive samples takes two different values:

   For samples that have an even index $n$, the following sample always is -90° further, no matter which value $s_{n/2}$ has

   For samples that have an odd index $n$, the phase advance is -90° if $s_{\lfloor n/2\rfloor}=s_{(n+1)/2}$, and +90° in case the bits differ.

   Considering the phase difference of consecutive samples as the instantaneous frequency, we thus see that we get an alternation of a constant $-\frac \pi2$ between samples and either $0$ or $-\frac \pi2$, carrying the information. I'm not even convinced any more this fulfills the criteria for being MSK.

   Since we'd probably want spectral symmetry, we'd employ the usual tricks to shift the signal in frequency by half the sampling rate (i.e. "multiplication" with $[-1\,+1\,-1\,+1\,\dots]$).

Oooh and I'd have an approach for the Gaussian pulse shaper, but that simply applies

   the fact that repeated convolution of something with itself converges against a Gaussian – same reason the CLT works.

 That means I can just cascade a set of moving average filters – all of which work without multipliers to get an approximated Gauss filter.

Sadly, that doesn't inherently solve the issue

   that my MSK approach above basically only allows 1 sample per symbol – but we might at least recreate the inter-symbol dependence characteristic of GMSK trellises.


Converting from stereo to mono by averaging

There are indeed rare instances of recordings in which the L&R channels are in opposite phase. This was done in the late 70s and early 80s to create a cheap "spatial sound" effect (like the "surround" button on low-end hi-fi systems of the time), and on such recordings, computing the average of the left and right channels gives a null output. Some stereo effects (like chorus or reverb) would also make a recording sound less lively or flatter after a mono conversion, but this doesn't really qualify as noise. All these are rare problems, and averaging the L&R channels is the standard way of converting from stereo to mono. In the past I've worked on a system that handled large catalogues of audio files, and my approach was to check whether the energy of the average of the channels was below a certain fraction of the energy of the individual channels ; and if so, used the difference instead of the sum as a mono conversion. Ideally you could compute the cross-correlation between the left and right channels to find the lag at which they are maximally correlated, and delay one channel by the lag before the mixing so they are maximally in-phase with each other, but this is overkill for just handling the 15 or so tracks out of 100,000 that had a massive phase inversion problem - to give you an idea of the rarity of the problem...

The most obvious error is converting from 44kHz to 22kHz by skipping every second sample. This will alias all frequencies in the 11kHz - 22kHz range, and will make the downsampled recording sound "metallic", "crunchy" or "gritty". A proper way of downsampling a recording is to first apply a low-pass filter to remove all frequencies above half the target sample rate ; then decimate (skip every second sample).


We want to do inverse DFT:
$$
x(k) = \frac{1}{\sqrt{N}}\sum_n y(n) e^{j 2\pi kn/N},
$$
but can only do forward DFT. So let's try to make it "look" like a forward DFT. Note that we want a minus sign in the complex exponential. We can introduce it via complex conjugation:
$$
\overline{x(k)} = \frac{1}{\sqrt{N}}\sum_n \overline{y(n)} e^{-j 2\pi kn/N},
$$
That is starting to look like a forward DFT now! So the recipe is:

Complex conjugate the given sequence that we want to inverse DFT
Calculate its forward DFT
Calculate complex conjugate of the result.

That gives you the inverse DFT of the original sequence.

(Note that I cheated a little bit on the scaling factor to make it easier!)


Simulating spatial perception - with cosine wave based sound - in MATLAB

Try this

DELAY = 0.5e-3; % half a millisecond for starters.
t = [ 0 : 1 : 2000 ]';   % Time Samples, COLUMN vector
f = 500;    % Input Signal Frequency 500Hz
fs = 44100; % Sampling Frequency
x =  cos(2 * pi * f / fs * t);  % Generate Sine Wave
y =  cos(2 * pi * f * ( t/fs - DELAY));  % Generate Sine Wave
Y = [x,y];

A bunch of problems with your code

Your variable t was a row vector so you generated a mono output signal
The units of your DELAY variable are unclear. I have written it in a way so that it is in actual seconds.


For the sine wave $\sin ((\pi/4)t)$ , the minimum sampling frequency is greater than 1/4 Hz. That amounts to a sampling period less than 4 seconds. In your case, sampling period of 1 second is therefore a valid choice.


Apply band pass filter (BPF) before calculating the power of a signal (PSD)?

You can switch to the frequency domain with regards to Parseval's theorem. I guess you can transform your problem via discrete-time fourier transformation.


Phase information of a signal when creating a spectrogram

Yes, you are correct. If you take a DFT of square wave and only look at the amplitudes, doing an inverse DFT but using different or random phases for the sine components, it does not look like square wave in time domain any more. But it will have a matching spectra. Kind of like two racecars on a circular track, they might always have same velocity, but without knowing the initial position, or even the direction they are going, it will be impossible to know if they side by side or on the opposite side of track.


Why is the high-pass filter result in a discrete wavelet transform (DWT) downsampled?

Assume the original sampling rate is $F_s$ and that the filters are perfect brickwall Halfband filters. After the low pass filtering the frequency content is from $0 \rightarrow \frac{F_s}{4}$. Because of the new bandwidth, the sampling rate can be reduced to $F_s/2$.

For the high-pass filter, the frequency content lies in the range $\frac{F_s}{4} \rightarrow \frac{F_s}{2}$. There is no content in the range $0 \rightarrow \frac{F_s}{4}$. We can downsample by a factor of 2 in this case as well. The downsampling causes the frequency content in the $\frac{F_s}{4} \rightarrow \frac{F_s}{2}$ to be aliased into the range $0 \rightarrow \frac{F_s}{4}$, but because there is no frequency content in this range, due to the high-pass filter, we have not lost any information.


Low-pass elliptic filter design: equation for calculating required minimum order

The estimate is just that: an estimate. If the estimate yields $N=2.4$ then you're pretty safe with a $3^{\rm rd}$ order filter. If the estimate yields $N=2.95$, then a 3rd order might work, but there is a high likelihood that a $4^{\rm th}$ order is required. In cases where the estimate is just under a whole number, the way to be safe is to either increase the estimated order by 1, or to design both a $N$ and $N+1$ order filter to see if both meet your requirements.


It looks like a straight line is not a great approximation in either case. One thing which is unclear to me is that the y-axis in the top plot is $\sigma^2/v$ whereas in the second plot it is just $\sigma^2$ (as far as I can see), and not $\sigma^2/v^2$ as I expected.

What you need to figure out is what you should actually expect. When you say that $\sigma^2$ should depend on $v^2$ then this doesn't say much. From what you do it looks like that you expect those two quantities to be proportional:

$$\sigma^2=\alpha v^2$$

because then you would ideally have a constant ratio $\sigma^2/v^2$. Is this the case?

Note that if you model $\sigma^2/v$ as a linear function of $v$, then you have

$$\sigma^2/v=mv+q$$

which is equivalent to

$$\sigma^2=mv^2+qv$$

So you actually model $\sigma^2$ by a quadratic function of the variable $v$.

To determine the best fit, it seems natural to compare the normalized approximation errors of the different fits. Given that $d_i$ are your measured data points, and $v_i$ are the discrete velocities corresponding to the data $d_i$, the approximation error can be defined as

$$\epsilon=\frac{\sum_i \left[d_i-f(v_i)\right]^2}{\sum_id_i^2}$$

where $f(v)$ is the approximation function, which could be linear as you've suggested in your question, i.e. $f(v)=mv+q$.


Clarification for "Digital Signal Processing" by Steven W. Smith related to complex representation of sinusoids

Phasors are useful for the analysis of (real-valued) linear time-invariant (LTI) systems. A phasor is a complex number
$$C=|C|e^{j\phi}=A+jB\tag{1}$$
which represents a sinusoidal signal:
$$x(t)=\textrm{Re}\big\{Ce^{j\omega t}\big\}=A\cos(\omega t)-B\sin(\omega t)=|C|\cos(\omega t+\phi)\tag{2}$$
That's already the end of the story. Eq. $(2)$ is the generally accepted definition of a phasor. The definition in the book is either a mistake or an idiosyncrasy making things unnecessarily confusing.


Plotting Hilbert Spectrum from Hilbert Huang Transform

The formula for the Hilbert spectrum in your question does not make much sense. Its right-hand side it is not a function of frequency, instead it is simply the sum of the intrinsic mode functions (IMFs), i.e., it represents the original signal. The Hilbert spectrum is the sum of the individual Hilbert spectra of each IMF:

$$H_j(\omega,t)=\begin{cases}a_j(t),&\omega=\omega_j(t)\\0,&\text{otherwise}\end{cases}$$

where $a_j(t)$ and $\omega_j(t)$ are the instantaneous amplitude and the instantaneous frequency (IF) of the $j^{\,th}$IMF, respectively.

Apart from that, note that your computation of the instantaneous phase and frequency will not work. For the phase you have to use arctan2 instead of arctan, and for computing the IF from the instantaneous phase, you need to unwrap the phase. However, it is generally better to compute the IF directly from the signal, as discussed in this answer.


16 bit PCM Audio Sounds Like Bit Crushing

In case anyone has the same problem, I want to tell that the above code works really well. My problem about the distortion of audio caused by wrong buffer size.

But while searching for a solution, I figured out that everyone says that using union may cause undefined behaviour in different CPU structures. It is said to use bitshift operations instead. So I decided to change the code for conversion as below:

Input Conversion

int16_t sample = 0;
    for (int k = 0; k < BYTE_RESOLUTION_IN; k++) {
    sample |= input_buffer[i + k] << (k * 8);
}

Output Conversion

int16_t sample_int16 = (int16_t) temp;
    for (int k = 0; k < BYTE_RESOLUTION_OUT; k++) {
    output_buffer[i + k] = (sample_int16 >> (k * 8)) & 0xFF;
}

And this conversion is made for Little-Endian system. If you have Big-Endian system, you have to change the k value properly.

Have a nice day


My doubt is, when the sound signal is reflected back, it is greatly attenuated and distorted due to many factors.

That's true!

  Do they use MFCC(which I think is only good for audible range)

I don't see how MFCC come into play here, at all. Radar/sonar applications typically depend on some kind of comparison between TX and RX, and a power spectrum only observation can't help you at all here.

  for identifying the echo or are there any other better techniques?

yes, there are – for example, simple cross-correlation would do.


How do optical illusions work (photo mosaic) from a signal processing perspective?

What are the signal processing principles behind our eye merging the many tiny images into a large coherent one?

From what I've read (Source 1: Front-End Vision and Multi-Scale Image Analysis, Source 2: Information Visualization), our eyes apply something very similar to a set of gabor filters to the visual inputs we get. So if you, i.e. focus on Mona Lisa's eye, your brain processes gabor filters at the scale roughly the size of the eye, and such a filter removes details at smaller scales.

  Are there etiological explanations of optical illusions in general - by this I mean can it be explained why image inputs get parsed in 'wrong' ways by the eye/brain, why does the human parser get confused; what transfer function is in play?

Regarding optical illusions in general: There is tons of research. The two books I've mentioned above quote some of it, but I don't know enough about it to compress it into a good answer here.

But regarding this special case, there's actually a very elegant explanation in Front-End Vision and Multi-Scale Image Analysis, p208, which I'll try to repeat here: Take a "natural" sample image. Split the image into small blocks. Multiply each block with a window function (e.g. a gaussian).

For example, I'll take the "Lena" image in Mathematica:

sample = ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"];

split it into 16x16 pixel blocks:

blocks = Flatten[Partition[ImageData[sample], {16, 16}, {5, 5}], 1];

multiply each block with a gaussian window:

window = Table[Exp[-(x^2 + y^2)/25], {x, -7.5, 7.5}, {y, -7.5, 7.5}];
blockTimesWindow = Map[window * # &, blocks];

the blocks look something like this:

Next I'll use SVD to find the principal components of these blocks:

{u, s, v} = SingularValueDecomposition[Flatten /@ blockTimesWindow, 10];
Image[0.5 + Partition[#, 16], ImageSize -> 100] & /@ Transpose[v]

and visualize them:

These are the principal components of small parts of my test image. Or in other words: any small part of my test image is (approximately) a weighted sum of these principal components, plus noise. And it's not just this test image. Take any photograph of a natural scene, the result will look very similar.

As you can see, the principal components look suspiciously like gabor filter kernels or gaussian derivatives. That explains why filters like these are a good choice to analyze natural images: in a natural image, any small part is (almost) the weighted sum of a few gabor filter kernels plus noise. And it explains why evolution would choose this kind of low-level processing.

Now your image is not a natural image: Neighboring pixels are not correlated the way they would be in a natural image. But our eyes still use processes that were optimized for natural scenes, where differences at smaller scales don't carry as much important information, so they are filtered out. That explains why the illusion works.


How to efficiently compute a sliding FFT/DFT with a variable window?

When using the sliding DFT instead of an FFT, you can actually have a different window length per bin. This is because each bin is computed separately anyway. The only constraint is that each window length must be a multiple of its bin wavelength (or half of its bin wavelength if you modify the transform a little). For example, you can do a constant-Q sliding DFT by using the same factor for each bin, and spacing the bins appropriately. More complex variants are possible but require a lot of tweaking, depending on what you want to do with the result. The different windows should be aligned at their center.

I once found a paper about this (after coming up with essentially the same thing myself); it's called "Sliding with a constant Q" and can be googled. For resynthesis, alternatively to the method described in that paper, the trivial resynthesis algorithm for the sliding DFT (essentially just summing everything) also works for variable window lengths.


Phase shift and phase spectrum terms in multidimensional signal

For a real signal, the magnitudes of the Fourier transform bins tell what are the amplitudes of the sinusoidal components that an infinite tiling of the signal frame can be decomposed into. The phases of the bins encode translation of the sinusoids along their direction, like so:

The translation can be decomposed into a series of component translations one along each of the axes. The axes of an image are spatial axes, and in video the third axis is temporal. The component translations have magnitudes inversely proportional to the frequencies indicated by the multidimensional index of the bin. For example, you can easily spatially translate an image in the frequency domain by adding to all phases (in a polar form representation of the bins) a function of form $af_x + bf_y$ where $a$ and $b$ are constants that give the amount of translation and $f_x$ and $f_y$ are the vertical and horizontal frequencies (bin indexes).

Translation of a sinusoid of zero frequency is a bit of an oxymoron. With real FFT its phase is fixed to zero.


Symbol timing synchronization using a high sampling rate

The primary issue with this approach is that it will appear to work fine in your design if you do not consider noise, but in actual application will fail under low or even moderate SNR conditions unless you start doing what an appropriate tracking loop response will do by design and that is average out the sample to sample noise to provide your best estimate of the correct timing location; that which you would not be able to properly estimate on any one given sample. So yes as long as you do the proper filtering/averaging and timing estimation and correct your timing estimate based on current and past updates you will be able to do your timing recovery- but this is what the timing recovery loop is doing and very efficiently.  For any such design that you pursue for timing (and carrier recovery, and data demodulation) be sure to test and benchmark your approach using corner case SNR conditions AND worst case symbol clock offset between the transmitter clock, signal as received, and receiver clock!

To add to this since not detailed in my other links is the Gardner Timing Error Detector (TED) implementation suitable for QPSK and QAM implementations, where the symbol $y_n$ shown is the complex I,Q output. I have not implemented this for OQPSK specifically, but do not see why simply offsetting the I and Q samples half a symbol would not work. The TED produces an average level that is proportional to timing error and will have significant pattern noise on any one given output so is assumed to be averaged as part of the overall loop design:

Note the the Gardner TED has a relatively high capture range for carrier offsets (see Location of Matched Filter for more details on that), so assuming the estimation of the carrier is within this range, the timing recovery can be done ahead of carrier recovery as shown in this QAM receiver implementation (two samples per symbol for timing recovery, which once acquired is decimated to one sample per symbol for carrier recovery. Also depending on the waveform pulse shaping used and timing error determination the inputs to the timing recovery may be better before or after matched filtering (hence the dashed lines in the block diagram showing that either may be used).

Below is a graphic that compares an approach of changing the frequency of the sampling clock to the ADC in the receiver vs interpolation and resampling. The latter is more prevalent in modern software radio solutions.

(It is unfortunate that the two implementations are in reverse direction flow and I apologize for that; it is the figures I had readily available)

The following graphic demonstrates the concept of interpolation and resampling :

A polyphase filter approach as publicized by fred harris would do this interpolation and retiming efficiently, so that you never actually need to interpolate to a higher rate than 2 samples per symbol. Note that you also only need three filter implementations in the poly-phase filter bank: a present and the two adjacent. After changing from present to adjacent, the filter coefficients can be updated from a ROM table. (The reason for keeping the three is so that all filters needed are properly loaded with the waveform history prior to switching two its output).

For further details on polyphase retiming and implementation of polyphase filters see these links:

Emulating a Variable Delay

How to implement Polyphase filter?


Choosing a Sampling Rate and a Cutoff frequency

The question states

  … generate … digitally … same sample rate to generate a 20 kHz tone.

So, you want to digitally create a 20 kHz tone, too.  You hence need to have a sample rate that is sufficient for 20 kHz.


How to apply an anti-aliasing filter before downsampling

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


You tell the phase demodulator that the signal has a maximum phase deviation of $\pi/2$, but then you actually use a signal with a phase deviation of $5\pi/2$ (due to the amplitude of the message signal). This doesn't make sense.

Note that even the most ideal phase demodulator cannot distinguish between $\cos(2\pi f_ct+\phi(t))$ and $\cos(2\pi f_ct + \phi(t)+2k\pi)$, $k\in\mathbb{Z}$.


How can we show that $|H(e^{j\omega})|=|H(e^{-j\omega})|$?

Your statement only holds true for a real sequence $h[n]$.
The frequency response of $h[n]$ equals to
$$H(e^{j\omega}) = \sum_{n=-\infty}^{\infty}h[n] e^{-j\omega n}$$
and if $h(n)$ is real, the conjugate symmetry condition holds
$$H(e^{-j\omega}) = \sum_{n=-\infty}^{\infty}h[n] e^{j\omega n}=\big(H(e^{j\omega})\big)^*$$
So we can derive that
$$\big|H(e^{j\omega})\big| = \big|H(e^{-j\omega})\big|$$
and
$$\arg\{H(e^{j\omega})\} = -\arg\{H(e^{-j\omega})\}$$
where $\arg\{\cdot\}$ represents phase response.
We can also derive that the real part of the frequency response is even symmetric and the imaginary part is odd symmetric.
$$\Re\{H(e^{j\omega})\} = \big|H(e^{j\omega})\big| \cos\big(\arg\{H(e^{j\omega})\}\big) = \Re\{H(e^{-j\omega})\}$$
$$\Im\{H(e^{j\omega})\} = \big|H(e^{j\omega})\big| \sin\big(\arg\{H(e^{j\omega})\}\big) = -\Im\{H(e^{-j\omega})\}$$


I have the same book as you. :-)

The [1 -1] filter is a a simple differentiator, also known as, the Haar mother wavelet. You need to understand convolution to understand what he is saying. On the first level, you convolve your signal with [1 -1]. Then on the next level, you convole your signal with [1 0 -1], etc. That is all the CWT is doing. (Look at convolution first). If you understand that well, you can understand CWT very easily.


The given system is not time invariant due to the term $\cos(\pi t / 3)$.

Given a system definition of the form:

$$y(t) = 3 x(t) + \cos( \pi t/3) $$ it's easy to see that

$$ y(t-d) = 3 x(t-d) + \cos(\pi (t-d) /3) \neq T\{x(t-d)\} = 3 x(t-d) + \cos(\pi t/3)$$

hence the system is time varying...


You can solve this using scipy.fftpack (sfft) instead of np.fft, because the sfft implementation can be directly used on 2-dimensionnal arrays so you don't have to do it in a convoluted way (ba dum tsss).

In your code header, add :

import scipy.fftpack as sfft

Then, to calculate the fft and shift the spectrum, use :

Z_fft = sfft.fft2(Z)
Z_shift = sfft.fftshift(Z_fft)

The obtained spectrum is then nicely arranged for image display :

plt.figure(4)
plt.imshow(np.abs(Z_shift))

Also, the way you are constructing the circle seems overly complicated, you can take advantage of python's syntax using boolean syntax :

X,Y = np.meshgrid(x,y)
Z = X**2+Y**2
R = 1
Z[Z > R] = 0;Z[Z > 0] = 1;

What I did here is using the meshgrid you created to calculate the circle equation result, in the form of Z. Then I test Z for all values above R, put them to zero, and the other values to 1.

My code is below :

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.fftpack as sfft # you have to import this package

plt.close('all')

numb_steps = 100
x = np.linspace(-10, 10, numb_steps)
y = np.linspace(-10, 10, numb_steps)

#constructing the disk in a simpler manner
X,Y = np.meshgrid(x,y)
Z = X**2+Y**2
R = 1
Z[Z > R] = 0;Z[Z > 0] = 1;

# Here, use sfft routines instead of np.fft
Z_fft = sfft.fft2(Z)
Z_shift = sfft.fftshift(Z_fft)
FreqCompRows = np.fft.fftfreq(Z.shape[0],d=2)
FreqCompCols = np.fft.fftfreq(Z.shape[1],d=2)
FreqCompRows = np.fft.fftshift(FreqCompRows)
FreqCompCols = np.fft.fftshift(FreqCompCols)

S,D = np.meshgrid(FreqCompRows, FreqCompCols)

# plotting your disk signal image
fig = plt.figure(1)
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y,Z)

# plotting "3D" fft non-shifted
fig = plt.figure(2)
ax = fig.gca(projection='3d')
ax.plot_surface(S, D, np.abs(Z_fft))

# plotting images of FFT's
plt.figure(3)
plt.imshow(np.abs(Z_fft))

plt.figure(4)
plt.imshow(np.abs(Z_shift))


Spectrogram with square or non-square magnitude of STFT: power vs. magnitude

In the image you gave, the right-hand scale is in decibels (dB). So, essentially, the square turns into an affine scaling in the logarithm domain, which essentially yields the same image, at least the same relative dynamic range. So, in that case, it does not matter.

Aside, $x \to x^\alpha$ transformations with $x\ge 0$ (a positive signal, a magnitude spectrum, a spectrogram) are customary when signals have different orders of magnitude that cannot be caught by human eyes or ears: $\alpha > 1$ puts emphasis on the largest components,   $\alpha < 1$ relatively enhances smaller values, which is useful when you want to detect low-amplitude signals.
Depending on whether you want to capture audio fingerprints, or perform subtle audio forensics, one or the other should be chosen. Such $\alpha$-power transformation are also useful in your other spectral subtraction question, What is the name of this very simple spectral subtraction technique?


How to Map CWT to Synchrosqueezed wavelet transform?

Let me explain the intuition briefly.  The authors of the paper you've cited assume that the signal $x(t)$ can be written in the form
\begin{align*}
x(t) &= \sum_{k=1}^K a_k(t) \exp(2\pi\mathrm{i} \phi_k(t)),
\end{align*}
where the $a_k$ denote instantaneous amplitudes, the $\phi_k$ denote instantaneous phases and the $\phi'_k$ denote instantaneous frequencies. The continuous wavelet transform enables the visualization of the $\phi'_k$ in the time-scale plane, but the visualization is necessarily blurry due to the Fourier uncertainty principle.  The idea is that the peaks of the continuous wavelet transform can give you an idea of what the $\phi'_k$ values are, but the energy is spread out in the time-scale plane, so the exact values of $\phi'_k$ may not be easy to obtain.  The goal of the synchrosqueezing transform (SST) is to partially undo the blurring by a frequency reassignment, which can be explained as follows.

First of all, let's note that each scale $s$ corresponds to a natural frequency $\xi$ satisfying the relation $s = c/\xi$, where $c$ is the center frequency of the mother wavelet $\psi$.  Now, suppose the time $u$ is fixed.  If $\xi = c/s$ is close but not exactly equal to an instantaneous frequency $\phi'_k(u)$, then the coefficient $W(u,s)$ will have some nonzero energy (i.e., $|W(u,s)|^2 > 0$).  The idea of synchrosqueezing is to move all this energy away from the frequency $\xi$, and reassign its frequency location closer to the instantaneous frequency $\phi'_k(u)$.  Then, we can arrive at a time-frequency representation where the energy is more closely concentrated around the instantaneous frequency curves.

So, for each scale $s$, we compute the frequency reassignment  $\omega(u,s)$ of the wavelet transform coefficient $W(u,s)$ by the formula you cited above.  As Daubechies, Wu, and Lu discovered in the paper you've linked, it turns out that $\omega(u,s)$ is often a good approximation to the instantaneous frequency curve $\phi'_k(u)$ for values of $(u,s)$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$. (In fact, when $\phi'_k(u)$ is a constant function, $\omega(u,s)$ exactly equals $\phi'_k(u)$ for values of $s$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$!)

Then, the computation of the SST in the continuous setting is as follows. First, fix a time of interest $u$. Next, compute the frequency reassignment $\omega(u,s)$ for all scale values $s$. Then, for each frequency of interest $\eta$, we compute the SST $T(u,\eta)$ by adding up all values $W(u,s)$ where the reassigned frequency $\omega(u,s)$ equals $\eta$. This means

\begin{align*}
T(u,\eta) &= \int_\mathbb{R} W(u,s) \delta(\eta - \omega(u,s)) s^{-3/2} ds,
\end{align*}

where $\delta$ is the Dirac delta. (To ensure the convergence result for SST, the authors use a formulation for $T$ which includes an approximation to $\delta$ rather than $\delta$ itself, but in computational practice I've never found this approximation to be necessary.) The multiplication by $s^{-3/2}$ is necessary for reconstruction purposes (which you can read more about in the paper).

In practice, one is limited to the discrete setting where we have only finitely many possible frequency bins $\eta_\ell$. In this case, for a fixed frequency $\eta_{\ell_0}$, one finds all the values of $\omega(u, s)$ which are closer to $\eta_{\ell_0}$ than to any other frequency bin $\eta_\ell$. We can more explicitly formulate the SST in the discrete setting as follows.

For simplicity's sake let's assume the possible frequency values $\eta_\ell$ are uniformly spaced by a distance $\Delta \omega$. So, in the discrete setting, the formula above becomes
\begin{align*}
T(u, \eta_\ell) &= \sum_{s: |\omega(u,s) - \eta_\ell| < \Delta\omega/2} W(u,s) s^{-3/2},
\end{align*}
where I haven't bothered to give an explicit notation that shows that the values $u$ and $s$ are discrete.

By doing this reassignment, you end up with a representation $T$ which is sparser than $W$, and hopefully very sharply concentrated about the curves $\phi'_k$. (In practice, for visualization purposes, it's better to reassign not the original coefficient $W(u,s)$ but the magnitude-squared $|W(u,s)|^2$, or even just the constant number $1$. This is because the coefficients $W(u,s)$ are complex-valued, and summing these coefficients may not actually leave you with a very large energy. However, if you additionally want to use the reconstruction formula given in the paper to compute the $x_k$ from the SST, then you need to sum the $W(u,s)$ in order for that formula to work.)

By the way, I am leaving out the fact that we generally throw out coefficients $W(u,s)$ that do not fall above a certain threshold.  This is because the computation of $\omega(u,s)$ is not necessarily accurate for small $W(u,s)$.


Frequency response of a rolling linear regression

This is a statistics question and should not be asked here. But I will try to answer the question as best as I can (it won't be the best answer). First off, I don't see what the phenomenon you are observing has to do with frequency response or LPFs. Low pass filters makes data smoother, which is not what you have. Instead you have a situation, where estimate for one set of data differs greatly from a other set of data. The OLS estimator also is not an MA process, the formula being cov(x,y)/var(x). The filter described in the paper is not a linear time invariant filter.

Indeed, if you regress a line on white noise process or linearly growing process, the regression line should not change slope significantly as more data is gathered. Instead as the amount of data grows you should be more certain of that the coefficient is right, generating greater variance for the rolling process, but on average the slope should be right.

What you have here is not a process that satisfies Gauss Markov assumptions. As can easily be seen the expected value of the error is not zero, and there is heavy auto correlation. This means that the error term is correlated with the independent variable, producing a biased estimate. Besides that the data is not randomly selected.

Most importantly, the data is increasing as a function of time. This generates a so called spurious correlation. Since Y increases faster than X as time passes, the beta becomes upwardly biased. You can model autocorrelation with an MA process, but this is not really important...

I am sure statistics stackexchange (crossvalidated) would provide you with much more mathematically precise answer and a better explanation of what the phenomenon. One thing fix that might help is to difference both of the variables, which changes the autocorrelation structure of the error. This can be equally viewed as a (1-1) MA process applied to both of the variables, which is a comb filter at a high frequency, producing a high pass filter below the nyquist frequency. But the data has such a non-standard structure that I wouldn't be so sure that the fix will help you.


Relation between samplingrate and frequency

you have 9 $\times$ 44100 = 396900 samples. I assume you did a single FFT and you should take the first , 198450 complex points for your plot after taking the magnitude square.  The last 198450 points out of the FFT are a mirror image of the first 198450 points so you can neglect them. The first kept  point corresponds to $f=0$ and the last kept point to $f=22050$ Hz.  which makes each bin  .1111 wide


why transmitting a long pulse its spectrum gets spread?

this is John BG

1.- GFSK is not same as single tone

According to the CC1101 specs you refer to, the C1101 uses GFSK.

GFSK varies the varrier frequency, therefore the spectrum cannot be a delta.

Call it a 'shifting delta' but it cannot be same spectrum as just a tone.

2.- you want the carrier frequency to vary

ISM are crowded bands, the carrier has better chances to deliver if varying than if constant,
against jamming, interferrers, and any power level that neighbouring transmitters may radiate.

3.- not even ASK

You can still set modulation to MSK, OOK, or even ASK but even with ASK

Even with PKT_FORMAT[1:0] to Asynchronous mode, CRC disabled, Manchester disabled and PKTCTRL0 setting infinite packet length, PKTLEN is going to cut the packets to ensure syncing takes place.

This is done to avoid persistent propagation of errors along same packet that might generate uncertainty
to possibly overlap next packets.

4.- Temperature drifts, manufacturing tolerances

Such small and really low power devices may be subject to temperature drifts and manufacturing tolerances.

The packet header is needed at a high enough repetition to avoid transmitter and receiver missing frequency lock.

Hope it helps, John BG


What's the point of using Morphological processing concerning images?

Many operations can be performed on images for enhancement and restoration. Sometimes, they seem mere handcraft and tinkering.

To make more sense of such workflows, and provide it with solid roots, it is useful to consider that  images are classes of abstract objects, and operations as well-defined structural actions on the former. Depending on premises or assumptions on the image shape (is it binary, continuous-valued, what do we expect from filtering, what elements are meaningful, etc.), several constructions can be derived.

Standard linear filtering techniques in images assume that images are composed of linear superpositions of components. This yields convolution-based methods, and treating images as elements of a linear vector space (the underlying structure), with a strong "continuous" underlying theory. This can be related to Fourier or harmonic analysis, such a strong influence of least-squares method.

But image formation is often non-linear, for instance because of occlusion or saturation. So you have a collection of non-linear tools, with homomorphic, min, max or median filters for instance, or stack filters for bit-wise versions. Plus,  they can be discretized in both values and space.

Mathematical morphology is another comprehensive framework, based on another underlying mathematical structure, lattices/set theory and related operations (erosion, dilation), with a strong discrete background (especially on image values).

Both have been proved recently connected in some ways, through the Cramer transform, see for instance An Explanation for the Logarithmic Connection between Linear and Morphological System Theory. Another interesting reading is 1989's Fourier analysis, mathematical morphology, and vision by C. Ronse.


Detect keypoints in images without corners

The computer vision system toolbox in Matlab includes Feature Detection, Extraction, and Matching. The line, edges, histogram of oriented gradients are all included other than the corner features. And there is also an example on point feature matching in the toolbox (the names are called corners but they are actually not).


Discrete Fourier Transform of PMF and PDF

A random variable (RV) doesn't need to be continuous in order to have a PDF. The PDF of a discrete RV, i.e. a RV that takes on only discrete values, is of course not a continuous function but is only different from zero at the possible discrete values of the RV.

An example would be the outcome of a fair die. It takes on the values 1 to 6 with equal probability $1/6$. It's PDF is given by

$$p(x)=\frac16\sum_{i=1}^6\delta(x-i)\tag{1}$$

where $\delta(x)$ is the Dirac delta impulse. The Dirac delta impulses are necessary in the PDF (as opposed to the PMF) because probabilities are computed from the PDF by integration. Since for any interval $[a,b]$ containing the value $0$

$$\int_{a}^b\delta(x)dx=1$$

the PDF given by (1) integrates to $1$, as desired. And the probability of the RV $X$ being inside an arbitrary interval $[a,b]$ is

$$P[a\le X\le b]=\int_a^b p(x)dx$$


Automotive Audio Question: IIR Filters and minimum phase

Two minimum phase IIR filters with different magnitude responses generally will have different phase responses. But two minimum phase IIR filters with exactly same magnitude responses will exactly have same phase responses. (How can two systems with exactly same transfer functions have different magnitude/phase responses?). If the first filter is a minimum phase filter and the second filter has exactly same magnitude response but different phase response then the second filter is not minimum phase.

In contrary, two linear phase FIR filters with different magnitude responses may have exactly same phase responses when the required conditions are fulfilled. (E.g same delays, same symmetric types, and no negation).