The model of noise in Image denoising problem

Theory behind it: the additive Gaussian noise  can be related to (Johnson–Nyquist) thermal noise in standard images. For the  denoising of more generic multivariate images, chapter 2.2.2 provides some references on why global noise present in acquired data can be realistically modelled by an additive zero-mean spatially white Gaussian noise.

The shot noise can be modeled as a Poisson process, that behaves, for large numbers, like a Gaussian noise.

Other types of noises are observable in images: salt-and-pepper, quantization (more uniform), grain. In some image modalities, noise can be not stationary, multiplicative. They can be modeled by stochastic  processes in general.

Gaussian noise is also widely used because it enjoys relatively tractable theories, related to least-square estimation. A very practical reason.

Note, however, that most images have 8-bit channells, with integer values in $[0,\ldots,255]$. hence, adding a Gaussian noise yields non-integer negative or values above 255, which is not  fully realistic: natural sampled  noise is  likely to be quantized in the  admissible values.


How can I generate a sine and cosine function using an exponential function?

If you have defined a sampling frequency then you should use it in defining your time grid "nn", otherwise you are just incrementing in digital steps.

The grid "nn" is basically a digital grid rather than representing true sampling based on your defined sampling frequency. Also the range of your time grid should be in mutliples of the time period of sinusoid (it could be anything but if you want to define it in that terms then define it in terms of time period of sinusoid and not the frequencies itslef).I have modified your code below to view 5 cycles of the sinusoid. Hope it helps

fc = 4092; % carrier frequency Hz

fs = 16368; % sampling freq Hz

nn = [0 : (1/fs):5*(1/fc)] % array with numbers from 0 to fc

y1 = exp(i*2*pi*fc*nn); % exponential equation


Why are polynomial filters for gradients in the opposite sign?

If you look at the formula above the table you linked to, you see that the output values are computed like

$$Y^{\prime}_j=\sum_{n=-4}^{4}D_ny_{j+n}$$

In Matlab this could simply be implemented by a vector product (for one output value $Y^{\prime}_j$, where the vector $y$ contains the input values around the current index $j$):

yd = D'*y;

On the other hand, convolution of these sequences is defined by

$$Y^{\prime}_j=\sum_{n=-4}^{4}D_ny_{j-n}$$

(note the '-' sign in the index of $y_j$!). So if you want to use the 'conv()' function in Matlab, you need to flip the coefficients in the vector $D$.


Is it possible to reconstruct an image using its magnitude only?

Re-construction using just the magnitude (and an assumed phase of zero) will only work for exactly symmetric images.  The phase is needed to make the top look different from the bottom and left look different from the right side of the image, etc.

That's because all the cosine waveform DFT basis vectors are exactly symmetric around the center of an FFT or DFT window.  Thus, unshifted, they can only reconstruct a symmetric result.  You could shift the phase of each basis vector randomly (or 90 degrees to make them all sinewaves) but, without knowing the phase shifts of the original image's DFT (or guessing right), you would just get another garbage result.


Why does adding a time delayed version of a signal to itself create a filtered signal?

When you delay a signal by $T$ seconds and add it to the signal itself,
you are cancelling out or nulling the signal component at frequency
$\frac{1}{2T}$ Hz since that signal component will have changed phase
by exactly $\pi$:
$$\begin{align}
\sin\left(2\pi\frac{1}{2T}t + \theta\right) +
\sin\left(2\pi\frac{1}{2T}(t-T) + \theta\right)
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
+ \sin\left(2\pi\frac{1}{2T}t + \theta  - \pi\right)\\
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
+\sin\left(2\pi\frac{1}{2T}t + \theta\right)\cos(\pi)\\
&\ \hspace{0.2in}-\cos\left(2\pi\frac{1}{2T}t + \theta\right)\sin(\pi)\\
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
-\sin\left(2\pi\frac{1}{2T}t + \theta\right)-0\\
&= 0.
\end{align}$$
A similar thing happens at odd multiples of $\frac{1}{2T}$ Hz also.
For nearby frequencies, the cancellation is not as complete, and
of course, at even multiples of $\frac{1}{2T}$ Hz, the signal
component is doubled in value instead of being cancelled.
Similarly, if the delayed signal is reduced in amplitude,
cancellation is not complete at $\frac{1}{2T}$ Hz etc.

To summarize, the signal is being filtered because different
frequencies are being passed through with different gains.

If you want the frequency-domain explanation, the transfer
function $H(f)$ of the system is the Fourier transform of what Matt's
answer gave as the impulse response, viz.
$$\mathcal F\left[\delta(t) + \delta(t-T)\right] = 1+\exp(-j2\pi fT)$$
which is a nonconstant function of $f$ (in fact, $|H(f)|$ varies
sinusoidally from a maximum of $2$ to a minimum of $0$ as discussed
above), and so $Y(f)=H(f)X(f)$ is
not a scalar multiple of $X(f)$. Filtering!


The Hermitian symmetry is used to obtain a real-valued time-domain signal. It is a special case of OFDM called discrete multitone (DMT). It exploits a property of the discrete Fourier transform (DFT), namely that the DFT of a real-valued signal has Hermitian symmetry. The motivation is usually the channel: if the signal shall be transmitted over a low-pass channel without additional radio-frequency upconversion (which requires an oscillator at the transmitter and a local oscillator at the receiver) the transmit signal must be real-valued.

Given an OFDM system with $N$ subcarriers the following conditions must hold
$$
c_{-k} = c^*_k\\
c_0\in \mathbb{R}
$$
where $c_k$ is the complex value of subcarrier with index $k\in[-N/2\ldots N/2-1]$. This can be achieved by mapping $N/2$ complex symbols to subcarriers $0$ to $N/2-1$ and then assigning the respective complex conjugate value to the according subcarrier with negative index. This would normally take place after the mapper and before the inverse DFT (IDFT) block.


Hint: It's a system with a damped oscillation of a single frequency sine wave, which also starts with an arbitrary phase. That should give you amplitude curve (damping), oscillation component (sinusoidal wave) and phase (starting point of the wave). If you can figure out the general form this system takes mathematically, it's a matter of plugging in the right values for those three.

Another hint:
Don't overthink the system order. Think about what a first order system could look like. Do you need to go further? What could a second order system look like? Do you need to go further, etc. I guarantee you will stop pretty early.

Perhaps a more important hint:
You're given a step response, not an impulse response, but there's a way to go from one to the other (what is it). Once you find the expression to describe this curve, going to impulse response should be very easy.

One last hint:
Depending on how your subject is taught, this may help or this may hinder your attempts, so take this one with a grain of salt. If it confuses you, ignore it. My Linear Systems professor said this and made us repeat this out loud in class almost every lecture: Solutions to linear systems are always exponentials!


How do I find the transfer function in the frequency domain?

Your final answer looks good, but the derivation and notation are a bit mixed.  Remember the capitalized $H()$ refers to the frequency domain transfer function, so writing $H(n)$ with the time-domain index $n$ doesn't make sense.  Also, it is the time-domain impulse response that is a delta function, not the frequency domain transfer function, so $H(f)=-\delta(f)$ is not right.

The problem can be solved very easily:

$y(n)=-x(n)$

Take the Fourier Transform of both sides:

$Y(f)=-X(f)$

Divide both sides by $X(f)$:

$\frac{Y(f)}{X(f)} = H(f) = -1$

In the time domain:

$h(n)=-\delta(n)$


How poles are related to frequency response

One thing that really helped me understand poles and zeros is to visualize them as amplitude surfaces.  Several of these plots can be found in A Filter Primer.  Some notes:

It's probably easier to learn the analog S plane first, and after you understand it, then learn how the digital Z plane works.
A zero is a point at which the gain of the transfer function is zero.
A pole is a point at which the gain of the transfer function is infinite.
Often there are zeros or poles at infinity, which aren't always included in descriptions of the transfer function, but are necessary to understand it.
The frequency response in the S plane happens along the jω axis only.

The origin is 0 Hz, or DC, and the cutoff frequency of filters increases radially away from the origin.  Putting a pole at any point along a circle at a certain distance from the origin will produce the same cutoff frequency.
To increase the cutoff frequency of a filter, move the poles radially outward.
To increase the Q of a biquad filter, move the poles along the circle towards the jω axis, which keeps the cutoff frequency constant, but increases the effect that the pole has on the frequency response, making it more "peaky".

If a zero appears on the jω axis, then the frequency response will drop to zero at that frequency; if you input a sine wave at that frequency, the output will be 0.
If a pole appears on the jω axis, then the impulse response is an oscillator; any impulse will cause it to ring forever at that frequency.  Impulses have finite energy, but the response of the filter has infinite energy, so it has infinite gain.

A simple example is an integrator H(s) = 1/s:

This function equals 0 when s is infinite, so it has a zero at infinity.
This function equals infinity when s is zero, so it has a pole at zero.

In other words, it has infinite gain at DC (the step response of an integrator is forever-increasing), and the gain decreases as frequency increases:

Moving the pole away from the origin, along the imaginary axis into the left hand of the S plane, makes the gain at 0 Hz on the jw axis finite again, and now you have a low-pass filter:


What is difference between terms $X(j \omega) ,X(\ e^{j \omega })$ and $ X(\omega)$?

they should have different subscripts.

bilateral Laplace Transform:

$$ \mathscr{L} \left\{ x(t) \right\} \triangleq X_\mathcal{L}(s) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-st} \ dt $$

continuous Fourier Transform:

$$ \mathscr{F} \left\{ x(t) \right\} \triangleq X_\mathcal{F}(\omega) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j\omega t} \ dt $$

it's pretty easy to see that

$$ X_\mathcal{F}(\omega) = X_\mathcal{L}(j \omega) $$

likewise in discrete-time, the bilateral $\mathcal{Z}$ Transform:

$$ \mathcal{Z} \left\{ x[n] \right\} \triangleq X_\mathcal{Z}(z) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] z^{-n} $$

and the Discrete-Time Fourier Transform (DTFT):

$$ \mathcal{DTFT} \left\{ x[n] \right\} \triangleq X_{\small{DTFT}}(\omega) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j\omega n} $$

it's also easy to see that

$$ X_{\small{DTFT}}(\omega) = X_\mathcal{Z}(e^{j\omega}) $$

now often the context is understood and everyone knows that $X(e^{j\omega})$ is the Z transform evaluated at $z=e^{j \omega}$, which is the DTFT and we don't bother with the clunky subscript notation.  we might leave the space for subscripts for another function, like if we have a vector of signals as with a state-variable system representation in control theory.


How does direction finding work when using GMS mobile station and additional handset?

The fact that this thing only has a single antenna port and a loudness knob:

This will just look for a transmitter in the band of interest, and you'll have to turn around a directive antenna until the tone gets loud.

That's absolutely low-tech, it can't discern between GSM phones logged into the same network (and also active), and the whole setup will only work if the straightest path to the phone is also the one that delivers the most energy – which is not the case for e.g. a lot of the indoor and urban mobile channel scenarios that we use nowadays.


Resampling of signal with non uniform sampling frequency

Without specific constraints on the data/noise properties or sampling assumptions, smoothing splines could be helpful. Indeed, constraining the curve to pass exactly through the given points could be too harsh.

One example of such a toolbox in Matlab is SPLINEFIT with several examples:

  Direct spline interpolation of noisy data may result in a curve with
  unwanted oscillations. This is particularly bad if the slope of the
  curve is important.  A better approach is to reduce the degrees of
  freedom for the spline  and use the method of least squares to fit the
  spline to the noisy data.  The deegres of freedom are connected to the
  number of breaks (knots),  so the smoothing effect is controlled by
  the selection of breaks.

Matlab also has Cubic smoothing spline (csaps).

One possible reading:

Flexible Smoothing with B-splines and Penalties, Paul H. C. Eilers; Brian D. Marx, 1996


What does it mean that the normalized power of a received equals one?

In mundane words: the concept is very close to that of orthogonality. For two components of $x$, say  $x_i$ and  $x_j$, if they are different ($i \neq j)$, the average outcome of their product is zero (or $E(x_i x_j^H)  =0$). If   $i = j$, one assumes that the energy is normalized to one (or $E(x_i x_i^H)  =1$).

Resultingly, the product matrix is a unit matrix: one for the diagonal ($i = j$), and zero outside.


Lattice structure realization of an all-pole system

Rearrange into a difference equation form and your answer will become clearer!


What is the type of noise in this image?

Looking at it, my thought is that the top image has been printed out (onto poor quality paper), then later scanned.  The artefacts would then be a mixture of paper texture and poor printing.  It could also include ink transfer from the facing page (if part of a book or newspaper, and/or bleed-through from something printed on the back of the paper.


How can I compute a log-spaced power spectrum?

I found a paper that addresses this question directly:

Michael Tröbs and Gerhard Heinzel, "Improved spectrum estimation
 from digitized time series on a logarithmic frequency axis," in
 Measurement, vol 39 (2006), pp 120-129. (Or free pre-print.)

The first few figures in the paper nicely illustrate the problem that this algorithm solves, and the references contain a useful bibliography of other approaches (constant-Q transform, tempered Fourier transform, a survey article, etc).

Their approach is not to re-bin the output of an existing FFT-based power spectrum estimation, but to only compute the discrete Fourier transform at the (logarithmically-spaced) frequencies of interest.  For each frequency to be estimated, they basically implement Welch's algorithm, but with a transform length (and hence also, number of averages) specifically chosen for each frequency. The computation of each frequency bin uses the entire time series, but segmented differently. The results have the desirable property that the resolution (bin width) is a smooth function of frequency, and the results can be calibrated as either a power spectral density or a power spectrum.

Matlab implementation here: https://github.com/tobin/lpsd

Disclosure: The authors of this paper are at the same institution as me.


Narrowband DOA Beamforming techniques on Wideband signal

My suspicion is that you can only perform phase shift beamforming techniques on narrowband signals and will not work even if you are looking for a narrowband within a wideband signal. Am I correct?

No; the correlation coefficient of your reference signal to your received signal should, given additive, phase-uniform, white noise, have an expectation value of the phase shift the individual receiver sees, no matter how wide your band is.

Note that this is only true if your transmitter/receiver system can be modeled as linear phase for the bandwidth we're talking about, but without loss of generality, radar/sonar/sound systems are designed to pretty much have exactly that characterestic.


The problem is that you can't use a Dirac delta impulse as an input, which would be the theoretically correct thing to do. Since you only have a discrete impulse, which has some finite value, and since Simulink probably uses something like a zero-order hold to make a continuous-time signal out of that impulse you actually use two different input signals: a step and a smoothed step.


How is angle of arrival estimated?

The estimation of angle in linear array antenna is based on the relative phase variation between the target location and antenna elements. This phase variation in spatial dimension between antenna elements will provide you with the location of the target (in terms of angle).

Thus, the basic approach to obtain an angle through array antenna is to perform Fast Fourier Transform (FFT) over the antenna elements. This conventional method with provide you with the angle of target but with very coarse angular resolution (around 2/N radians) (angular resolution using FFT over linear array)

For obtaining the angle with higher resolution, we can use CAPON beamformer or Super resolution algorithms such as MUSIC, ESPIRIT etc.

In order to answer your question,

  Algorithms such as Capon and Bartlett use the position of the source
  for calculate the auto- correlation matrix which is further used to
  calculate pseudo spectrum

This is not entirely true, as the auto correlation matrix in spatial dimension is the input to the CAPON beamformer and till that point the algorithm does not know the location of the target.

The angle of arrival of the target is determined by taking possible set of direction of arrivals and then suppress the other direction which is not under consideration. One has to scan through all possible direction of arrivals (using array manifold matrix) and then with the spectrum estimate the direction of arrival (normally choosing the highest peak).

For your kind reference:

Algorithm for CAPON is mentioned below:

Take auto correlation (A) (among antenna elements)
Take inverse of auto correlation matrix (Ainv)
Create array manifold matrix (M)
Obtain capon spectrum using = 1/ ((hermitian(M) * Ainv * M), where hermitian operation is the conjugate transpose.
Finding the peak spectrum location (provides angle of arrival)

  How is the angle of arrival being estimated if the source position is
  required to run the algorithm in the first place?

So we do not need to know the source dimension, but we have to find it among the possible source position/direction (obtaining the maximum likelihood of the actual direction with possible directions)


Frequency response and impulse response of black box

If the system is continuous-time then you can't simulate an impulse (Dirac delta).

Generally people take a number of approaches:

Instead of using an impulse, find the step response of the system of interest.  This can be related to the integral of the impulse response as the unit step is the integral of the unit impulse.
Generate a set of spaced sinusoids and measure the output phase and amplitude of each.  The transfer function is found by interpolating these discrete frequency points.
Generate a swept sine wave (chirp) signal over a known frequency range and deconvolve the output to find the transfer function.
Generate a pseudorandom noise sequence, measure the response, and use it to estimate the transfer function.


Find the frequency response if i have the magnitude response?

The frequency response of a system can be represented in polar format, in which the magnitude and phase response are considered separately:

$$
H(\omega) = |H(\omega)| \angle H(\omega)
$$

With this representation, it should be clear that the magnitude response alone is not sufficient to characterize the full frequency response of a system; you have to know (or assume) what its phase response is also.


What is differences between structured image and textured image

Imho there is no "mathematical" description of either terms.

Both (texture and stucture) can be viewed as a repeating pattern. In the case of texture, the pattern might get a more random perturbation (think at an image of sand or grass).

In the case of a structured image you might get less random and more systematical perturbation (Think of a roof picture with tiles) Sometimes more info can be extracted from structured images. Imagine that you can actually get some 3D info by examining the way the pattern varies along one direction or another.

Of course I could be completely wrong and the "structure image" has a concrete definition for the author that you quote.

Edit:

I found two examples where 3D info is extracted from textural information. Chapter 9.5 "Shape from Texture" from "Introductory Techniques for 3-D Computer Vision" by Emanuele Trucco, Alessandro Verri (1998 ed) and Chapter 12.4 3D Cues Available in 2D Images from "Computer Vision" by L. Shapiro and G. Stockman (2000 ed)

There, you have the idea of deterministic texture vs stochastic texture. Also read Jedi Commymullah answer for some other insights.

Images are free:

http://pixabay.com/en/roof-roofing-house-roof-186391/
http://pixabay.com/en/grass-grassy-stalks-green-275986/


How have common PRBS polynomials been chosen, and what are they?

The number of maximal-length sequence generators increases rapidly with the degree and so tabulating all of them is impracticable for large degree. A standard table of all irreducible binary polynomials of degree 16 or less (and a few selected polynomials of degrees up to 34) takes up 17 pages in in an Appendix of Peterson and Weldon's Error-Correcting Codes, (2nd ed.), MIT Press 1972 and they are much more densely packed than the prolix notation used in the OP's reference.  As to which ones are in common use or which have desirable properties, it depends on what is meant by desirable; mathematically, they are all equivalent, practically, there might be good reason to choose one over the other. One preferable feature might be small number of nonzero coefficients (which reduces the number of XOR gates (or instructions in software implementations).  Another might be smaller overall complexity of implementation: the primitive polynomial of degree 8 used in the Reed-Solomon code implementation in the NASA system was carefully chosen to minimize the overall complexity of the decoder (and no, it is not the first one in the Peterson&Weldon table).


Where 2-D DFT is necessary in image processing applications instead of 2-D DCT?

Any place you would like to filter an image with anti symmetric filter.


Choosing minimum sampling frequency of a filtered signal

As you say, you should select a sampling frequency at least at twice the largest frequency you wish to keep.

What you'll do is called downsampling: going from a higher to a lower sampling frequency. This process usually starts by filtering, so there is no energy above the new Nyquist frequency, and any aliasing is negligible.

Many numerical programs, such as Matlab, include the filtering step in commands such as downsample and resample. In any case, you should always look at the spectrum before and after downsampling, and make sure the process meets your requirements.


Is there a relation between abs(fft(u.^2)) and abs((fft(u)).^2

Is there a relation?

Generally no. Squaring is a highly non-linear operation and completely changes the amplitude statistics and spectrum of the original signal. A simple example: if your input consists of a few sine waves of different frequency, the output will be all positive and the spectrum will contain lines are all possible combinations of sum and difference of the original frequencies. The bandwidth will also double which in the discrete domain will often result in aliasing.

has a power law +1

I have no idea what that means.


How to prevent Octave Jumping in ACF of Speech Signal?

Octave jumps down can occur if there is more noise (or other modulation) around a sub-harmonic of the pitch frequency, than at pitch.

Octave jumps upwards can occur if the timbre contains strong overtone or harmonic sequences and weak or nearly missing spectrum at or around the fundamental pitch frequency.  Noise, or a non-flat frequency response in the acoustic environment or channel, can also add power to an overtone peak.

Windowing effects in the ACF can add noise which can increase the probability of the above estimation errors.

Humans might be infering octave by using using preceding or following pitch, timbres, and transients to error correct spurious jumps in the time-local spectrum.


In my opinion, you could use a couple of easiest and more optimized solutions, unless you need only that filter and you want to avoid third party or extra library/files:

SigPack: http://sigpack.sourceforge.net/index.html
Kfr: https://www.kfr.dev/

Cheers.


In QPSK, the bit information is embedded in the phase of the carrier signal. For example, for bits 00, the phase shift is 45degree. This means the carrier signal $x(t) = cos(2\pi f_c t + \phi)  0 \lt t \lt T$, the phase shift $\phi = \pi / 4$. So $x(t) = cos(2\pi f_c t + \pi / 4) = cos(2\pi f_c t) cos(pi/4) - sin(2\pi f_ct)sin(pi/4)$. The Inphase component corresponding to $cos(2\pi f_c t)$ has factor of $cos(pi/4)$ which is what you see in the I component. Similar logic for Q component corresponding to $sin(2\pi f_c t)$. If the bits are 10, then I component will have factor of $cos(3 \pi/ 4)$ which will be $-\frac{1}{\sqrt{2}}$ so you will see phase shift of 180 degree due to the -1 factor.


Zoom an image with respect to a given image

Since you want the image size and resolution to remain constant, the term optical magnification seems apt there.

You could then say

  Please use 20x optical magnification to identify said feature.

Please be careful that the term magnification is different from zoom.

Two different cameras with the same 5x optical zoom might produce images with different optical magnification depending on their optics.

Zoom refers to the ratio of the max to min focal length of your camera, whereas,
Magnification refers to the ratio of projection size to the real subject size.

Given a particular magnification, you could use your camera focal lengths and apertures to come up with at zoom metric for that particular camera to produce the desired results.

See this for a better explanation.


A few overview questions about FFT. Relationships between input, buffer size, and expected output?

So, I'm not an expert at all, but I feel that I've perhaps learned enough to be somewhat helpful to anyone who ends up here with little introduction to FFT.  If any of this is wrong, see the "P.S" at the bottom.  tl;dr I'd rather share than leave this question unanswered entirely.

A few things are important to note: Parseval's theorem, bin frequencies, and normalization:

Parseval's theorem is something to the ends that the energy of the signal will be the energy of FFT bins (bins being each frequency/magnitude pair).  So, what this means in practice (I think) is that a small buffer size will make the bin magnitude of a sine wave frequency vary based off the amount of high amplitude wave is in the buffer.  It won't be consistent.

Bin Frequencies are the frequency of each bin from your FFT.  They aren't going to be the exact numbers you are looking for without a bit of work, but they should be close enough to give you some info.

Normalization is something that helps with all this.  For my project, I found the max bin magnitude, averaged it with the max from the last cycle, and then divided everything by it.  The result is that my visualization is showing something much more like what my ears feel they would see.

Now for the Q's (assuming all signals are normalized to -1 < x < 1):

Magnitudes of ~200 +-100, in the bins near the signal wave freq
Same, but magnitudes of ~300 +-50.  The larger buffer limits cases of more or less high amplitude chunking (sin(x) at [0 to Pi/4] vs [Pi/4 to Pi/2])
Magnitudes of ~150 +-100, in the bins near the signal wave freq.  The sweeping spreads out the bins a little, but the small buffer size (relative to sweep speed) helps.
Magnitudes of ~50 +-100, in the bins touched by the signal wave freq.  If the sweep is larger or faster expect the magnitudes to be even more spread out.  Note: the more bins touched, the lower the magnitudes overall.
Magnitudes of ~5 +-5, in all bins which happen to be strongly represented in this small buffer chunk
Magnitudes of ~1 +-1 in all bins.

So, it seems there is probably some form of log function which would help make the perceived levels of all these sounds not so distant when graphed by magnitude, but, overall, using normalization has worked in my project to pick out the most noteworthy note tones.  (btw notes and tones are anagrams)

P.S.  If any of this is wrong, please do make an answer covering and correcting these topics.  I will not be marking this as the accepted answer for at least a month because I'm not experienced enough to be an authority on this subject.  That being said, this is what I think may be true based off the results I've been getting, and it's better than nothing for a newcommer.


find cross correlation lag by optimization

This paper has an algorithm that seems to try to do what you're reqeuesting.  I don't have time right now to look at implementing it, but may do in a few days.

The paper uses the signal model:

and the RLS algorithm to modify a filter:

I'm not sure whether this qualifies as any different, as it seems to use the peak of the autocorrelation as its delay estimate.


Role of Microphone sensitivity in building an SPL meter or dB meter

This is all about dynamic range since you can always compensate for sensitivity with a proper pre-amp. Dynamic range is limited at the low end by the self noise of the microphone (or the quality of the electronics) and the high end is limited by the distortion of the microphone or the clip points of the electronics.


What are some free alternatives to SIFT/ SURF that can be used in commercial applications?

Both SIFT and SURF authors require license fees for usage of their original algorithms.

I have done some research about the situation and here are the possible alternatives:

Keypoint detector:

Harris corner detector
Harris-Laplace - scale-invariant version of Harris detector (an affine invariant version also exists, presented by Mikolajczyk and Schmidt, and I believe is also patent free).
Multi-Scale Oriented Patches (MOPs) - athough it is patented, the detector is basically the multi-scale Harris, so there would be no problems with that (the descriptor is 2D wavelet-transformed image patch)
LoG filter - since the patented SIFT uses DoG (Difference of Gaussian) approximation of LoG (Laplacian of Gaussian) to localize interest points in scale, LoG alone can be used in modified, patent-free algorithm, tough the implementation could run a little slower
FAST
BRISK (includes a descriptor)
ORB (includes a descriptor)
KAZE - free to use, M-SURF descriptor (modified for KAZE's nonlinear scale space), outperforms both SIFT and SURF
A-KAZE - accelerated version of KAZE, free to use, M-LDB descriptor (modified fast binary descriptor)

Keypoint descriptor:

Normalized gradient - simple, working solution
PCA transformed image patch
Wavelet transformed image patch - details are given in MOPs paper, but can be implemented differently to avoid the patent issue (e.g. using different wavelet basis or different indexing scheme)
Histogram of oriented gradients
GLOH
LESH
BRISK
ORB
FREAK
LDB

Note that if you assign orientation to the interest point and rotate the image patch accordingly, you get rotational invariance for free. Even Harris corners are rotationally invariant and the descriptor may be made so as well.

Some more complete solution is done in Hugin, because they also struggled to have a patent-free interest point detector.


MATLAB and Radon Transform-cannot obtain a sinogram

There are two points that stand out from the provided code:

The determination of detector_size

You seem to be trying to track the "widest line of sums" (?). But we know that this is the diagonal and it occurs at 45 degrees in square images. Therefore, we can "pre-allocate" our detector to have that maximum size anyway.

The rotation of the image

The code seems to be along the right track, taking the "image grid coordinates" (i,iy) and rotating them to (x_img,y_img) but:

There is no centre of rotation in the calculation (?). This usually is taken to be at the centre of the image.
While some rotation seems to be taking place, the data is still read from the original locations of A (involving some scaling too) (?).

There are two approaches that you can take here:

Rotate the image
Rotate the "detector" (but do it properly)

The first approach is dead simple but a bit wasteful in terms of memory and computation. It goes a little something like this:

Given a $M\times N$ image ($I$), construct a $D \times D$ image ($A$) where $D = \sqrt{M^2 + N^2}$
Now, centre $I$ into $A$. This basically translates $\left(0,0\right)$ of $I$ to $\left(\frac{D-M}{2},\frac{D-N}{2}\right)$ in $A$.
For each  angle $\theta \in \Theta$, rotate the image (using imrotate) to produce $A_r$. Now, note here that depending on your frame of reference, you might have to transform $\theta$. That is because, if you sum the columns of the rotated image, then your detector is positioned at $+ \frac{\pi}{2}$. But that is just a simple linear transformation to align your results with convention. Also important to note here that it is worth rotating $A$ afresh each time rather than keep rotating $A_r$ by increments, which would progressively distort the image due to successive applications of interpolation.
Once you have $A_r$, produce the sum of the "detector" by adding up its columns.

And that is it.

It is dead simple because you don't have to handle the rotation of the image which, by itself, is not difficult but it involves interpolation which can get messy. The rotated image grid will never land on integer pixel locations and so you are going to have to apply some form of interpolation that goes beyond nearest neighbour. For linear interpolation (for example), this would involve determining the patch of the 4 nearest pixels to the target location, then obtaining the distances to those pixels and finally creating the weighted sum.

It is "wasteful" in CPU because $A$ is bigger than $I$ and all pixels are rotated even if their contribution to a sum is in fact zero.

It is "wasteful" in memory because $A$ is bigger than $I$ and for large differences of $M,N$, large square images are produced.

The immediate improvement you can apply to this is to rotate the "scanning grid" (as you are trying to do now) but when you obtain a slice off of the image to centre it in the "detector slice" which would have to be $D$ pixels wide anyway. In this case of course, you still have to handle interpolation. But in this approach, any coordinates landing outside of the original image are simply skipped.

Hope this helps

EDIT:

The rotation is still not performed correctly:

% Load the image
Q = imread("testImage.jpg");
% Resize the image, in case it is too big
Q = imresize(Q,0.05);
Q_size = size(Q);
% Find its diagonal
D = ceil(sqrt(sum(Q_size.^2)));
% Create a square image of size equal to the diagonal
I = zeros(D,D);
I_size = size(I);
% Center the image in the square
offset = ceil((D-Q_size)./2);
% Copy the original image across
I(offset(1):offset(1)+Q_size(1)-1,offset(2):offset(2)+Q_size(2)-1) = Q;

I contains the padded image.

Using imrotate:

% Theta angles to obtain sums at
theta = 0:11.25:180;
% Initialise the output. One row for each angle, one column for every summing "sensor".
R1 = zeros(length(theta), D);
for k=1:length(theta)
    R1(k,:)=sum(imrotate(I,theta(k),'nearest','crop'));
end;

figure;imshow(R1./max(max(R1)))

The above finishes very quickly. Notice nearest as interpolation in imrotate and crop so that imrotate does not add yet more padding.

Not using imrotate:

% Rotating the detector
theta_rad = (theta./180).*pi; %This simply converts theta (from above) in radians
% Initialise the output. One row for each angle, one column for every summing "sensor".
R2 = zeros(length(theta_rad), D);
% THE CENTRE OF THE IMAGE
I_centre = round(I_size./2);
% FOR EVERY ANGLE
for k=1:length(theta_rad)
    % m,n scan the image left to right top to bottom...BUT...
    for m=0:(I_size(2)-1)
        for n=0:(I_size(1)-1)
            % THEY GET REMAPPED TO U,V
            % Remap m,n to their rotated versions.
            % I do the rotation in polar coordinates here
            rho = sqrt(sum((I_centre-[n,m]).^2));
            phi = atan2(n - I_centre(1), m - I_centre(2));
            % This point now maps to
            u = round(I_centre(2) + rho .* cos(phi-theta_rad(k)));
            v = round(I_centre(1) - rho .* sin(phi-theta_rad(k)));
            % Check that the point is inside the image
            if (u>0 && u<I_size(1) && v>0 && v<I_size(2))
                % NOTICE THAT I IS READ FROM THE REMAPPED LOCATIONS
                R2(k,m+1) = R2(k,m+1) + I(v,u);
            end;
        end;
    end;
end;

figure;imshow(R2./max(max(R2)))


Interpolation of estimated frequency between fft bins: what about phase?

If you do an fftshift before the FFT (to center the phase reference in the data frame), you can use the same interpolation methods on the real (even) and imaginary (positive odd) components in the complex FFT result, and thus get interpolated phase estimates (referenced to the original data frame’s center) between FFT result bins.


OFDM demodulation, how to set the sampling rate?

Your receiver will consist of, apart from synchronization blocks, an FFT. You want that FFT to generate one output bin for each of the subcarriers. Therefore, you should make use of the relationship between the input sample rate of a DFT to its output bin spacing:

$$
\Delta f = \frac{f_s}{N_{FFT}}
$$

where $\Delta f$ is the subcarrier spacing and $f_s$ is the sample rate (both in Hz). Assuming you know the subcarrier spacing and the FFT size, then it's straightforward to choose the sample rate:

$$
f_s = N_{FFT} \Delta f
$$


Mathematically, what you probably want is the autocorrelation function. That function is a function on a function over time or a sequence, which gives you info about how similar your signal is to itself when shifted by shift $\tau$.

If your signal is periodic with period $\tau$, then it will be very similar to itself when shifted by $\tau$, $2\tau$, $3\tau$, and so on. So, you'll get $\tau$-periodic peaks in your autocorrelation function

A very common thing to do is to look at the Power Spectral Density (PSD), which is, in fact, the Fourier transform of the autocorrelation function. You can, however, in your sampled case, calculate the PSD directly by Fourier-transforming (Matlab function fft) your sequence and finding the abs² of it. The peaks correspond to periodicities.


Confusion on using FFTW 3.3.4 (C version) for real and even data sets

If the imaginary part is very small you are probably not doing anything wrong.  This is a common occurrence due to floating point quantization noise.  As long as you know the result should be real and the imaginary part is small, just drop the imaginary part.


Comparing Welch and FFT power spectrum density

scipy.signal.welch will, by default, detrend the data by subtracting the mean of each segment. This DC level suppression will only affect bin zero of the resulting spectrum. This is very useful since the DC offset can be so much larger than the noise that it wouldn't fit on a reasonably scaled log plot.


Identifying steps in walking accelerometer

It would be simple to make an acquisition of a single step and identify the phases of the step. For that, you shall consider the geometry of your setup, and therefore answer these:

The accelerometer has (or not) a relative speed to the ground equal to ZERO?
The accelerometer measurement axis is aligned to gravity?

My guess is that, when one foot is on the ground, the accelerometer shall measure only gravity for a short period. Therefore, I would look for segments where the jerk is close to ZERO in the vertical axis, AND the magnitude of it is close to -1g.

EDIT:
All this new information you have given us, like sensor close to chest, gravity, alignment, etc, demands new dynamic models to represent the measurements you have there. But any model, with respective solution, can have counter-examples to demonstrate that the model is inadequate at all. For example, if the user carries a heavy rucksack, there is a new model to measure the acceleration close to the neck. This is different from the acceleration tied to the feet. That said, the body shall damp the actual acceleration at the feet to a filtered and delayed measurement at the neck. Nevertheless, the readings shall have the same low frequency components as at the feet. This is more less the same I see at your graph, the one with blue lines. The blue lines are as periodic as the low freq. components of the red line, but for a phase. You should consider to review your necessity of precise instant of touch on the floor in accordance with your application.


What is going wrong with my IDFT program?

I suspect reason is that my IDFT program takes into account every leakage frequency.

Also in complex numbers multiplication I do not take into account imaginary unit.

I apologize for the errors.

My answer to my question is the following C++ program. It performs DFT and IDFT on 20 input values in array, and outputs the same values. Sample rate for calculate input signal from existing sinusoids is 71114 Hz, and N = 20. Existing sinusoids are two:

1. Amplitude = 995,574, frequency = 20247 Hz, phase = 0.34 π.
2.    "      = 492,356      "     = 3818 Hz,    "   = 0.15 π.

Sample rate for DFT and IDFT is 45940 Hz, and N = 20. I point out that existing frequencies are not multiples of sample rate by N. So, there is leakage but is not problem. Any one can change values. These in program are random in limits. Any one can use first program in my question, for calculate input x(n) from any sinusoids.

     #include <iostream>
     #include "math.h"

     int main()
     {
     //Unchanged variables.
     float pi = 3.141592653589793; //pure number.
     int n = 0; // sample.
     float Re = 0; //pure number.
     float Im = 0; //pure number.
     float xn_cos = 0; //pure number.
     float xn_sin = 0; //pure number.
     float cosine = 0;
     float sine = 0;
     float sumReal = 0;
     float real = 0;
    int m = 0;
    int no = 0;

    //Sampling frequency and samples of segment.
    int samplFreq = 45940; // sample/second.
    int N = 20; //sample. Number of input samples.

   // Testing grid.
    float testFreq = 0; // Hz.
    float testUnitAngle = 0; // rad/sample.
    float testAngle = 0; // rad.
    float lowTestFreq = 0; //Hz.
    float highTestFreq = samplFreq; //Hz.
    int M = N; // pure number. Number of DFT output samples.
    float testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.
    float xn[] =
{ 1.09595e+06,
 635609,
-544799,
 640870,
 1.40537e+06,
-142213,
-392826,
 1.00018e+06,
 307386,
-1.16555e+06,
-211430,
 518555,
-998272,
-1.20759e+06,
 390434,
 26244.3,
-1.18203e+06,
 25149,
 1.08298e+06,
-200068,};
   int No = N; // pure number. Number of IDFT output samples.

  // IDFT including DFT.
  for (no = 0; no <= No - 1; no++)
    {
            sumReal = 0;
            for (m = 0; m <= M - 1; ++m)
            {
                    testFreq = lowTestFreq + m * testFreqStep; // Hz.
                    testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

                    Re = 0; Im = 0;

                    // Start DFT.
                   for (n = 0; n <= N - 1; ++n)
                    {
                            testAngle = n * testUnitAngle;

                            xn_cos = xn[n] * cos(testAngle);
                            xn_sin = xn[n] * -sin(testAngle);

                            Re += xn_cos;
                            Im += xn_sin;
                    }
                   // End DFT.

                   cosine = cos(((no * 2 * pi) / N) * m);
                    sine = sin(((no * 2 * pi) / N) * m);
                    real = Re * cosine - Im * sine; // Product of complex
//numbers (Re, Im i) * (cosine, sine i) = (Re * cosine + Re * sine i + Im i *
//cosine + Im i * sine i) = (Re * cosine + Re * sine i + Im i * cosine - Im *
//sine), for, by definition, i * i = -1. Then, first and last product in
//parenthesis are real numbers. Other two are imaginary and their sum is
//practically zero. For the values in array are between -26 and 25.
//More over because  input values are real, output should be real.
                    sumReal += real;
            }
            std::cout << sumReal / N << "\n";
    }
    return 0;

}


Measuring signal phase difference between receivers

Preface

First of all: Arduino might really not be the tool of choice here; with a 1MS/s max, and some something-in-between-8-and-16-effective-bits ADC, there's hard limits on what you can detect. I will explain the theory of how to do better estimation than just sample-wise, but I won't address these limits here.

Furtherore, I assume you

do real-valued sampling with the built-in ADC (because that is not a simultaneous two-channel/IQ ADC), meaning that you probably have some IF receivers that convert a RF bandpass signal to a directly sample-able IF signal (or one that you can undersample)
keep true to the filtering requirement of Nyquist's sampling theorem,
use the actual received signals, and not just e.g. the time of their amplitude peak, and
have knowledge of DSP basics and some understanding of the DFT.

Confirmation of your Calculation

So, your 300m calculation is correct if you can only compare two sample streams only by shifting one of them integer amounts of samples. And with this, I assume you compare using a cross-correlation function (ie. you find the product of one reception and the other, shifted by all reasonable full-sample delays, and look for a maximum).

But you don't have to do that; you can either

interpolate in time domain and calculate the cross-correllation in time domain, and find the argmax of that
compare the two sequences in frequency domain and estimate the time difference directly.

Interpolation in Time Domain

Since we know that you've adhered to Nyquist sampling theorem (you did, didn't you?), we want to keep it that way. To properly interpolate by a factor of $N$, you just pad every sample you get with $N-1$ zeros, and then low-pass filter the result to the original band (ie. apply an $\frac 1N$th band-filter). That gives you a $N$ times higher resolution.

However, it also means that your correlation  just went in complexity from $L^2$ ($L$ being the length of your sample sequence) to $N^2L^2$, which is really unfortunate on an Arduino. So you'd have to do that on a PC, probably.

Frequency Domain Techniques

You're probably aware of the time-shift property of the Fourier or the Discrete Fourier Transform (DFT). For the latter it states that

$$
\mathcal F\{x[t-\tau]\}(f) = e^{-j\frac{2\pi}{L} \tau f} F\{x[t]\}(f)
$$

So let's name the signal that the "closer" receiversees $x_0$, and let's define the other receiver sees, $x_1$, to be $\tau$ "late"; it directly follows that
$$\begin{align}
\mathcal F\{x_1[t]\}(f) &= \mathcal F\{x_0[t-\tau]\}(f) \\
&= e^{-j\frac{2\pi}{L} \tau f} F\{x_0\}(f) \\
&\implies\\
\frac{\mathcal F\{x_1[t]\}(f)}{\mathcal F\{x_0[t]\}(f)} &=e^{-j\frac{2\pi}{L} \tau f}\\
&=: \varphi(f)
\end{align}$$

I introduced the "auxillary" function $\varphi(x)$ here – so if you look at that,

$$\varphi(x)=e^{-j\frac{2\pi}{L} \tau x}\text,$$

you'll notice that it's a complex sinusoid with frequency $\tau$!

So, by estimating the frequency of that, you get $\tau$,  your time delay. Notice that the variance of a frequency estimate reduces with the number of samples you observe – in other words, the longer your DFT (you'd probably use an FFT), the higher your frequency resolution, and the smaller your delay granularity, and the better your Time-Difference of Arrival estimate.

Remarks

Sampling Rate of $\frac1{1\text{ µs}}= 1\text{ MHz}$: That sounds too high, actually. Where did you get that number? I ask because I have my doubts: what are you going to do with that amount of data inside an Arduino? Spit it out via the serial port? Certainly not, that's much slower. Process it on-Arduino? Don't think so, unless I underestimate the processing power of an ATMega seriously. Store it in RAM and use it later? Well, at 1MS/s, how long will your RAM last?

Correlation in time Domain is seldom done directly, but often by (conjugate) multiplying the DFTs of the sequences, and then reverting back to time domaind – an FFT has a complexity of $\mathcal O(N\log N)$, and thus the whole "transform both sequences – multiply – inverse transform" has complexity of about $\mathcal O(3N\log N+ N) "<" \mathcal O(N^2)$. When you write the fast correlation down, you'll notice similarities to the frequency-domain approach above, if you apply the shift theorem.

Variance of the $\tau$-estimator can't be pushed down indefinitely with real-world signals, since well... real-world signals do have finite length.


How to find analytic description of filtered signal

Assuming your algebra is correct, there is nothing wrong with complex numbers in the Laplace domain. Complex poles reflect an oscillatory signal - there is no oscillation at all if all poles are strictly real, just exponential growth or decay. The variable $s$ is itself a complex number, and the transfer function is defined over all $s$, i.e. the entire complex plane.


Feedforward disturbance signal compensation

1 - What transfer function are you trying to implement? The controller transfer function or the "plant" transfer function?
2 - You are trying to simulate some kind of inverter/active rectifier, right? Inverters are non-linear, you can simulate them with only linear blocks if you want but you will not be able to test your controller with the non-linearities that are present in real-life.
3 - What are P1 and P3 ?
4 - You kind of compensate your cross-coupling terms in the second image the $\omega Li_q$  and  $\omega Li_q$ terms. The cross-coupling terms are compensated to improve the response of your controller, but they are not technically needed.
5 - In a real-life implementation you don't have to compensate the cross-coupling terms but compensating the cross-coupling makes your control system settle faster if I recall correctly.


What is the relationship between angular frequency and normalized angular frequency

Again there is no wrong or right here. In the Alan Oppenheim's Discrete-Time Signal Processing book, the notation is as follows:

when there are only continuous-time signals we use $\omega$ for radians per second frequency.
when there are only discrete-time signals we use $\omega$ for radians per sample frequency
when both types of signals are present, (as in sampling), we use $\Omega$ for the continuous-time radians per second frequency, and $\omega$ for the discrete-time radians per sample frequency.

Furthermore the relationship between the two frequencies because of sampling normalization is:

$$ \Omega = \frac{\omega}{T_s} $$

or equivalently

$$ \omega = \Omega ~ T_s $$

However, your instructor seems to prefer the opposite notation for frequencies. That's why most literature seems to be reciprocal of your instructor's.


Why is $\delta[an]=\frac 1 a \delta[n]$ in discrete time?

It's actually good that you couldn't prove it because the claim is wrong. Note that in discrete time we have $\delta[n]=1$ for $n=0$ and $\delta[n]=0$ for $n\neq 0$. So if $a\neq 0$, you simply get $\delta[an]=\delta[n]$ (note that $a$ must be an integer for this to make sense). If $a=0$ then $\delta[an]=1$ (for all $n$).


How to sketch/plot this discrete time signal?

Begin with defining the base signal as $x[n] = n ~u[n]$ and recognize the output as a linear compositions (superpositions) of the base signal as:
$$y[n] = x[n+5] + 2 x[n] + x[n-10]$$

Now what you have to answer is the result of translations and scalings applied on the basic signal $x[n]$. This ability is a simple yet an indispensibly necessary skill for a clear understanding and successful applications of all subsequent signal processing theory. Therefore you must train yourself enough for this task.

Incidentally, the signal $x[n+5]$ is an advanced (shift left) version of $x[n]$ whereas $x[n-10]$ is a delayed (shift right) version.

Graphically the result is obtained by summing the three components drawn on the on the same plot. A more detailed approach would use the piecewise definitions of the components to define the sum. It would proceed as follows:

Let $x[n]$be defined as:
$$ x[n] = n ~\text{u}[n] =
\begin{cases}
n , \text{ for } n \geq 0 \\
0 , \text{ for } n < 0 \\
\end{cases}
$$

Then you can define the three superposed signals as:

$$
\begin{align}
x_1[n] = x[n+5] = (n+5) ~u[n+5] &= \begin{cases} n+5, &\text{ for } &n+5\geq 0 \\ 0, &\text{ for } &n+5< 0 \\ \end{cases} \\
x_2[n] = 2x[n] = 2n ~u[n] &= \begin{cases} 2n,~~~~ &\text{ for } &2n \geq 0 \\ 0,~~~ &\text{ for } &2n < 0 \\ \end{cases} \\
x_3[n] = x[n-10] = (n-10) ~u[n-10] &= \begin{cases} n-10, &\text{for} &n-10 \geq 0 \\0, &\text{for} &n-10<0 \\ \end{cases} \\
\end{align}
$$

Modifiy the intervals:
$$
\begin{align}
x_1[n] = x[n+5] = (n+5) ~u[n+5] &= \begin{cases} n+5, &\text{ for } &n \geq -5 \\ 0, &\text{ for } &n<-5 \\ \end{cases} \\
x_2[n] = 2x[n] = 2n ~u[n] &= \begin{cases} 2n,~~~~ &\text{ for } &n \geq 0 \\ 0,~~~ &\text{ for } &n<0 \\ \end{cases} \\
x_3[n] = x[n-10] = (n-10) ~u[n-10] &= \begin{cases} n-10, &\text{for} &n \geq 10 \\0, &\text{for} &n<10 \\ \end{cases} \\
\end{align}
$$

Therefore $y[n] = x_1[n] + x_2[n] + x_3[n]$ is:
$$
y[n] =
\begin{cases}
0 , &\text{ for }     &-\infty<n<-5 \\
n+5 , &\text{ for } &-5 \leq n <0 \\
3n+5 , &\text{ for } &0 \leq n <10 \\
4n-5 , &\text{ for } &10 \leq n <\infty \\
\end{cases}
$$


For any convolution matrix (even truncated), your diagonal entries would be $h_0$; hence, the characteristic polynomial, setting $b_0 = h_0-\lambda$,

\begin{align}|H-\lambda I| = \begin{vmatrix} b_0 & 0  & 0\\h_1 & b_0 & 0\\h_2 & h_1 & b_0 \end{vmatrix} &= b_0\,\begin{vmatrix} b_0 & 0\\h_1 & b_0 \end{vmatrix} - 0\,\begin{vmatrix} h_1&0\\h_2&b_0\end{vmatrix} + 0\,\begin{vmatrix} h_1 & b_0\\h_2&h_1 \end{vmatrix} \\ \\
&=b_0\,\begin{vmatrix} b_0 & 0\\h_1 & b_0 \end{vmatrix}\\
&=b_0\left(b_0\,b_0 - 0\,h_1 \right)\\
&=b_0^3
\end{align}

Little induction is necessary to see that for any square matrix of dimension $N\times N$ with an lower left structure and a constant $d$ in each element of the diagonal, the characteristic polynomial is $(d-\lambda)^N$; the Eigenvalues are roots of that polynomial, and hence there's but one Eigenvalue (with $N$ multiplicity): $d$, which in your case is $h_0$.

So, since not all zeros of $H(z)$ are defined by your first filter tap, no, you can't find all zeros of $H(z)$ through Eigenvalues.


What are the minimum knowledge required to begin the study of Digital Signal Processing?

I personally believe that the approach of studying something just in order to "feel ready" to study something else later on is not efficient and tends to overwhelm the student. (This of course only applies because I assume that you know what a computer scientist should know anyway: basic math, and - most importantly - how to think clearly and logically.) I could give you a whole list of mathematical topics which are important for doing DSP, but this would just paralyze you and would make you feel that there is so much you need to know that you won't even know where to start. Note that DSP has become such a huge topic in itself that there is hardly anyone who can claim that he knows even the basics of all its different sub-fields.

So what do I recommend to you after all? Get a good book on DSP and simply start studying it. If you come across problems that you cannot solve yourself, dig deeper by searching the internet, asking questions on this site, etc. In this way you will learn all necessary basics that you might be lacking. But the point is that only in this way will you find out what it actually is that you are lacking, so filling the gaps in your knowledge will be most efficient.

One good (and free!) book I would like to recommend to you is Introduction to Signal Processing by Orfanidis. It contains a lot of necessary basics and it is very practical. I think you will learn a lot in quite a short time. Also take a look at this online-book.

The answers to this question might also be relevant for you.


Why the curvature of straight line in curvature scale space is not zero?

My guess is you're calculating the derivatives wrong, or you're applying the whole thing to the wrong input data. I can show you the expected output for each step graphically, so you can compare it to what you get. I've implemented this in Mathematica; I'll post the code at the end, but I don't know if it'll help you much - Mathematica has built in functions for e.g. convolutions, that you would have to write yourself in C++.

First the input data for a rectangle would look like this:

That's 150 points along a rectangular path, using linear interpolation between the corner points.

Next, the Gaussian and it's derivatives (using $\sigma=3$ for illustration):

If you convolve these with the X and Y coordinates separately, you get this:

The first row is just for illustration - you don't actually need the convolution with the Gaussian itself, only the two derivatives.

At the sides of the rectangle, the first derivatives (second row) are constant and the second derivatives (third row) are 0. (You might want to write the results you get to a text file and look at them e.g. using Excel to check this.)

Then, using the curvature formula you quoted, you get this:

Along the straight sides of the rectangle, the curvature is 0 as expected. At the corners, the curvature of a rectangle is infinite; however, what we're calculating here isn't the real curvature, but the curvature of a Gaussian-smoothed rectangle, with filter size $\sigma=3$:

To check this result, I'll draw a circle with radius $r=\frac{1}{\kappa }$ at each point:

Mathematica code:

(* interpolate between corner points *)
corners = {{-2, -1}, {2, -1}, {2, 1}, {-2, 1}, {-2, -1}};

dist = Rescale@Prepend[Accumulate[Norm /@ Differences[corners]], 0];

int = Table[
   Interpolation[Transpose[{dist, corners[[All, i]]}],
    InterpolationOrder -> 1], {i, 1, 2}];

pts = Outer[#1[#2] &, int, Range[0., 1 - 10^-9, 1/150]]\[Transpose];

cornerIdx = Most[dist*Length[pts] + 1];

Row[{
  ListLinePlot[pts, AspectRatio -> Automatic, Mesh -> All,
   ImageSize -> 400, PlotLabel -> "Points",
   Epilog -> {Red, PointSize[Large], Point[pts[[cornerIdx]]]}],
  Spacer[50],
  ListLinePlot[Transpose[pts], PlotRange -> All,
   GridLines -> {cornerIdx, {}}, PlotLabel -> "X/Y Coordinates",
   PlotLegends -> {"X", "Y"}, ImageSize -> 400]
  }]

(* Gaussian & derivatives *)
\[Sigma] = 2;

g0 = Table[
   1/(Sqrt[2 \[Pi]] \[Sigma])*Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g1 = Table[-x/(Sqrt[2 \[Pi]] \[Sigma]^3)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g2 = Table[(-\[Sigma]^2 + x^2)/(Sqrt[2 \[Pi]] \[Sigma]^5)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];

ListLinePlot[{g0, g1, g2}, PlotRange -> All,
 PlotLegends -> {"Gaussian", "1st derivative of Gaussian",
   "2st derivative of Gaussian"}]

(* Convolve x and y with gaussian&derivatives *)
derivatives =
  Outer[ListConvolve[#1, #2, Round[Length[g1]/2] + 1] &, {g0, g1, g2},
    Transpose[pts], 1];

Grid[MapThread[
  ListLinePlot[#1, PlotRange -> All, PlotLabel -> #2,
    GridLines -> {cornerIdx, {}}] & ,
     {derivatives,
    {
    {"X (smoothed)", "Y (smoothed)"},
    {Overscript[X, \[Bullet]],
     Overscript[Y, \[Bullet]]}, {Overscript[X, \[Bullet]\[Bullet]],
     Overscript[Y, \[Bullet]\[Bullet]]}}}, 2]]

{{xs, ys}, {dx1, dy1}, {dx2, dy2}} = derivatives;

(* calculate curvature *)
\[Kappa] = dx1*dy2 - dx2*dy1/((dx1^2 + dy1^2)^(3/2));
ListLinePlot[\[Kappa], PlotRange -> All,
 PlotLabel -> "Curvature \[Kappa]", GridLines -> {cornerIdx, {}}]

ListLinePlot[Transpose[{xs, ys}], Mesh -> All]

(* osculating circle animation *)
animationIdx =
  Select[Range[Length[xs]],
   Abs[\[Kappa][[#]]] > 10^-2 || (Mod[#, 5] == 0) &];

Monitor[frames = Table[Rasterize@Row[
      {
       ListLinePlot[Transpose[{xs, ys}], Mesh -> All,
        ImageSize -> 400,
        PlotLabel -> "Osculating circle r=1/\[Kappa]",
        Epilog -> Module[{r, pt, n},
          (
           pt = {xs[[i]], ys[[i]]};
           r = Clip[1/\[Kappa][[i]], {-1, 1}*10^4];
           n = Normalize[{dy1[[i]], -dx1[[i]]}];
           {Red, PointSize[Large], Point[pt], Circle[pt - n*r, r]}
           )]
        ],
       ListLinePlot[\[Kappa], PlotRange -> All,
        PlotLabel -> "Curvature \[Kappa]",
        GridLines -> {cornerIdx, {}}, ImageSize -> 400,
        Epilog -> {Red, Line[{{i, 0}, {i, 10}}]}]
       }], {i, animationIdx}], i];


Impulse response measurement : practical questions

If this had to be the approach to measure the impulse response, here are my further suggestions specific to the qeustions:

How to determine the appropriate excitation signal duration?

The equivalent noise bandwidth of each FFT bin if not windowed is $1/T$ where $T$ is the duration of the signal in time prior to taking the FFT (note zero padding will not change this, but windowing will). So this will be the equivalent to the resolution bandwidth on network and spectrum analyzers, and will effect the resolution of the resulting transfer function in frequency. Resolution is not equivalent to number of samples, since we can zero pad in time which will result in more samples in frequency (over the same frequency duration) but this will just be a smoothing or interpolation of the same frequency domain result and will not help us further distinguish a closely spaced variation in the transfer function that may actually exist. To be able to measure such a tight variation in frequency, we need a longer time sample. Ultimately this is the expected length of the impulse response.

Do you recommend excitation signal windowing?

This depends on the excitation signal and the dynamic range of the resulting frequency transfer function. It would generally be a good idea but significantly complicates the processing; for example if you did end up using a swept sine as the excitation signal, the window will be significantly attenuating all frequencies represented at the start and end of the sweep, thus reducing the SNR and therefore quality of the results in the frequency domain transfer function for  those frequencies. Windowing decreases the frequency resolution of the main lobe of each DFT bin's equivalent filter at the benefit of significantly decreasing the sidelobes and thus improving overall dynamic range. Please see this post further detailing windowing considerations as this is what would all apply here:
How to calculate resolution of DFT with Hamming/Hann window?
Why would one use a Hann or Bartlett window?

How to determine appropriate output signal measurement duration?
Before applying fft, how much zero padding do I need to add to
excitation signal and output signal?

Both these questions affect the zero padding since we need the FFT's to be of the same duration prior to computing the frequency domain transfer function. The output signal will be the duration of the input signal plus the length of the impulse response. We already established that the input signal should be the duration of the input response, so in this case it would be zero padded out an additional amount. There is no harm in having a longer duration sequence, these are the minimums so if I were to do this it would likely be an iterative approach where I would evaluate both the sampling rate and signal duration from the frequency domain transfer function and time domain impulse response respectively: In the frequency domain, if the transfer function has not sufficiently attenuated by the time we reach half the sampling rate, aliasing will occur and the sampling rate should therefore be increased. In the time domain if the impulse response has not sufficiently attenuated by the time we reach the end of the duration of the resulting answer, time domain aliasing will also occur and the overall time duration of the signal should be increased.  To summarize with consideration to signal duration and zero padding: The signal duration must be longer in time than the impulse response duration. This signal must also be zero padded out further by the duration of the impulse response duration to allow for the output to fully settle. The output duration will match the input duration and no further zero padding will be required. The sampling rate must also be confirmed to be sufficiently high enough by reviewing the resulting frequency transfer function.

How to avoid zero division when excitation signal max frequency <
f_s/2 ?

Since you are creating the excitation signal you will know precisely if any frequency is 0, so this would be trivial to avoid. What I would be tempted to do personally is decide what dynamic range my solution needs to cover for the resulting frequency domain transfer function and then replace any frequency domain zeros with a very small number that I know is much smaller than that. However having an excitation signal that had lower energy in any given bin would be a poor choice for determining the frequency transfer function, which is required to be accurate to determine the impulse response. Granted in the design of this you would want to have an excitation signal with a bandwidth that sufficiently exceeds the bandwidth of the transfer function, in which case we could have its energy taper off at the higher frequencies that we know are well beyond the transfer function's bandwidth--and then it that case we would simply ignore (zero) the output beyond that frequency where the result is meaningless (and we wouldn't want it to be applied to the IFFT).
OTHER APPROACHES
There are a lot of drawbacks and challenges to determining an impulse response using this approach similar to the pitfalls in designing filters using the frequency sampling approach detailed here:
Difference between frequency sampling and windowing method
Alternate methods that would be more straightforward would be to slowly sweep a frequency tone in the time domain with a rate less than the resolution bandwidth of the measurement desired while comparing the amplitude and phase of the output and input directly in time (as done in a network analyzer). Alternatively but similarly to step through each frequency and compare the static results for each one.    Another approach would be to use a PRN sequence as the excitation tone or any other pseudo-random white noise source (random number generator) and use the Wiener-Hopf equations to determine the frequency domain transfer function which results in the least-mean-square solution for the transfer function. This is further detailed at this post with practical MATLAB code specific to doing this:
How determine the delay in my signal practically


Notch Filter - Bandwidth vs pole / zero location vs Q factor

In this paper by Nehora, "A Minimal Parameter Adaptive Notch Filter With
Constrained Poles and Zeros" IEEE Trans ASSP, vol 33, No 4, 1985 link

The Bandwidth is given as
$$BW=\pi(1-\rho) $$
where $\rho$ is the radius of the pole. The zero is assumed to be on the unit circle.

Excerpt from the relevant section of that paper:

And the references 1 through 3:

Another paper gives the same bandwidth and it is listed as the 3dB bandwidth - see (Yong, Kechu, Xiaofan, "A gradient algorithm for complex adaptive notch filters" - Intl Conf on Signal Processing, 2004) link.

Yet another paper with bandwidth details is by Regalia,"A Complex Adaptive Notch Filter" IEEE Sig Proc Letters, Vol 17, No 11, 2010, Link.
Note that Regalia's notch filter is based on All-pass filters and so it has a different equation for bandwidth.

The relevant section from the Regalia paper:


What is a window function in DSP and why do we need it?

The FFT computes the coefficients for a periodic, ever repeating signal. In practical applications, signal segments are finite and not periodic. Just applying the FFT to a segment of a signal essentially wraps the end of the segment around to the start, generating a jump discontinuity. Such jumps result in an undesirable background artifacts in the FFT amplitudes. To make the wrap-around smooth, one can fade the signal to zero at both segment ends by multiplying it with a window function.

If the signal is essentially smooth, then even a simple linear fade to zero may not enough, since also kinks in the signal lead to some background artifacts in the amplitudes, however much smaller than those resulting from jumps. The choice of the windowing function is thus a balance between background artifact reduction and computational simplicity, and also of the length of the segment.

Additional considerations are required for overlap-add schemes, where the shifted window functions have to add to one.


If $\mathbf{w_k} $ was temporally correlated, it would have memory and would be equivalent to an augmented state.  In other words, it can be written as a state variable with its own independent noise term. You would be able to rewrite the equation to the same form you originally showed.  It is the most basic form of the model.

The symmetry of $Q_k$ is because $E\{ ab\}=E\{ba\}$ for any terms $a$ and $b$ of the vector $\mathbf{w_k}$

Edit to address comments.

Let $\mathbf{x}=\mid x_0,x_1, \dots x_{N-1} \mid^T$ be a random Gaussian distributed random vector.
$$
E\left\{ \mathbf{x} \mathbf{x}^T \right\}=\mathbf{R} \quad \text{Is typically Full rank and also typically} = \sigma^2 \mathbf{I}
$$
when the elements of $\mathbf{x}$ are independent.
There are some exceptions to full rank, like like a linear dependence of some of the elements of $\mathbf{x }$ such as $x_i=x_j$ of more generally $x_i=\sum_{k=0, k\ne i}^{N-1} \alpha_k x_k$ .

The easiest example is to consider $E\{ x_i x_i \} = E\{ x_i^2\}$  .  A Gaussian random variable can have zero mean but a variance $\sigma^2$. The matrix $\mathbf{R}$ would have nonzero diagonal terms, and is Full Rank.

The Matrix $\mathbf{R}$ isn't required to be diagonal and has nonzero off diagonal terms when $E\{x_i x_j\} \ne 0$ for some $i$ and $j$.  For your problem, your book will likely show the covariance of the evolution of the the state vector and that will most likely have off diagonal nonzero elements.
$$
E\left\{ \mathbf{x}\right\} E\left\{ \mathbf{x}^T \right\} \quad \text{Is rank one}
$$


FM modulation audio 50 Hz – 16 kHz MATLAB

From a practical point of view, the maximum value that your microphone will send you is limited by the input ADC clipping level. The clipping level is probably stated in terms of the analog reference voltages that the audio codec uses. Moreover, the analog value that's converted into digital may be further normalized or scaled by the audio codec firmware before sent to the calling application.

In the old PC WAV audio standard, 8-bit channels sent a maximum value of integer 255 with a bias of 128. The 16-bit channels used a signed format with maximum corresponding to $2^{15}$. Note that this value is typically mapped to some range and converted to floating point values for further processing such as done in MATLAB. So be aware of that.

The standard DSP diagnosis procedure recommends a test with a known or predictable input signal to observe the mathematical range of the resulting input samples.


LowPass Filter & Window *Img & Code* - Frequency Response has no lobe roll off

You need to zero pad the FFT so that it reveals more of the frequency response (increasing the resolution of the plot). I corrected your code (while taking fft, I took 1024 point FFT for windows as well as h_lpf)

%% Low Pass Filter
Fs=10*10^3; cutoff = 1000; %100 Hz
fc = cutoff / Fs ;
b = 800/Fs ;%Transition bandwidth of 800 Hz
M = 51 ; %Length of Filter
n = 0:1:M-1;
blackman_window = 0.42 - 0.50 * cos((2*pi*n)/(M-1)) + 0.08*cos( (4*pi*n) / (M-1));
hamming_window  = 0.54 - 0.46 * cos((2*pi.*n)/(M-1));
hanning_window  = 0.50 - 0.50 * cos((2*pi.*n)/(M-1));
blackman_window_ft = fftshift(fft(blackman_window,1024)/M);
hamming_window_ft  = fftshift(fft(hamming_window,1024)/M);
hanning_window_ft  = fftshift(fft(hanning_window,1024)/M);

sinc_func = sin(2*pi*fc.*(n-((M-1)/2))) ./ (2*pi*fc.* (n-((M-1)/2)));
sinc_func((M+1)/2) = 1 ;
sinc_func_ft = fftshift(fft(sinc_func,1024)/M);

h_lpf = sinc_func .* blackman_window;
h_lpf = h_lpf ./ sum(h_lpf); %Normalise
H_lpf = fft(h_lpf,1024);
H_lpf = fftshift(H_lpf);

figure(1);cla;clf
subplot(3,2,1);hold on;
plot(blackman_window,'DisplayName','Blackman')
plot(hamming_window,'DisplayName','Hamming')
plot(hanning_window,'DisplayName','Hanning')
plot(sinc_func,'DisplayName','Sinc')
legend
title('Windows and Sinc Impulse Response');hold off

subplot(3,2,2);hold on;
plot(20*log10(abs(blackman_window_ft)),'DisplayName','Blackman')
plot(20*log10(abs(hamming_window_ft)),'DisplayName','Hamming')
plot(20*log10(abs(hanning_window_ft)),'DisplayName','Hanning')
plot(20*log10(abs(sinc_func_ft )),'DisplayName','Sinc')
legend
title('Window and Sinc Frequency Response (dB)');hold off;

subplot(3,2,3);hold on;
plot(h_lpf);
title('Low Pass Filter Impulse Response');hold off;

subplot(3,2,5);hold on;
plot(abs(H_lpf),'DisplayName','Magnitude');
plot(real(H_lpf),'DisplayName','Real');
plot(imag(H_lpf),'DisplayName','Imag');
legend
title('Low Pass Filter Frequency Response');hold off;

subplot(3,2,6);
plot(20*log10(abs(H_lpf)));
title('Low Pass Filter Frequency Response (dB)');


Spectral Property of Bandlimited Quantized Signal

I guess this probably just a mistake in your analysis code.
Quantization noise is white and for a noise signal it's uncorrelated to the original signal, so spectrum of the original noise doesn't matter.
I did repeat your steps and saw exactly what I expected: The quantization noise is white and the spectrum of the quantized signal follows the original signal at low frequencies and is white at high frequencies (since it's dominate by the quantization noise).
Potential errors in your code

Clipping the signal
Improper spectral estimation, i.e. wrong  choice of window and/or overlap. Not windowing or properly during spectral analysis. Either window manually or just use PWELCH instead.

Here is a MATLAB code snippet that shows this
%% script to generate quantization noise
nFFT = 16384;
fs = 48000; % sample rate
fc = 4000;  % cut off for low pass filter

% create noise, filter, and normalize
x0 = randn(10*nFFT, 1); % 10 frames or so
[b,a] = butter(4,fc*2/fs);
x0 = filter(b,a,x0);
x0 = x0./max(abs(x0));

% quantize to 8 bit
xq = round(2.^7*x0)*2^(-7);

% calculate quantization noise
noiseSignal = xq-x0;
% power spectrum
yAll = [x0 xq noiseSignal];
psd = pwelch(yAll,hamming(nx));
freqAxis = fs*(0:nFFT/2)'/nx;
% graph
clf; figure;
plot(freqAxis,10*log10(psd));
title('Power Spectra');
ylabel('Level in dB');
xlabel('Frequency in Hz');
legend('Original','Quantized','Noise');
grid ('on');


Should I (and if yes, how can I) assure channel-dependant temporal accuracy in a filterbank?

Using different filter lengths for different channels is basically what the Wavelet Transform does. I think it basically addresses this issue. Wavelet transforms divide the spectrum up into spectral bins of unequal width. There are efficient options for computing the Wavelet Transform (a tree structured filter bank approach usually).

In the linked Wikipedia Article, you can see a picture comparing the time-frequency tiling of the STFT to that of the Wavelet Transform. It is directly applicable to your question.

I realize that I addressed the "How can I" part of your question, but not the "Should I" part... I'm not sure I can answer the "Should I" part unless I am given more information regarding the target application. So, I will briefly say that Wavelets have found widespread use in compression and image recognition. They are often helpful for multiresolution analysis which (in rough terms) is when one wishes to find something of a certain shape but has no idea how large it will be (for example a pulse of a certain shape that may have been compressed or expanded significantly). If you don't require these types of properties, then a uniform filter bank is often the way to go.


How to find lattice coefficients of all pole IIR filter with filter coefficient $k_m$ is 1?


Showing frequency and amplitude after doing an FFT on a signal

A strong amplitude response at 0 Hz simply means that you have a very strong DC offset. In other words, it just means that the mean of your signal is not 0.

If this is the only problem you have, then all you really need to do, is remove the mean of your signal. In other words:

vp_sig_orig = vp_sig_orig - mean(vp_sig_orig);


Spectral density of signal written as quadratures

This answer is related to my answer to this question of yours. The short answer is that in general the random process $f(t)$ has no power spectral density because it is not wide sense stationary (WSS), even if $I(t)$ and $Q(t)$ are.

I will first introduce some notation. I use the following definition of the autocorrelation function of a (possibly complex) random process $X(t)$:

$$R_X(t,\tau)=E\{X(t+\tau)X^*(t)\}$$

where $E\{\cdot\}$ denotes the expectation operator, and $^*$ denotes complex conjugation. Note that this definition with the complex conjugate is consistent with the usual definition of the inner product of complex functions.

The cross-correlation of two (possible complex) random processes $X(t)$ and $Y(t)$ is defined by

$$R_{XY}(t,\tau)=E\{X(t+\tau)Y^*(t)\}$$

Consequently, the cross-correlation of $X(t)$ and its complex conjugate $X^*(t)$ is given by

$$R_{XX^*}=E\{X(t+\tau)X(t)\}$$

If $X(t)$ is wide-sense stationary (WSS), then its autocorrelation function only depends on the time difference $\tau$ and is written as $R_X(\tau)$. Similarly, if two processes $X(t)$ and $Y(t)$ are jointly WSS, their cross-correlation function only depends on the time difference $\tau$ and is written as $R_{XY}(\tau)$.

For WSS processes, the power spectrum is defined as the Fourier transform of their autocorrelation function:

$$S_X(\omega)=\int_{-\infty}^{\infty}R_X(\tau)e^{-j\omega\tau}d\tau$$

Similarly, for two jointly WSS processes, their cross-spectral density is given by the Fourier transform of their cross-correlation function:

$$S_{XY}(\omega)=\int_{-\infty}^{\infty}R_{XY}(\tau)e^{-j\omega\tau}d\tau$$

Now for the answer of the question. To see that $f(t)$ is generally not WSS, define a complex-valued WSS process $X(t)$ by

$$X(t)=I(t)+jQ(t)\tag{1}$$

with real-valued and jointly WSS processes $I(t)$ and $Q(t)$.
The process $f(t)$ is then given by

$$f(t)=\Re\{X(t)e^{j\Omega t}\}=\frac12(X(t)e^{j\Omega t}+X^*(t)e^{-j\Omega t})\tag{2}$$

The auto-correlation function of the real-valued process $f(t)$ is

$$R_f(t,\tau)=E\{f(t+\tau)f(t)\}\tag{3}$$

Plugging (2) into (3) gives (after some straightforward manipulations)

$$R_f(t,\tau)=\frac12\Re\{R_X(\tau)e^{j\Omega\tau}\}+
\frac12\Re\{R_{XX^*}(\tau)e^{j\Omega(\tau+2t)}\}\tag{4}$$

From (4) it is clear that the auto-correlation function of $f(t)$ depends on $t$, and, consequently, the power spectral density of $f(t)$ is not defined. However, if $R_{XX^*}(\tau)=0$ and if $E\{X(t)\}=0$ then $f(t)$ is WSS and has a power spectral density.

The cross-correlation of $X(t)$ and $X^*(t)$ is given by

$$\begin{align}R_{XX^*}(\tau)&=E\{X(t+\tau)X(t)\}\\&=E\{(I(t+\tau)+jQ(t+\tau))(I(t)+jQ(t))\}\\&=\ldots\\&=R_I(\tau)-R_Q(\tau)+j(R_{IQ}(\tau)+R_{IQ}(-\tau))\end{align}\tag{5}$$

where $R_I(\tau)$ and $R_Q(\tau)$ are the auto-correlation functions of $I(t)$ and $Q(t)$, respectively, and $R_{IQ}(\tau)$ is the cross-correlation of $I(t)$ and $Q(t)$. From (5) it is obvious that the condition $R_{XX^*}(\tau)=0$ is satisfied if

$$\begin{align}R_I(\tau)&=R_Q(\tau)\\
R_{IQ}(\tau)&=-R_{IQ}(-\tau)\end{align}\tag{6}$$

is satisfied, i.e. if $I(t)$ and $Q(t)$ have the same auto-correlation function (and thus the same power spectrum), and if $R_{IQ}(\tau)$ is an odd function, the latter implying $R_{IQ}(0)=0$, i.e. $I(t)$ and $Q(t)$ are uncorrelated when sampled at the same instant.

The auto-correlation function $R_X(\tau)$ can be expressed in terms of $R_I(\tau)$, $R_Q(\tau)$, and $R_{IQ}(\tau)$:

$$\begin{align}R_X(\tau)&=E\{X(t+\tau)X^*(t)\}\\&=\ldots\\&=R_I(\tau)+R_Q(\tau)-j(R_{IQ}(\tau)-R_{IQ}(-\tau))\end{align}\tag{7}$$

If the conditions (6) are satisfied, this becomes

$$R_X(\tau)=2(R_I(\tau)-jR_{IQ}(\tau))\tag{8}$$

Pluggin (8) into (4) gives for the auto-correlation function of $f(t)$

$$R_f(\tau)=R_I(\tau)\cos(\Omega\tau)+R_{IQ}\sin(\Omega\tau)\tag{9}$$

(assuming the conditions (6) are satisfied). From (9) the power spectral density of $f(t)$ can be expressed by the power spectral densities of $I(t)$ and $Q(t)$, and by the cross-power spectral density of $I(t)$ and $Q(t)$:

$$S_f(\omega)=\frac12[S_I(\omega-\Omega)+S_I(\omega+\Omega)]+\frac{1}{2j}[S_{IQ}(\omega-\Omega)-S_{IQ}(\omega+\Omega)]\tag{10}$$

which is the final result. Note that $S_f(\omega)$ is of course real-valued because due to the second condition in (6), the cross-power spectral density $S_{IQ}(\omega)$ is purely imaginary and, consequently, the second term on the right-hand side of (10) is real-valued.


What are the mathematical tools that I can employ for automatic music composition?

There is an open source on the similar project you may be interested: http://sourceforge.net/projects/amc/
The program generated chords and bass lines from .wav or .midi input.


How to compare two speech recognition softwares?

NIST has a web page that provides a good detailed discussion of how to score ASR software.

http://www.isip.piconepress.com/projects/speech/software/tutorials/production/fundamentals/v1.0/section_04/s04_03_p01.html


Differences in log scales correspond to ratios of the underlying values.  Saying you are within .01 on a $\log_{10}$ scale means the true value is within $10^{-.01}$ and $10^{.01}$ of the stated value, ratio wise.  That calculates out to 0.9772 and 1.0233.  Since the value is so close to 1, an interval on the ratio scale is close to an interval on a difference scale, so you could also say the results are within about 2.3%.

I haven't read the paper.

Hope this helps.

Ced

Followup:

So, on a whim I looked up "centibel" and found https://en.wiktionary.org/wiki/centibel.  A "cB" if you will.  Not too many other references though.

I should have also looked at the OP's math link.  The correct answer was already there.


How to estimate a transfer function from a magnitude-only frequency response?

One approach would be to use the frequency-domain least-squares (FDLS) method. Given a set of (complex) samples of a discrete-time system's frequency response, and a filter order chosen by the designer, the FDLS method uses linear least-squares optimization to solve for the set of coefficients (which map directly to sets of poles and zeros) for the system whose frequency response matches the desired response with minimum total squared error.

The frequency response of an $N$-th order linear discrete-time system can be written as:

$$
H(\omega) = H(z)|_{z=e^{j\omega}}
$$

where $H(z)$ is the system's transfer function in the $z$ domain. This is typically written in the rational format that follows directly from the system's difference equation:

$$
H(z) = \frac{\sum_{k=0}^{N}b_kz^{-k}}{1 + \sum_{k=1}^{N}a_kz^{-k}}
$$

The frequency response is therefore:

$$
H(\omega) = \frac{\sum_{k=0}^{N}b_ke^{-jk\omega}}{1 + \sum_{k=1}^{N}a_ke^{-jk\omega}}
$$

Rearrange the above to get:

$$
\sum_{k=0}^{N}b_ke^{-jk\omega} - H(\omega) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega}\right) = 0
$$

This equation is linear in the $2N+1$ unknown system difference equation coefficients $b_k$ and $a_k$. Given a desired frequency response $H(\omega)$, we would like to find coefficients that meet the above equation exactly for all values of $\omega$. For the general case, that's hard. So instead, we will search for a set of coefficients for a system whose frequency response approximates the desired response at a discrete set of frequencies.

In order to solve for an appropriate set of coefficients using the linear least squares method, we generate an overdetermined system of equations in those unknowns. To generate those equations, choose a collection of frequencies $\omega_m \in [0, 2\pi), m = 0, 1, \ldots , M-1$ (where $M > 2N+1$, and often $M \gg 2N+1)$. For each frequency, substitute the corresponding value of $\omega_k$ into the above equation to yield:

$$
\sum_{k=0}^{N}b_ke^{-jk\omega_k} - H(\omega_k) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega_k}\right) = 0
$$

The values $H(\omega_k)$ are obtained by sampling the desired frequency response at the chosen frequencies $\omega_k$. After generating the system of linear equations, the least-squares solution for the system's coefficients $b_k$ and $a_k$ (and therefore its poles and zeros) is easily obtained. If you substitute those coefficients back into the equation for $H(\omega)$ shown above, it should (hopefully) yield a function that is close to the template frequency response that you started with.

This technique has a few advantages:

Any arbitrary complex (magnitude and phase) frequency response can be used as the template. If you only have a magnitude constraint, you could just pick a phase response, such as linear phase.
It can be used to design both FIR and IIR filters; for an FIR realization, just remove the $a_k$ coefficients from the above.
The technique is very simple to implement and is easily parameterizable based upon the desired system order.
While there may not be a good way to estimate a priori what the required system order is to meet your design constraints, it is simple to iteratively increase the order $N$ until some selected error metric is met (such as peak error, total squared error, or deviation within a specific band).

You could extend this method a bit to use weighted least-squares optimization if needed; this would allow you to specify regions of the frequency response whose approximation error is weighted more than others. This allows you to more tightly control passband/stopband areas while allowing more slop in "don't-care" areas.


Your code has some troubles, and perhaps also your theoretical understanding. Let me put here a mini summary of observing a practical sine wave and its frequency spectrum using FFT function of MATLAB / Octave.

Assume that there's a continuous-time ideal infinitely long sinusoidal wave with frequency $\Omega_0$ given as:

$$ x(t) = A \cos(\Omega_0 t) , \tag{1} $$

Let's sample this infinetely long signal within a finite duration of $[0 , T_s]$ using a sampling rate of $F_s$ or the sampling period of $T_s = \frac{1}{F_s}$, so that $t = nT_s$, and  we obtain $L$ samples of the signal :

$$ x_L[n] = x[n]w[n] = A \cos(\omega_0 n) ~~,~~n = 0,1,...,L-1 \tag{2}$$

where $\omega_0 = \Omega_0 T_s$ and $w[n]$ is a rectangular window, used to represent the finite duration of observation.

$$ w[n] = \begin{cases} { 1 ~~~,~~~~ 0 \leq n \leq L-1 \\ 0 ~~~,~~~ ow }\end{cases} \tag{3}$$

DTFT (Discrete-time Fourier transform) modulation property can be used to analyse the frequency domain representation of $x_L[n]$ :

$$
X_L(e^{j\omega}) = \frac{1}{2\pi} X(e^{j\omega}) \star W(e^{j\omega})  \tag{4}
$$

where $W(e^{j\omega})$ is the DTFT of the rectangular window :

$$
W(e^{j\omega}) = e^{-j \omega (\frac{L-1}{2}) } \cdot  \frac{ \sin(\omega L /2)}{\sin(\omega/2)} \tag{5}
$$

A plot of this function is shown below :

Again we also we know that the DTFT of $x[n] = A\cos(\omega_0 n)$ is :
$$
X(e^{j\omega}) = A\pi \delta(\omega + \omega_0) + A\pi \delta(\omega - \omega_0)  \tag{6}
$$

Hence we see that the DTFT of the finite length observation $x_L[n]$ is

$$
X_L(e^{j\omega}) = \frac{A}{2} [ W(e^{j(\omega+\omega_0)}) +  W(e^{j(\omega-\omega_0)})] \tag{7}
$$

Now, when you take $M$-point FFT of the $L$-point sequence $x_L[n]$, then you are effectively calculating the $M$-point DFT (discrete Fourier transform) which is equal to $M$ uniform samples of DTFT $X_L(e^{j\omega})$  as given by

$$
X_L^M[k] = X_L(e^{j\omega})|_{w = \frac{2\pi}{M} k} ~~~, ~~~k = 0,1,...,M-1 \tag{8}
$$

which yields the complex frequency vector that you compute using FFT function:

$$
X_L^M[k] = \frac{A}{2} [ W(e^{j(\omega_k+\omega_0)}) +  W(e^{j(\omega_k-\omega_0)})]  ~~,~~k = 0,...,M-1 \tag{9}
$$

So with an M-point FFT of $x_L[n]$ we calculate the samples of the shifted versions of $W(e^{j\omega})$.

An example plot of $|X_L^M[k]|$ for $L =  16$ , $M = 1024$, $A = 5$ and $\omega_0 = 0.2 \pi$ is shown.

Note that the peak values (which are $A \cdot L /2$ when $\omega_0$ meets some criterion, see below) are obtained for some $k=k_0$ and $k = M-k_0$. Consider the first peak (to the right of the plot here) and it can be shown that its value is

$$ |X_L^M[k_0]| = |0.5 A [  W(  \frac{2\pi}{M} k_0 - w_0) +  W(  \frac{2\pi}{M} k_0 + w_0) ]| \tag{10}
$$

Where the integer DFT index $k_0$ is found to be $k_0 = \text{round}( \frac{M \omega_0}{2\pi} )$. So the peak value of the spectrum is in general not equal to $0.5 A \cdot L$ which you are expecting, which would only be the case when

$$  \frac{2\pi}{M} k_0  =  w_0 \tag{11}$$ for some integer $k_0$ which requires that $\frac{M \omega_0}{2\pi}$ is also an integer, or equivalently

$$ \omega_0 = \frac{ 2\pi }{M}k_0 \tag{12}$$ ,in that special case $|X_L^M[k_0]| = 0.5 A L$ will happen. Otherwise, you will have two complex numbers added and their magnitude is computed.

Note that by selecting $M = L$ when $\omega_0$ meets the criterion (i.e, $ \frac{M \omega_0}{2\pi}$ is an integer), then one can get what is an illusion of lollipop plot indicating only two nonzero samples at the frequencies of $\pm \omega_0$ with magnitude 0.5 $A \cdot L $ and all zeros. But this is a conseqeunce of frequency sampling and not to be taken unnecessarily into too far conclusions.

Here is a MATLAB / OCTAVE code to produce the results

clc; clear all; close all

% Define the sampling grid
Fs = 1000;             % sampling frequency in Hz
Ts = 1/Fs;             % sampling period in seconds

ti = 0;                % initial time (s)
tf = .15;               % final time   (s)
t =  ti :Ts : tf;  % sampling time index t = n*Ts

L = length(t);          % number of samples in x[n] = A*sin(2*pi*f0*n*Ts)
M = L;                  % Number of DFT points to be computed

% Define the sampled signal
A = 0.5;               % amplitude of sine wave

% Signal Frequency Type-1 :
%f0 = 100;              % sine wave frequency in Hz.
%w0 = 2*pi*f0*Ts;        % set an arbitrary frequency
%x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.

% Signal Frequency Type-2:
w0 = 2*pi/M*floor(L/3); % set a special frequency so that DFT PEAK = A/2
x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.

% Calculate DFT of x[n] using FFT function
k0 = ( M*w0/(2*pi));

figure,stem(linspace(-pi,pi,M),fftshift(abs(fft(x,M))/L))
title(['Magnitude of X_L^M[k] for L=',num2str(L),' , M=', num2str(M),...
      ' , A=',num2str(A),' , \omega_0 =', num2str(w0), ' , k_0 = ' num2str(k0)]);
grid on


Basics of frequency spectrum for onset detection using FFT

Here are some thoughts concerning your questions:

I guess you have to find the relevant frequency range by experiment, because, depending on the instrument, a lot of signal energy can be in the harmonics.
If your sampling frequency is 44.1 kHz and you use a 1024 point FFT your frequency interval is about 43 Hz, because the 513 points you mentioned are between DC and half the sampling frequency.
You could indeed average some frequency bins, or simply use a shorter FFT length.
I would say your methodology is pretty reasonable. Have you checked the literature? You're probably not the first one the do onset detection.


DSP processor: some beginner's questions

I'v have some experience with TI C55 and C6000 series DSPs, I'll give you some points, hope to be helpful.

  Programming each chip that came from a different company requires different softwares or big differences in commands?

Absolutely! Although the code body might be the same, but since every company has its special processor architecture, their compilers would be different. Even two DSPs of a same company might be too different, for example Texas Instruments C5000 processors are too different from C6000, since their memory architecture is different. I mean, you can't simply use the same segment of code for both.

  Is there a company making both DSPs and Bluetooth?

Usually most of the DSPs are only Microprocessor and a few essential Interfacing (like I2C, SPI,...) are considered for them, however some are more sophisticated System on a Chip designs might include Bluetooth. Basically you can buy and add a Bluetooth module to any DSP you want.
https://www.sparkfun.com/products/12574

  If you have some tips regarding the chips I will be happy to read them :) (easiest way to begin, where I should go into the architecture and what is less importent and can look at it from a bigger picture)

It depends on your programming level and those calculations you've mentioned above, however, I think the best chip to begin with is MicroChip DSPIC, which is simply a PIC micro associated with some DSP blocks (like circular conv, FFT ,etc), I think your last choice should be TI's DSPs or AD's Shark DSPs, they are too complicated and frustrating to begin with.


RADAR Video Signal and Matched Filtering

The simple answer is the matched filter increases the signal to noise ratio of the reflected return signal.

As the signal propagates from transmitter, the power that impinges on the target is proportional to $1/r_{target}^2$.  The reflected energy back to the transmitter is also proportional to $1/r^2$, so the received energy is proportional to $1/r^4$.

The received back echo power is very small and the matched filter maximizes the deflection of the echo with respect to the noise in the receiver.

This is a simple explanation. There are many types of RADAR, and there are additional ways to make the echo stand out with respect to the noise of the system and the clutter echos present in all real systems

Those additional ways that are specific to the different kinds of RADAR systems to make the echo discernible. There are various ways to modify a matched filter for a specific kind of applications but it would be hard to say they were alternative to a matched filter. The underlying optimality characteristic of a matched filter justifies few alternatives.


What is the function of init in Active Contour Model

Active contour is just refining your initial contour.
You could play yourself with the model :)


The only quick reference I can find is this one, which states:


What are the differences between notch filter & Goertzel?

The Goertzel algorithm allows you to sample the DTFT of a signal with slightly lower complexity than direct computation (still $\mathcal{O}(N^2)$ to produce the DFT but with a smaller coefficient). I think of it as being advantageous in two ways:

If you know the frequency you are looking for (like in telephone dialing systems)
If you only want a small number of frequencies so that computing the FFT and discarding those you don't want is overly computational

Of course if the frequency you are looking for is not a multiple of $\frac{2\pi}{N}$ where $N$ is the signal support, then the FFT won't get it anyways (sans zero padding tricks).

Notch filters, are, well, filters. So their complexity is straightforward. An advantage in some cases is that they can have a non-impulse frequency response (but still sharp over a small concentrated band) so that if you are only sure of a specific region the signal you are trying to detect is in, this is a good way to do some "averaging".

As an interesting aside, Goertzel's original paper doesn't mention directly applications to Fourier analysis (but rather just people who wanted to program trigonometric series), and is less than a single page long)


Phase spectrum: $\mathrm{arg}\mathcal{F}(\omega)$ or $\mathrm{-arg}\mathcal{F}(\omega)$?

Both definitions are common. However, if you talk about the phase, then what is usually meant is the (frequency dependent) phase angle $\phi(\omega)$ of the Fourier transform $F(\omega)$:

$$F(\omega)=|F(\omega)|e^{j\phi(\omega)}\tag{1}$$

The negative phase $\beta(\omega)=-\phi(\omega)$ is often called the phase shift or phase lag. It is commonly used together with the attenuation $\alpha(\omega)$, which is the negative natural logarithm of $|F(\omega)|$:

$$\alpha(\omega)=-\log (|F(\omega)|)\tag{2}$$

$$F(\omega)=e^{-(\alpha(\omega)+j\beta(\omega))}\tag{3}$$


Spectral images - How do they appear during ADC?

One way to represent AD conversion mathematically is to think of it as a multiplication of the signal by an impulse train, which is a sequence of dirac delta distributions spaced by the sampling period $T_s$.

With this in mind and remembering that multiplicating in the time domain is equivalent to doing convolution in the frequency domain, it's possible to see that the frequency representation of the digitalized signal will be the convolution of the frequency representation of the sampled analog signal with the frequency representation of the impulse train.

As it turns out, the frequency representation of an impulse train spaced by $T_s$ is an impulse train spaced by the sampling frequency $f_s$. Also, the convolution of a dirac delta with any function $g$ is just $g$ itself. By displacing the delta function by $f_s$, we also displace the result of the convolution by $f_s$.

So, the spectrum of the digitalized signal will be an infinite sum of copies of the spectrum of the analog signal distant of each other by $f_s$. Since the spectrum is bilateral, if your analog signal has frequency components above $f_s/2$, those components will end up summed from two different replicas in the spectrum of the digital signal, and this is in fact the aliasing effect.