Position of poles and Stability in $z$ domain

Short Answer:
All the poles of a causal (right-sided) and stable LTI system must be inside the unit circle whereas all the poles of an acausal (left-sided) and stable LTI system must be outside the unit circle.
Explanation:
First, consider the following right-sided (causal) exponential sequence $h[n] = a^n u[n]$ whose $\mathcal{Z}$-Transform is:
$$H(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - a z^{-1}} ~~~,~~~ |z| > |a|.$$
As can be seen for a right-sided (causal) sequence the ROC of $H(z)$ extends outward from the pole to the infinity.
Now, consider a causal LTI system with right-sided impulse response $h[n]$ whose transfer function can be expressed in terms of ratio of polynomials:
$$H(z) = \sum_{n=0}^{\infty} h[n]z^{-n}  = \frac{P(z)}{Q(z)} = K \frac{\prod_{k=0}^{M-1} (1 - d_z z^{-1})}{\prod_{k=0}^{N-1} (1 - d_p z^{-1})} $$
Where $d_p$ are the poles of the transfer function at which $H(z)$ diverges. Hence, associated with the $\mathcal{Z}$-Transform $H(z)$, is a region of convergence (ROC) which is singly connected region that does not include any poles inside. Which implies that ROC for the causal system extends outwards from the largest pole.
On the other hand, the BIBO stability of LTI system requires that $|y[n]|<\infty$ for all bounded inputs $|x[n]| < B$, which translates to $h[n]$ be absolutely summable as follows:
$$ |y[n] = \sum_{k} h[k] x[n-k] | < \infty \implies \sum_{k=0}^{\infty} |h[k]| <  \infty$$
This means, for stable system, at $|z| = 1$, the $\mathcal{Z}$ transform
$$ H(z) = \sum_{k=0}^{\infty} h[k] z^{-k} < \sum_{k=0}^{\infty} |h[k] z^{-k}| = \sum_{k=0}^{\infty} |h[k]|  <  \infty$$
converges! Thus for a BIBO stable LTI system, its ROC must include the unit-circle $|z|=1$ inside.
Combinig these two concepts yields: for a causal LTI system to be stable its ROC must include the unit circle and it should extend outwards  from the largest pole to infinity. The corollary of this fact is that the largest pole of the stable and causal LTI system is inside the unit circle. Hence we conclude that all of the poles of the stable and causal systems are inside the unit circle.
An opposite complementary conclusion occurs for the stable but acausal LTI systems (whose $h[n]$ includes left-sided terms like $-a^nu[-n-1]$) for which the ROC is inwards from the smallest pole, and must include the unit circle (for stability), and therefore that all of its poles must be out of the unit circle; i.e., it's smallest pole is larger than $|z|=1$.


What is the meaning of the delay in a FIR filter

In mathematical terms, time discrete systems are are most often expressed in terms of difference equations. A delay is simply $y[n]=x[n-1]$ so the difference equation for your FIR filter would be $$y[n]=\sum_{i=0}^{M}b_{i}\cdot x[n-i]$$.

See https://ccrma.stanford.edu/~jos/filters/ for a good primer


wavelet analysis normalized to have unit energy

Take the mean of the original signal and add it to the reconstructed signal. IIRC wavelets have no DC response, and so the reconstructed signal will have 0 mean.


Yes, you should be able to see the component of the piano note. But there are quite a few pitfalls on the way. I would suggest you experiment a little to get a 'feel' on how fft works first. Use GNU Octave or Matlab to do simple experiments, and then you can plug in your real data, and see if it works. For example, in Octave I can easily generate a sine wave of 440 Hz and do the FFT:

t = linspace(0, 2*pi, 44100);
s = sin(440.0 * t)              % This is the 440 Hz tone
f = abs(fft(s))                 % Take the fft (magnitude)
plot(f)
xlim([0,1000])                  % Show only the 1000 first bins

Then you can read your samples and do the same. At least, this way you
exclude insecurity about the fft function you're using. You could also
save the generated sine, and feed it to your four1 function, and check
if the result is correct.


These filter types are well covered here http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

However these are all single biquad IIR implementations. To turn this into an FIR, you could simply calculate the impulse response of these filters and window to a finite size with the desired amount of accuracy. Another option is to design the filter as IIR, sample the frequency response and then do any of the many FIR fitting techniques.

In general you will find that FIR filters are not well suited for this and the number of taps required will heavily depend on the relationship between the corner frequency and the sample rate and also on the Q of the filter.


How do I deepen my voice (I'm a female) to sound like a man in a video for a school project?

I know some audio software that does a good job:

Audacity (free software) - effect -> change pitch (set the pitch to down, lower down the semi-tones and the percent change)
Adobe Audition (paid software) - there is a plugin called pitch shifter (lower down a little bit the semi-tones and cents)

Both are audio software which means you will need to export your audio from the video clip and because there are other people talking in the video you will need to select the parts only were you are talking and then apply the effect.


Need help to understanding concept of particle filter

The following papers are good resources:

Gordon, N. J.; Salmond, D. J. and Smith, A. F. M. (1993). "Novel approach to nonlinear/non-Gaussian Bayesian state estimation". IEE Proceedings F on Radar and Signal Processing 140 (2): 107–113.
Arulampalam, M.S.; Maskell, S.; Gordon, N.; Clapp, T.; (2002). "A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking". IEEE Transactions on Signal Processing 50 (2): 174–188.
Doucet, A.; Johansen, A.M.; (December 2008). "A tutorial on particle filtering and smoothing: fifteen years later". Technical report (Department of Statistics, University of British Columbia). (Available here)

I'd recommend you read at least the first one.

A particle filter works on a model specified by 2 equations (a hidden Markov model):

1) The state equation
$x_{k+1} = f_k (x_k,w_k)$ which tells you how the state (minimal set of variables to describe a system) evolves as being driven by noise.

2) The measurement equation $y_k = h_k (x_k, v_k)$ which maps the state to what you observe.

Given the observations $y_0, \ldots, y_n$, estimate the distribution of $x_n$ is the problem the particle filter tries to solve. Once you have the approximated distribution of $x_n$ (which is the distribution of the particles), you can estimate $x_n$ (say, by averaging the particles).

The particle filter essentially starts with a bunch of samples (called particles), evolves the state by running each particle through the state equation and re-samples the particles based on the observation you see in order to make the distribution of particles consistent with the observations.


Fractional Frequency and negative Frequencies

On figures 1, 2 and 3, you are plotting 2 curves which you call real and imaginary parts. These can be interpreted as the complex exponential

$$
x[n] = e^{jn\theta_0}
$$

with $\theta_0$ being 0, $2\pi/8$ and $2\pi/2$. These are not cosines, but complex exponentials.

Now, on figure 4, you only draw the sine part (what you call Imaginary), and not the cosine part. So this is why you don't find a cosine part!! You forgot to draw it. It has nothing to do with aliasing (which you interpret correctly).


There a few commercial algorithms that do exactly that (Dolby Prologic, DTS Neo 6, Lexicon Logic 7, Bose Videostage, etc.). If your project has the funds, you can simply try to license one of those.

The inner workings of these algorithms are rather complicated and typically a mixture of time domain and frequency domain feature extraction, some steering logic, and some real time steering mechanics (time variant mixers, sub band filters, etc.)

Simple signals like your example are normally easy to tease apart, however almost all music is  partially correlated and overlapping both in time and frequency. For starters it's difficult to define what actually SHOULD be in the center and the Left and Right. The algorithm has to make a lot of "artistic" decision to make sure there are no steering artifacts and it meets the design intent over a wide variety of source material.


This reference provides specific implementations for three error correcting BCH generators with a minimum hamming distance of 7

BCH codes

Note, also refer to Peterson's Table of Irreducible Polynomials which lists all primitive and irreducible polynomials over GF(2):


Rajeev Bansal in the appendix to this article explains the $2D/\lambda$ rule:

How does one arrive at the various criteria for the far-field zone?
Basically, the criteria are guidelines to the boundaries where the
fields start to approximate the “ideal” assumed characteristics. First
of all, only the radiation ($1/r$) terms remain significant; higher
order terms fade away. Second, in the far-field zone, the angular
field distribution becomes independent of the distance [5]. Third,
only transverse field components remain, and the ratio of the electric
and magnetic field components approaches the free-space impedance, 377
ohms [6 ,7]. Finally, for a receiving antenna, the incoming wave-front
is nearly planar across the aperture. In fact, the $r = 2 D^2/\lambda$
formula corresponds to a phase error (due to the curvature of the
actual spherical wave-front) of no more than $\pm 22.5$ degrees across the
aperture, as compared with the ideal plane wave-front [6, 8].
References:

IEEE Standard Dictionary of Electrical and Electronics Terms, IEEE standard 100-1984, 1984.

C. Paul, K. Whites, and S. Nasar, Introduction to Electromagnetic Fields, 3rd ed., McGraw-Hill, 1998.

Reference Data for Radio Engineers, 6th ed., SAMS, 1981.

G. Smith, An Introduction to Classical Electromagnetic Radiation, Cambridge, 1997.


In principal, you get $\text{Ws}=\text{J}$ for energy and $\text{W}$ for power, but as you are calculating these measures from a wav-file, you do not get any unit. What you get is just a ratio to the maximum possible value of a wav-file. A square wave with amplitude one has a level (power) of $0\text{ dBFs}$, which means dB Fullscale. Your signal has roughly $-75\text{ dBFs}$, so that would be a valid result. If you want to get a result in actual physical units, you will need to know several things:

the sensitivity of the sensor (e.g. a microphone), that was used for audio recording. This is usually given in $\text{dBV/Pa}$
the calibration of the A/D converter used, usually given in form of the overload point in $\text{dBV}$, meaning that a sine wave with this level just drives the ADC to fullscale ($0\text{ dBFs}$), any increasing of the voltage leading to clipping.

With this information, you could calculate back the real power the actual sound event had at the point of the microphone.


Algorithm for detecting whether a signal stored in rolling buffer is Increasing, Decreasing or staying Constant

What sort of algorithm would effectively notice this?

What you are describing, is an online algorithm that uses a small set of "past" measurements to make an inference, in this case, about the slope of the signal.

Since you are keeping track of the current beginning and end of the rolling buffer, it is essentially like having access to a linear buffer of N past samples.

So, provided that the phenomenon you are trying to study (your signal) evolves slower than the length of the buffer you have chosen, you can apply a plain simple slope estimation on the values of the buffer and then use the slope estimate value (and its residuals) to decide if the signal is mostly heading downwards or upwards or remains constant. The most simple "classifier" here, would be a plain simple threshold one, possibly with the addition of hysteresis as it has been suggested already. More complex "classifiers" could be monitoring both the first and second derivative and returning additional information on whether the signal also changed direction within a buffer.

If you want to be getting one new estimate every time a new sample is received, you would have to look for iterative (or "recursive") methods, (for example, iterative least squares).

Hope this helps.


why does the spectral envelope of human speech not change w.r.t. pitch when taking a Fourier transform?

The spectral envelope, which determines where the formant frequencies are, is determined by the shape of the mouth, tongue, lips, and nasal coupling.  That is independent of the pitch, which is dependent on the tension or stiffness of the vocal cords.


Solution verification for this Fourier series problem

You made a small mistake in the final step. The second term in parentheses should be $2((-1)^k-1)$ because you have two terms $e^{-j\pi k}$. So your final result should be
$$X[k]=\frac{3}{2\pi jk}\big(1-(-1)^k\big)\tag{1}$$
With $(-1)^k=e^{-j\pi k}$, Eq. $(1)$ can be rewritten as
$$\begin{align}X[k]&=\frac{3}{k\pi}\frac{1-e^{-jk\pi}}{2j}\\&=\frac{3e^{-jk\pi /2}}{k\pi}\frac{e^{jk\pi /2}-e^{-jk\pi /2}}{2j}\\&=\frac{3e^{-jk\pi /2}\sin(k\pi/2)}{k\pi}\tag{2}\end{align}$$


Sampling with an alternating phase offset

Bottom Line
Given a sampled tone at normalized radian frequency $\omega_1$ and real amplitude $A$:
$$x[n] = A e^{j\omega_1n} \label{1} \tag{1}$$
If every other sample is shifted by a constant phase $\phi$, the result in the Fourier Transform will be two tones $x_1$ and $x_2$ (precisely so as the only non-zero values when the sampling rate is an integer multiple of $\omega_1$, otherwise $x_1$ and $x_2$ refer to the maximum bin in vicinity of the two peaks):
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
And $\phi$ is determined using the magnitude and phase of each of the two bins as follows:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$

Further Details
The first tone will be at the bin associated with the same normalized radian frequency $\omega_1$ as the original but shifted in phase by $\phi/2$ and an additional tone will be given due to the phase modulation at the bin associated with the normalized radian frequency $\pi + \omega_1$ and initial phase such that the new tone will always be in quadrature phase with $x[n]$ as $-j (-1)^n x[n]$. The amplitude of the original tone will be reduced according to $A\cos(\phi/2)$ while the amplitude of the additional tone will be $A\sin(\phi/2)$. Given the modified original tone as $x_1[n]$, and the new tone as $x_2[n]$ this relationship to $x[n]$ ahd $\phi$ is:
$$x_1[n] = x[n]\cos(\phi/2)e^{j\phi/2}= A\cos(\phi/2) e^{j(\omega_1n + \phi/2)}\label{2} \tag{2}$$
$$x_2[n] = -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2} = jA\sin(\phi/2)e^{j((\pi+\omega_1)n+ \phi/2)}\label{3} \tag{3}$$
So given this and not knowing $\phi$, we can determine $\phi$ from the DFT results as follows: Given $x_1$ and $x_2$ as the complex values for two bins where $x_1$ is the frequency bin for the expected tone and $x_2$ is the frequency bin for the interference tone due to the introduced distortion of an unknown fixed phase shift every other sample, with:
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
$x_1$ and $x_2$ will always be in quadrature, and the phase between the $x_1$ and $x_1+x_2$ will always be $\phi/2$. Further, when $x[n]$ has no phase offset as we had introduced in \ref{1}, then the sum  of $x_1+x_2$ will also always be real and the phase can be determined from twice the phase angle of $x_1$ alone. However, what is developed here is the more general case applicable to any starting phase offset for $x[n]$. This general case can be solved for using the following:
$$\phi = 2 \times \text{atan2}(|x_2|\sin(\phi_1-\phi_2), |x_1|)$$
Where $\text{atan2}$ is the 2-element arctangent. The term $\sin(\phi_1-\phi_2)$ will always be either +1 or -1 depending if $\phi$ is positive or negative, since the relationship of the two tones will always be in positive or negative quadrature, so the above can also be written as:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$
To extend this to a real waveform such as in the case of the OP, the original tone from the positive frequency axis would be used as $x_1$, and it's negative image would be used as $x_2$ as depicted in the graphic below:

Using a real waveform that is exactly at normalized frequency $\omega_1 = \pi/2$ would not be resolvable since the image from the added phase distortion and the primary negative frequency  component would both be on the same frequency location.

Derivation:
Relating \ref{1},  \ref{2} and \ref{3}:
$$x[n] = x_1[n] + x_2[n]$$
$$ = x[n]\cos(\phi/2)e^{j\phi/2} -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2}$$
$$ = x[n]e^{j\phi/2}(\cos(\phi/2) - j(-1)^n\sin(\phi/2))$$
$$ = x[n]e^{j\phi/2}e^{-j\phi/2(-1)^n}$$
$$ = x[n] \text{ for n even}$$
$$ = x[n]e^{j\phi} \text{ for n odd}$$
Intuition:
The above relationships are more intuitive to understand knowing that the Fourier Transform decomposes any single valued time domain analytic function into distinct phasors (as the basis functions) that each must rotate with time with a constant magnitude and at a constant rate. We can see the solution for how this would occur for the OP's problem by doing the decomposition geometrically.
First consider an arbitrary $x[n]$ as a rotating phasor on the complex plane in the time domain as shown on the left, where when sampled the normalized radian frequency corresponds to the phase shift from sample to sample. The Fourier Transform of this on the right is a single impulse in frequency, a positive frequency here given the counter-clockwise rotation of the phasor (positive increasing angle). The single impulse corresponds directly to the single spinning phasor.
x[n] with no distortion:

Now consider if every other sample was shifted forward in phase by a constant amount $\phi$:
x[n] with phase distortion:

To decompose this into phasors of constant magnitude and rate of rotation, the first phasor MUST be the first phasor rotated by $\phi/2$ since the error represents a modulation going back and forth by $\phi/2$ for every sample. And for this same reason the second phasor must be in alternating quadrature relationship such that it is modulating the phase back and forth from sample to sample (resulting in a normalize angular frequency of $\omega_1 - \pi$ when we consider it rotating opposite the original tone, or similarly within the range of $0$ to $2\pi$ at the normalized angular frequency of $\pi+\omega_1$.)
Fourier Decomposition of distorted x[n] into two tones:

And if we look at the relationship between any two samples at a given time instance in the decomposed waveform of two phasors, we see through the geometric relationship how the two must be always in quadrature with the cosine and sine relationship of magnitudes such that the normalized magnitude of the sum of the two phasors will be normalized (to one here) as in the original waveform. Consider the first two samples of the original waveform with the second sample offset by an additional $\phi$ as shown by the red phasors below, together with the first two samples of the decomposed tones showing how the quadrature relationship and relative magnitudes must be maintained:
First Sample of distorted x[n] and the Decomposed Tones:

Second Sample of distorted x[n] and the Decomposed Tones:

In summary, we shift the original tone by $\phi/2$ so that each sample it is midway between the sum of the two decomposed phasors. Then in order for the phase to jump back and forth by exactly $\phi/2$ from sample to sample, the second tone must be in alternate quadrature relationship from sample to sample (so will have rotated $-\pi + \omega_1$ since the original phasor modified with a phase shift of $\phi/2$ still rotates at the same rate as $+\omega_1$). Further for the sum of the two decomposed phasors to be on the unit circle (or whatever magnitude the original tone was), the magnitudes of each are related by the $\cos(\phi/2)$ and $\sin(\phi/2)$ since $\cos(\phi/2) + j\sin(\phi/2) = 1$.
This shows the relationship between the Fourier Transform components in the general case with a single complex tone, and taking the real of the above would result in two symmetric tones with the same relationships as is the case for the OP.  Further, to be precise, the descriptions above are specific to the Discrete Time Fourier Transform (DTFT) in which the time domain waveform is extended for all time and would be a continuous function in frequency with the exact solution shown as 2 weighted impulses in frequency corresponding to the two phasors in the time domain and zero everywhere else. To extend this derivation approach to the DFT, x[n] would be truncated in time (windowed with a rectangular window of duration $T$) which would result in convolving the impulses in frequency with a Dirichlet Kernel (aliased Sinc function) and then sampling the resulting frequency response at a spacing given by $1/T$ where $T$ is the window duration. When the original tone is on a bin center of the DFT (multiple of $1/T$), then the DFT would also result in just two tones for this case. For all other cases, multiple smaller tones would also appear (spectral leakage) but the same relationship would apply between the strongest tone closest to $\omega_1$ and the second tone at $\pi+\omega_1$, where the $N$ frequency bins given as integers $k$ from $0$ to $N-1$ are related to the normalized radian frequency $\omega$ as $\omega = 2 \pi k/ N$ .


Clarification on Discrete-time convolution and signal shifting

You're right that both results should be identical. The mistake lies in your conclusion that $y_1[n]=\infty$ in the first case, that's not true. You get

$$y_1[n]=(0.1)^nu[n]\sum_{m=0}^n3^m=(0.1)^nu[n]\frac{3^{n+1}-1}{3-1}=\ldots\tag{1}$$

Note that you must also include the unit step $u[n]$ in $(1)$ because the response is zero for $n<0$.


Check out the comment in this code:

// ----- Parameters and Variables for Tone Detector -----
// cos = (cos (2*PI*(DTMF_Freq/8000.0))) * 256*128;
#define DTMF_697Hz   27980        // DTMF Row Frequency
#define DTMF_770Hz   26956
#define DTMF_852Hz   25701
#define DTMF_941Hz   24219
#define DTMF_1209Hz  19073        // DTMF Column Frequency
#define DTMF_1336Hz  16325
#define DTMF_1477Hz  13085
#define DTMF_1633Hz   9315
#define DTMF_1394Hz  15014        // DTMF Row Frequency    2nd harm
#define DTMF_1540Hz  11583
#define DTMF_1704Hz   7549
#define DTMF_1882Hz   3032
#define DTMF_2418Hz -10565        // DTMF Column Frequency 2nd harm
#define DTMF_2672Hz -16503
#define DTMF_2954Hz -22318
#define DTMF_3266Hz -27472

It explains how the cos_fact is being calculated pretty directly. For example, for the first item in the list:

$$
\cos( 2 \pi \times 697 / 8000) \times 256 \times 128 = 27979.5762684 \approx 27980
$$


There are many intepretation of OFDM to answer your question. One way I find intuitive is looking at the resolution of (Discrete) Fourier Transform.

By capturing a signal during time $T$, and by assuming periodicity outside $T$, DTFT is able to provide a frequency resolution $1/T$. Thus if you choose a subcarrier spacing $\Delta f$, you will want a frequency resolution at least $\Delta f$. As a consequence, signal must be captured at least $T = 1/\Delta f$.

One may argue that I capture then $T > 1/\Delta f$. It is fine, but the extra capturing time may catch the next symbol and cause ISI. Assuming that in a good system design, e.g. with CP or zero guard time, so that ISI is avoided, the resolution of DTFT (and DFT/IFFT) $1/T < \Delta f$. The issue here is that if $\Delta f$ is not a multiple of $1/T$, the output of IFFT will be at frequencies $k/T, k \in \mathbb{Z}$ outside the subcarrier positions and this output will be sinc-interpolation of data.

So, at the end the requirement is $\Delta f=k/T, k \in \mathbb{N_+}$ or $T=k/\Delta f$. For the sake of economy, we choose the smallest $T$ possible at $k=1$ yeilding the famous thumb rule "subcarrier spacing is reciprocal of the symbol period".


Estimating Convolution Kernel from Input and Output Images

Take Fourier Transform of both original image and blurred image, then divide Fourier transform of blurred image by Fourier transform of original image. This will give you the Fourier transform of kernel (except at frequencies with zero value which lead to division by zero). Then take inverse Fourier transform and find the kernel in spatial domain.


Preparations to implement FFT on a wav file (code Java, not Library)

Simply taking the Fourier transform of a ten minute signal is unlikely to provide you with any useful information - I suspect what is actually of interest to you is the spectral information at intervals throughout the signal, e.g. 600 frequency plots at 1 second intervals, which if you so desire can be combined into a spectographic image for analysis. This can be achieved by performing multiple windowed Short Time Fourier Transforms (STFT) and storing the results.

Can I to split my file in order to make smaller [calculations]?

Yes, it is up to you what granularity of frequency information you wish to obtain. Small FFT lengths will afford you a better time resolution, longer FFT lengths will afford you a better frequency resolution. For purely visualisation purposes an FFT length of 1024 is generally sufficient.

Can I to fill with zeros to reach Pow(2, Integer) [for my FFT length]?

Yes, if you do not have sufficient samples to reach a power of two length for your FFT then you should zero pad the input. This will introduce the minimal amount of spectral error into your frequency domain results.


Advantages and disadvantages of multiplying sparse vector with sparse matrix

You are referring here to sparse signal sampling and Reconstruction. You do not require the sampling matrix to be sparse. Infact we require the sampling matrix to follow the restrictive isometry property. And sparse matrices do not follow this property. Random matrices do. So no point having a sparse sampling matrix


Does scaling property imply superposition?

The answer is "no" in general. Let us see the details.

The notion of linearity goes beyond systems. Let us start from vector, or linear spaces. They deal with two objects: vectors and "scalars". In a vector space, one can add vectors, and multiply a vector by a scalar. More formally, a vector space over a field $F$ is a set $V$ together with two operations, vector addition and scalar multiplication. Those operations in $V$ satisfy a number of axioms. The field $F$ possesses two operations as well, satisfying different axioms, namely field axioms, ruling the addition or product of scalars. Those two operations are often assimilated to addition and multiplication. But in general the addition/multiplication in $F$ are not the same as the vector addition/vector scaling in $V$. Indeed, if $f \in F$, $v\in V$, then $f.v \in V$, but is not anymore an element of $F$. The scalar has been cast to a vector. This plays a role in the standard confusion in defining linearity.

Now take for instance the field of rationals $F=\mathbb{Q}$ with the two standard operations. Take $V=\mathbb{R}$ the vector space of reals over $F$ with the standard real addition and real product. To relate this to the introduction, for instance, $1.\sqrt{2}$ is not rational anymore, although the multiplication seems to be the same, because you mix a scalar ($1$) and a vector $(\sqrt{2})$.

Suppose you have a system $S$ that provides the following output: if $v_1\in V$ is rational, $S(v_1) = v_1$; if $v_2\in V$ is irrational, $S(v_2) = -v_2$. For any rational $q$, $q.v_1$ is rational, and  $q.v_2$ is irrational. Resultingly, $S(q.v_1) = q.v_1 = q.S(v_1)$, and  $S(q.v_2) = -q.v_2 = q.S(v_2)$. Hence, $s(q.v) = q.S(v)$ for any $v\in V$, so $S$ verifies the scaling property.

However, $v_1+v_2 $ is irrational. Hence, $S(v_1+v_2) = -v_1 -v_2$, which is different in general from $S(v_1)+S(v_2) = v_1 - v_2$, take $v_1 = 1$ and $v_2 = \sqrt{2}$ for instance.

So scaling does not imply superposition (in your sense) in general.

But there exists somehow converse statements. In other domains, one sometimes calls the scaling "homogeneity", and with additivity we get  the superposition principle for a system $S$: $S(q_1.v_1+q_2.v_2) = S(q_1.v_1)+S(q_2.v_2)$. Keeping with standard fields, it can be shown that additivity implies homogeneity with rational scalars, by playing on $S(v+v) = S(v)+S(v) = 2.S(v)$ etc. To go beyond rational scalar often requires additional assumptions like continuity, which could be troublesome depending on the sets you choose, see for instance Why doesn't superposition imply linearity? Why is homogeneity needed?


Before I try to answer the question, let's see how the question relates to dispersion.

The solution space of a dispersive (linear) wave equation in $n$ dimensions is spanned by the basis functions

$$f_{\mathbf{k}}(\mathbf{x},t) = \exp( i\, \omega(\mathbf{k}) t )\exp( i\, \mathbf{k}\cdot \mathbf{x} )$$

where $\mathbf{k} \in \mathbb{R}^n$ is the spatial frequency vector or wave vector. The temporal frequency $\omega$ depends on the  wave vector, or just on its magnitude $|\mathbf{k}|$ if the dispersion is isotropic. If $\omega \propto |\mathbf{k}|$ we have the special case of no dispersion. We can also assume in general that $\omega\geq 0$ to simplify calculations.

The basis functions are orthonormal in the Dirac sense with respect to their spatial coordinate:

$$\langle f_{\mathbf{k}},f_{\mathbf{k'}}\rangle := \int_{\mathbb{R}^n} f_{\mathbf{k}}(\mathbf{x},t)\, f_{\mathbf{k'}}(\mathbf{x},t)\, d\mathbf{x} = \delta(\mathbf{k}-\mathbf{k'})$$

And the general solution of the wave equation is:

$$f(\mathbf{x},t) = \int_{\mathbb{R}^n} c(\mathbf{k})\, f_{\mathbf{k}}(\mathbf{x},t)\, d\mathbf{k} $$

If we now observe that general solution at $\mathbf{x}=0$, we get the relationship

$$f(0,t) = \Re\left(\int_{\mathbb{R}^n} c(\mathbf{k})\,\exp( i\, \omega(\mathbf{k}) t )\, d\mathbf{k}^n \right)$$

which looks similar to a Fourier transform of $c(\mathbf{k})$. If we can turn this into a proper Fourier transform, the inverse Fourier transform of the observation will give us the expansion coefficients $c(\mathbf{k})$ required for the reconstruction of the full spatial solution of the wave.

At this point, we need to make a few additional assumptions. The first one is that $n$ equals 1 so that the coefficients are also a function on $\mathbb{R}$. We also have to assume that the dispersion is isotropic and $\omega(|\mathbf{k}|)$ is differentiable and strictly monotonic.

This is not enough however. We also need to assume a direction of travel so that $\omega(|\mathbf{k}|)$ becomes an invertible function on a restriction of $k>0$ or $k<0$ respectively.

We pick a "right" travelling wave, so $k>0$. With this, we can change the integration variable with $dk=\frac{\partial k}{\partial\omega}\,d\omega$ to get

$$f(0,t) = \Re\left(\int_\mathbb{R_{\geq 0}} c(k(\omega)) \frac{\partial k(\omega)}{\partial\omega} \exp(i\,\omega t)\, d\omega \right)$$

where $k(\omega)$ is the inverse function of $\omega(k)$. Identifying

$$C(\omega) := c(k(\omega)) \frac{\partial k(\omega)}{\partial\omega}$$ for $\omega\geq 0$ and the complex conjugate for $\omega<0$, we can see that $f(0,t)$ is the Fourier transform of $C(\omega)$.

With the inverse Fourier transform we get

$$C(\omega) = \int_\mathbb{R}f(0,t) \exp(-i\,\omega t) dt$$

and consequentially

$$c(k(\omega)) = \frac{1}{k'(\omega)} \int_\mathbb{R}f(0,t) \exp(-i\,\omega t) dt$$

or simply

$$c(k) = \omega'(k) \int_\mathbb{R}f(0,t) \exp(-i\,\omega(k) t) dt$$

which we can express in terms of the inverse Fourier transform as

$$c(k) = \omega'(k)\,\mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right) $$

Using this theoretical background we can now reconstruct the spatial wave function of a free one-dimensional right travelling wave with open boundary conditions.

In the expansion of the general wave solution from above we can substitute the $c(k)$ to get:

$$f(x,t)=\int_\mathbb{R} \omega'(k)\,\mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right)\, \exp( i\, \omega(k) t )\exp( i\, k x ) dk $$

This can also be written as a Fourier transform in $k$:

$$f(x,t)= \mathcal{F_k}\left\{ \mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right)\,\omega'(k)\,\exp(i\,\omega(k)\,t)\, \right\}(x)$$

This translates to the following algorithm for the reconstruction of the spatial wave function at time $t$:

Take the inverse Fourier transform of $f(0,t)$
Evaluate the result at $\omega(k)$
Multiply the result with $\omega'(k)$
Evolve the result in time by multiplying with $\exp(i\,\omega(k)\,t)$
Take the forward Fourier transform with respect to k

You can implement this in the discrete domain by taking regular grids for the variables $\omega$ and $k$ and use the discrete Fourier transform.

As you can see, it's a little more complicated than what you described, but your idea was going in the right direction.


factoring poles / zeros: off by constant gain compared with textbook

Spectral factorization is not unique, unless you specify further constraints.

So let your given $H(z)H(z^-1)$ be (taking from your derivation) :

$$
\begin{align}
H(z)H(z^{-1})&= \frac{ \frac{5}{4} - \frac{1}{2}z - \frac{1}{2}z^{-1}  }{ \frac{10}{9}- \frac{1}{3}z - \frac{1}{3}z^-1 } \tag{1}\\
&= \frac{\left(1-\frac{1}{2}z^{-1}\right)(1-2z^{-1})}{\left(1-\frac{1}{3}z^{-1}\right)(1-3z^{-1})}
\end{align}
$$

Then all of the following are valid possibilities that satisf Eq.(1)

$$
H_1(z) =\frac{1-\frac{1}{2}z^{-1}}{1-\frac{1}{3}z^{-1}}
$$

$$
H_2(z) =\frac{1-2z^{-1}}{1-3z^{-1}}
$$

$$
H_3(z) =\frac{1-\frac{1}{2}z^{-1}}{1-3z^{-1}}
$$

$$
H_4(z) =\frac{1-2z^{-1}}{1-\frac{1}{3}z^{-1}}
$$

So you must have further constraints that make one of them as a unique answer. Such as minimum-phase, maximum phase, causal, stable, bothe etc...


Extended Kalman filter - why not use a 2nd order Taylor expansion?

There is nothing wrong with taking a higher order approximation than first order.
An issue in doing so is that you no longer have Gaussian statistics, so some of the standard Kalman filter equations (which are derived using Gaussianity) no longer hold strictly true.
The usual approach to dealing with (otherwise ignored) higher order terms in the standard EKF is to increase either the process noise or the measurement noise variances from what the model tells you they should be. This is like saying,

I know what the noise power is, but I'm going to say the second and higher order terms in my Taylor series expansion add to it so I'm going to increase the noise variance I use in my filter to cater for it.


What are the units of the product of two signals?

If the multiplier takes two voltages as input and returns a voltage as output, then there is necessarily a constant involved, with units of [1/V].

Take for example, AD633 (which was the first search result). The output is the product of the 2 inputs times a constant:

$V_{out} = \frac{V_1 \times V_2}{10V}$

So the output units are Volts.


First, use a timer and an ISR to get accurate timing (don't forget to configure the NVIC so that this timer interrupt takes over any other ISR that would be running). Only this will ensure a consistant sample rate. Little variations in timing would create noticeable degradations of audio quality. In particular, in your current example, unless the "filters" code has no branches, the time at which the SPI_Write occurs will be jittery depending on the path taken in the code.

For block processing, you need a "ping-pong" buffering scheme.

In your per-sample ISR:

SPI_Write(output[active_buffer][sample_index]);
SPI_Read(input[active_buffer][sample_index]);
++sample_index;
if (sample_index >= BUFFER_SIZE) {
  active_buffer = (active_buffer + 1) % NUM_BUFFERS;
  sample_index = 0;
}

In your main code:

while (render_buffer != active_buffer) {
  // Process samples from input[render_buffer][0] to input[render_buffer][BUFFER_SIZE - 1]
  // And write them to output[render_buffer][0] to output[render_buffer][BUFFER_SIZE - 1]
  render_buffer = (render_buffer + 1) % NUM_BUFFERS;
}

render_buffer and active_buffer are declared as volatile. Start with NUM_BUFFERS equal to 2. Higher values will introduce more latency, but might be necessary if the processing code has exceptionally long code paths.

The benefit of processing samples block by block is that some variables required by your DSP algorithm (such as filter state variables) can be kept in registers within the processing loop, saving load/stores. There are also situations in which you can benefit from batch load/stores. Finally, it makes it easier to schedule computations at "audio rate" (once per audio sample) and "control rate" (once per block).


Implementation of the filter function (Matlab) C/C++

FIR filtering is an operation whose implementation hugely benefits from utilization of architecture-specific instructions, compiler intrinsics, etc. Hence, if your only goal is to get a performant FIR algorithm, you should try to find a dsp library specifically for your target architecture. Those are most likely written in assembly language with a C interface and you can expect an orders-of-magnitude better performance than you would achieve with a generic C implementation of a transversal FIR.

If the latter is exactly what you're looking for - a generic FIR implementation in C - let me point you towards https://iowegian.com/. They offer a commercial filter design software called ScopeFIR. You can download and install a free trial version of it (no registration). In the program directory, you can find a source file "FirAlgs.c" that offers a variety of well-documented FIR implementations. It's written for MSVC but should be portable to basically anywhere without hassle.


Measuring frequency response with white noise signal

The question is about measurements in real world, not about simulations, am I right? If this is the case, there is an established technique and instrumentation for doing such measurements. The generic name of an instrument commonly used in laboratories and field testing for this purpose is Vector Network Analyzer. The article describes the operation principle in details and contains references for further research into the subject.
There is another kind of applications (audio environment acoustics, power grid control) where you need to measure frequency response of the system in real time (on the fly) -- or, admittedly, to do relatively instantaneous measurements of rapidly changing frequency response. In measurement of this kind, if the device under test (DUT) is linear, you can leverage the impulse response/frequency response duality, implementing pulse testing of your DUT. To be able to test a non-linear DUT, you can lower the crest factor (peak power/average power index) using a pseudorandom noise (as maximum length sequences that you mention) instead of single pulse stimulus.
But, to measure frequency response of "simple" networks "like an RC-filter", you need not go into such intricacies and can do with an ordinary lab equipment.


Recoverable continuous function after average sampling

A rectangular prefilter is a sinc(x) in the frequency domain, thus quite far from the text-book brick-wall filter presented in standard discrete sampling.

Still, if your (constructed) continous waveform contains no frequency components at >= samplingrate/2 and you can compensate for the high frequency roll off, I would think that near perfect reconstruction is possible.

By multiplying the spectrum of your continous waveform with that of the fourier transformed rect (sinc), you should be able to see how far from classical generic sampling you are.

Incidentally, I think that both audio A/D (time domain) and camera sensors (spatial domain) really does something close to a sample-period integrator/sample-and-hold (quite rect-like) prior to sampling. Presumably because this maximize the signal-related energy that can easily be harvested above the noise floor?

-k


You can use adaptive filters for this task. Feed the hardware with uniform white noise and collect the output. Then feed the adaptive filter with exactly the same noise sequence and use the output of the hardware you recorded as the desired response of the adaptive filter. The filter can then be trained with standard algorithm such as LMS. Of course, you will not get the exact transfer function, but the trained adaptive filter's transfer function should be close enough to serve as model of the hardware in question.


Music source separation: Incomplete notation in published IS-NMF/EM pseudo-algorithm?

The implementation in MATLAB is here:
https://www.irit.fr/~Cedric.Fevotte/extras/neco09/code.zip


To my eye, the equations in the slides look like $f(x) = a_1 x^1 + a_0 x^0 = a_1 x + a_0$. So it's an exponent on the $x$, not an index. As long as you are interpolating linear, the only exponents that you encouter are 1 (linear) and 0 (constant). When you go to quadratic, you'll see $x^2$ terms appearing. Now, $a_1$ is the coefficient of the linear term (your slope $m$) and $a_0$ is the coefficient of the constant term (the y-intercept $b$).

This should explain where the system of equations for $a_1$ and $a_0$ comes from.


Reminder what a permutation is: The exchange of indices in a sequence.

That's exactly what they're doing here.

Invertible in this context (only read the paragraph so far) means

"the multiplication with $\sigma$ under mod $n$ is invertible, i.e. given the result of that multiplication and $\sigma$, you can caclulate the original value in all cases."

Let's do an example:

let's say we're  $x\in \mathbb Z
\setminus 13$, and you multiply by $\sigma = 7$, you get $x\sigma = 4 \mod 13$. You know that $x=8$, undisputably! (because 7·8  = 56, and 56 mod 13 = 4, and it's only that easy because this is a prime field).

That works for any possible result $xz\in \mathbb Z\setminus 13$, because $\sigma$ is coprime to 13; so, $\sigma$ is invertible for $(\cdot,\mathbb Z \setminus 13)$.

I picked that 13 because it's prime, and prime numbers are by definition coprime to any smaller integer.

Let's do an counter-example of that:

For FFTs, the length is usually not a prime number, but can be any natural integer. $n=2^m, m \in \mathbb N$ is a very common FFT length (because you can decompose it easily using simple butterflies); then, not every $n$ is invertible in the sense above: for example, if your multiplication $x\cdot4 = 0 \mod 32$, you can't say whether $x=0$, $x=8$ or $x=16$; it's not invertible.


What does the convolution mean, what is the convolution philosophy

I'm somewhat in the same situation as you @Jony. But I'm working currently on a project involving audio processing, and I had to read a lot about convolution and related topics in the last few days.

I'm not sure I gained an intuitive understanding of convolution--esp. convolution vs. correlation-- but, to quote another answer, "Convolution in the time domain is equivalent to multiplication in the frequency domain." So it becomes obvious it may have a usage in filtering applications.

Wikipedia hs also a couple of nice animation that shows you how the convolution of two functions is related to the area of their intersection:

Depending on your field of interest, this may give you some applications idea. As the last word, let me repeat I'm very new to this topic, so take this answer with a grain of salt--and I would gratefully stand to be corrected by more experienced users


To what extent can we see signals that fall between frequency bins?

zero-padding will not "improve resolution" in the sense that it would help you separate two peaks in the spectrum that are very close (too close) to each other.

but zero-padding, as well as interpolation in the frequency domain, does improve resolution in the sense that it help you to more precisely estimate the location of an isolated peak in the spectrum.  in fact, zero-padding in the time domain is equivalent to interpolating in the frequency domain using the sinc() function (actually the dirichlet function) as the interpolation kernel.


Decreasing Sample Rate in DFT (FFT) for Audio Analysis

I've read that to increase frequency resolution of FFT results one should dicrease sampling rate and increase window size (number of samples).

To increase frequency resolution, you increase the window size (increase the number of samples).  The sample rate does not change.

For simplicity, I'll pretend your signal is sampled at 800 Hz.

If you were to take 4 samples of the audio, and run an FFT on it, you would get 4 samples in the frequency domain.  They would represent 4 frequency bins:

[0 Hz (DC), 200 Hz, ±400 Hz (Nyquist), -200 Hz]

Each of these bins is 200 Hz wide.  If instead, you had taken 8 samples of the audio, and run an FFT on it, you would get 8 samples in the frequency domain, representing 8 frequency bins:

[0, +100, +200, +300, ±400, -300, -200, -100]

Each of these bins is 100 Hz wide.  Does that make sense?  So to get narrower bins, you need to take more samples of your audio.  But the maximum frequency in the FFT (the Nyquist frequency of ±400 Hz) does not change; it's the same regardless of the number of samples.

Note that you can get these values like this:

scipy.fftpack.fftfreq(4, 1/800)
Out[8]: array([   0.,  200., -400., -200.])

scipy.fftpack.fftfreq(8, 1/800)
Out[9]: array([   0.,  100.,  200.,  300., -400., -300., -200., -100.])


Statement 1, as I read the garbled syntax, is correct, assuming "without a" and "doesn't have" imply a "zero" (an answer) and not a "null" (no answer).

The zero bin, also known as the DC bin, is the straight up sum of the signal over the sample frame in an unnormalized DFT calculation.  If the sum is zero, the average is zero.  So statement 2 is correct.

The explanation comes straight from the unnormalized definition of the DFT:

$$ X[k] = \sum_{n=0}^{N-1} { x[n] e^{ -i2 \pi { k \over N } n } } $$

Substitute in k = 0:

$$ X[0] = \sum_{n=0}^{N-1} { x[n] } $$

What is the source of these two statements?  Maybe there is a larger context to take into consideration.

Hope this helps,

Ced


Meaning of AWGN complex and real noise under BPSK Decision Scheme

It doesn't matter if the channel is complex or real. It's actually always complex, but since you are using BPSK, only the real part matters. How to add it is a mathematical trick. Say you need to generate an SNR of x in linear scale, and your signals have unit energy, then the noise variance must be 1/x. In other words, you don't generate real and complex parts for the noise, because then the SNR will be $1/(2x)$, which is not correct. You generate the real and complex parts only if the signal is complex, at which case the total noise variance gives the desired SNR. Hope this helps.


Simultaneously minimize and maximize two cost functions

Ideas:

1) Find the minimum eigen vector of $R_1$ and assign this to $w$. This will minimize $ J_1$ , doesn't maximize $J_2$. But when $R_1$ and $R_2$ are positive definite or full rank matrices and computation is an issue, this is a decent solution.

2) Form a new objective $J_1 - J_2$ and minimize this with the given constraint

3) Try and formulate as max SINR problem with $J_1$ as the interference term and $J_2$ as the signal term and the given constraint, you should be able to apply the relaxation $w^Tw < 1$ to convert this to a convex problem. This is readily solved in literature ( simply google (Max SINR) optimization.


You should use the the standard formula:

s = randn(m, n) + 1i*randn(m, n);

And as pointed out by MBaz, the output should be scaled accordingly by $\frac{1}{\sqrt{2}}$

s = s/sqrt(2);

More on that topic you can find for example in a book Digital Media Processing by Hazarathaiah Malepati:

It should generate exactly what you want.


Using NMF for guitar pitch detection using fixed pitch template

The goal of non-negative matrix factorization (NMF) is to factorize (as the name says) a matrix, denoted by $\mathbf{X}\in\mathbb{R}^{(N\times K)}$ in the reference you mention, into the product of two matrices:
$$
\mathbf{W}\mathbf{H}\approx\mathbf{X},
$$
with $\mathbf{W}\in\mathbb{R}^{(N\times M)}$ the template, dictionary, or basis matrix and $\mathbf{H}\in\mathbb{R}^{(M\times K)}$ the activation matrix. The approximation sign ($\approx$) in the equation indicates that in many cases it is not possible to actually find a true factorization but rather an approximation.

Inspecting the above equation reveals that each column of $\mathbf{X}$ is approximated by a linear combination of columns of $\mathbf{W}$, with the coefficients of linear combination in $\mathbf{H}$. (Experimenting with some simple examples, e.g., with $N=3,\,M=3, K=5$ and $\mathbf{X}_{ij}\in\left\{0,\,1\right\}$, will surely help gaining insight into this fundamental behaviour).

In many cases it is desired, given a matrix $\mathbf{X}$, to determine both the dictionary matrix $\mathbf{W}$ and the activation matrix $\mathbf{H}$. Optionally, of course, only one can be estimated if the other is given.

Your goal, if I understand correctly, is to keep $\mathbf{W}$ fixed, estimate $\mathbf{H}$ and use the resulting activations to determine when which note has been played. That should work as a first attempt, I think.

After this coarse overview now for your questions:

You are right: to create the $\mathbf{W}$ matrix put an estimate of each note's power spectral density estimate (e.g., obtained with Welch's method) into a column of $\mathbf{W}$ (I'm not sure what you mean by "do the FFT for each, take its spectrum".).
Yes, the matrix $\mathbf{X}$ in the notation used in this answer ($V^{(k)}$ in the paper that you mention) contains the STDFT magnitude spectrogram of the full input recording.
Silence in the beginning of the recording will obviously have an influence on the first columns of $\mathbf{X}$ (or, respectively, $v_j^{(k)}$ for $j<\text{some number}$). However, what do you guess will happen with $\mathbf{H}$ if the first columns of $\mathbf{X}$ contain zeros or only very small values? Correct: the first columns of $\mathbf{H}$ will also receive zeros or very small values.

I hope this answer helps in clarifying the basic functionality of the NMF algorithm. A rather intuitive description of the NMF idea is also given in this famous paper.


Matlab: Implement filter with no denominator

It's an odd notation but simple enough: $z^4% simeply means "shift the signal 4 samples to the left". That means its non-causal and can't therefore be implemented directly in a physical system.

In practice, the workaround that works most of the time is to add overall "causality delay" to the system until everything is causal again. Whether this is applicable here or not, depends on the specifics of your application.


According to the documentation and my past experience with the ETSI GSM 06.10 "GSM Full Rate Speech Transcoding" standard, you should be careful in the requirements of each block of the system and the length of the segments after the pre-emphasis or segmentation blocks.
Probably you made it, but check the length of the signal before running the LPC function. And remember that we can obtain p/2 formants of the voice. GSM uses p=8 in the standard.

Cheers.


Dimensional analysis in the Morlet wavelet

Okay, I think I've understood it now: The number $t$ in $\Psi_\sigma(t)$ is a dimensionless number, because in the formula for the wavelet transform it ends up being the ratio between an actual time and a time-scale.


Improve spectrogram appearance using windowing functions MATLAB

Windowing is nothing more than element-wise multiplication of your signal by window function.

Let's assume that you want to apply Hanning window. In your case signal is stored inside myRecording(start_sample:end).

Create the window values vector of your signal length: win =
hann(length(myRecording(start_sample:end)))
Multiply by your signal and store it in variable: sig = win.*size(myRecording(start_sample:end).
Perform the FFT, but keep in mind that now you cannot simply normalize but your signal length (as you do within FFT function), but you should (very simply speaking) divide by sum of your window values: sum(win). Probably you will notice that in case of no window applied = rectangular window, that is also true.

Answering second part of your question, which window to choose? You can start from this great white paper (on figure 8 you have recomendations for different window types): Windowing: Optimizing FFTs Using Window Functions.

Just keep in mind that if you you increase side-lobe attenuation then you also make your main-lobe wider (worse frequency resolution). Hanning is a good start.

For more info you can refer to literature:

  Harris F.J. - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform
  Heinzel G. - Spectrum and spectral density estimation by the DFT

Also I've noticed in your code that you are normalizing by the number of samples and then taking only first half of your spectrum. That's perfectly ok, but you must also multiply amplitude by 2 in order to compensate for removing of half of your energy in components above $\frac{f_s}{2}$

In case you have more queries you can search for questions and answers regarding windowing and leakage - you will find plenty.


On the de-noising of time-frequency images

I don't have experience in this area but I see that it has been studied: Minimum entropy approach to denoising time-frequency distributions

  In this paper, we introduce an entropy based approach to denoising time-frequency distributions. This new approach uses the spectrogram decomposition of time-frequency kernels proposed by Cunningham and Williams. In order to denoise the time-frequency distribution, we combine those spectrograms with smallest entropy values, thus ensuring that each spectrogram is well concentrated on the time-frequency plane and contains as little noise as possible. Renyi entropy is used as the measure to quantify the complexity of each spectrogram. The threshold for the number of spectrograms to combine is chosen adaptively based on the tradeoff between entropy and variance.

Essentially your problem is one of signal/source separation; the additive unmixing of a bunch of structured signals. In order to proceed you need to model your signals. Obviously the one of interest is periodic and is centered about some frequency, so you need to estimate the period (along the x axis) and the center frequency (on the y axis). Then you can characterize the others (noise). For starters, it seems that they come in nice curves.

With a model in hand I would consult a book like Handbook of Blind Source Separation: Independent Component Analysis and Applications.


Generate signals with a particular variance and SNR

It depends on what you mean by SNR. It's a common joke in the DSP community to spell it out as "something to noise ratio", referring to the fact that there is no unique definition of SNR, so the term by itself means nothing.

Define it yourself and use it appropriately.

What's common is to define it as ${\rm SNR} = \frac{P_{\rm s}}{P_{\rm n}}$ where $P_{\rm s}$ is the power (variance) of the signal samples ($x_n$ in your notation) and $P_{\rm n}$ is the power (variance) of the noise samples.

Hence:

P_s = 1;   % target signal power
SNR = 15;  % target SNR in dB
P_n = P_s / 10^(SNR/10); % calculated noise power

s = randn(1,N)*sqrt(P_s);
v = randn(1,N)*sqrt(P_n);
y = a*s + v;

You don't need the for loop there by the way. You can use the awgn function, it should do something very similar. Again, there is no clear need for it though.

Note that this definition of the SNR is independent of a. Depending on your application this may be desirable or undesirable. In the latter case, you can include it in your SNR definition. Feel free to do so. Just define everything very clearly.

*edit: In response to your question: the SNR is really a ratio of signal power to noise power. So to change the SNR, you can change both. If you want the SNR to be 10 dB higher, you can decrease the noise variance by 10 or alternatively increase the signal variance by 10.

If you want to change the SNR while keeping the noise variance fixed (say, 1), you're welcome to do so. Almost the same approach:

P_n = 1;   % target noise power
SNR = 15;  % target SNR in dB
P_s = P_n / 10^(-SNR/10); % calculated signal power

s = randn(1,N)*sqrt(P_s);
v = randn(1,N)*sqrt(P_n);
y = a*s + v;

This also means looking at the noise variance itself doesn't tell you about the SNR conditions. If you assume it's 1 this can mean a low SNR (if the signal power is, say, below 1) or a high SNR (if it is much higher than 1).


MMSE - How to minimize a complex error with respect to a set of real parameters

A very common approach is to consider $X(k) \;\text{and} \;\hat{X}(k)$ as  elements vector space $C^R$, and consider the distance between the 2 vectors as a norm, which are real, positive or zero, and satisfy the triangle inequality.

So using vector-matrix notation with the vectors as column vectors:

$$ \text{error}^2 =(\mathbf{X}-\mathbf{\hat{X}})^H
(\mathbf{X}-\mathbf{\hat{X}})
$$
where $H$ is conjugate transpose.

Constraint is satisfied when imaginary part of $A_r$ is zero.

$$
\frac{1}{2}(A_r - A_r^*)=0
$$
which is appended to the objective as a lagrange multiplier(s)
$$
\text{error}^2+\sum \lambda_r \frac{1}{2}(A_r - A_r^*)
$$

These sorts of problems are easier when you use Brandwood derivatives,

D. H. Brandwood, "A complex gradient operator and its application in adaptive array theory," in Communications, Radar and Signal Processing, IEE Proceedings F, vol. 130, no. 1, pp. 11-16, February 1983.
doi: 10.1049/ip-f-1.1983.0003
Abstract: The problem of minimising a real scalar quantity (for example array output power, or mean square error) as a function of a complex vector (the set of weights) frequently arises in adaptive array theory. A complex gradient operator is defined in the paper for this purpose and its use justified. Three examples of its application to array theory problems are given.
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4645581&isnumber=4645575

This is the convention used in

Van Trees, Harry L. Optimum array processing: Part IV of detection, estimation and modulation theory. Vol. 1. New York, NY, USA: John Wiley & Sons, 2002.


How to Apply Statistical Algorithms of Signal Processing to Regulate Variation of a Curve?

You can regulate its second derivative which is the curvature.

Something like:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} \left\| x - y \right\|_{2}^{2} + \frac{\lambda}{2} \left\| D D x \right\|_{2}^{2}  $$

Where $ y $ is samples you have and $ D $ is the Derivative Operator (In Matrix Form). By applying it twice we're regulating the Second Derivative.
You can find a statistical explanation for this (Prior of Gaussian Distribution of the Second Derivative) or just call is Tikhonov Regularization.

This has a closed form solution:

$$ \hat{x} = \left( I + \lambda {D}^{T} {D}^{T} D D \right)^{-1} y $$


Designing an efficient curve-matching algorithm

The two most obvious things you can try are:

Fitting a Gaussian to your data and then clustering their parameters
Estimate the similarity of waveforms directly and then try to cluster that

Since you know that the return waveform conforms to a Gaussian, it is better to use a method that takes this into account.

So, basically, for every pixel time course, you will end up with a fitted $\mu, \sigma^2$. You can then use these two, plus any other target dependent feature to characterise the surface.

  These algorithms work somewhat fine but since the signal can be noisy and sometimes very far from a Gaussian shape I think and I hope better algorithms can be implemented.

The fact that you know that your data is supposed to conform to a specific model is a great deal of knowledge and you should try to take that into account as much as possible. "...sometimes very far from Gaussian..." doesn't always mean "not Gaussian". You might for example be getting two or more returns in the same time-course which means that you might have to apply something like a kernel density estimator or an iterative fitting process to pick up each individual return. For more information, please see this paper.

The second option is to estimate "similarity" of waveforms directly. In this case, given two time series $U,V \in \mathbb{R}^\mathbb{N}$, you are looking at a function $F(U,V,\Theta) \rightarrow \mathbb{R}$ which simply returns a number that is proportional to how "similar" these two waveforms are, given some additional parameters in $\Theta$. Most of the times, "similar" is another word for "co-varying".

The most obvious of the possible $F$s is the discrete cross-correlation. In that case $f_{xcorr}(U,V)$ returns a number between $-1..1$ which tells you how similar $U,V$ are. Another possible $f$ is Coherence. In this case, parameters include the frequency component of which you want the coherence calculated at (e.g. $f_{coher}(U,V,k)$ (for the $k^{th}$ component) which returns a number between $0..1$. Another possible $f$ is Mutual Information with is a slightly different take on "similarity" and more challenging to apply correctly on continuous (rather than digital) signals.

And then, on top of the metrics mentioned above you have variations, like wavelet based cross-correlation instead of correlation, LZ complexity instead of Mutual Information and so on.

The field of Similarity Metrics is vast and there is a lot of work happening in the processing of Electroencephalography (EEG) signals in Functional Connectivity (which basically involves a lot of those evaluations of similarity under very challenging conditions of waveform noise and artifacts). For more information, please see this, this and this.

These publications are not exactly on your topic but they use a great deal of similar methods and tools for estimating the similarity between two time-series which is what you are trying to do as well. Therefore, when reading those, try to extrapolate the terminology and then try to search in remote sensing related journals for similar things. I am sure that such a search would return much more relevant results adapted to the specifics of estimating similarity of the laser pulse returns.

Hope this helps.


Recover myself from Dilate and Erode for noise reduction

First of all, are you using gray-scale morphology?  If so, then I would threshold the map first, to get a binary image, and then apply binary morphology to fill in gaps and get rid of the noise.

Second, instead of erosion and dilation, try opening and closing, which are higher order morphological operations. Opening is erosion followed by dilation with the same structuring element.  Closing is dilation followed by erosion.

A good thing to try is to do an opening to get rid of the noisy specks, and then do a closing to fill in the gaps.  This will simplify your process to three operations: threshold, open, close.  So you would only need to tune three parameters: the threshold and the two structuring elements.


How to remove sum and frequency components in FFT spectrum during signal resampling.?

You should be using an appropriate low pass filter during reconstruction.


good free software (or code) to separate multiples sources in a single wav?

A great free toolbox is the Flexible Audio Source Separation. Their website contains some demo examples in Python, Matlab, and C++ code.

Another toolbox is the Toolbox for Blind Source Separation. However, this toolbox is not user friendly and you need to do some extra programming in C, Matlab or Python.


I had a similar question myself, so I created a html demo to test this.

I think you may find it useful in understanding what the FFT is actually doing since you can also read and edit the source code.

https://seanwasere.gitlab.io/FFTBinary

This demo is doing FFT analysis of binary data. 1s and 0s

In my demo I have 1024 bits, and I convert them into an 8 bit wav, and then analyse it using the browsers inbuilt audio api fft analysis function. Firefox browser works best in my case.

The 8 bit wav is

1 channel
samplerate 48000
1024 bits = 0.021333333333333332 seconds in length
the 1s and 0s are actually converted behind the scenes to 256s and 128s

In the user interface there are 4 panels

1st is the FFT result of the binary data drawn as a spectrum graph,
2nd is the FFT result drawn as a waveform. It is essentially triangular spikes or squares, depending on the sequence 0,1,0 or 0,1,1,0.
3rd panel is a view of the actual binary data that was converted into the waveform before FFT analysis.
4th panel is your options, you can

set to have a '1' every N bits by using modulus,
and/or play random 1s and 0s which I call binary white noise.
set the FFT buffer size, the lower numbers such as 128 give lower quality output but is much faster to compute.

The source code is MIT so you can play around with it.


BIBO Stability and the convergence of the frequency response of a system

The same happens with sinusoids: they have a Fourier transform even though they are not absolutely summable in the time domain. These are cases in which the Fourier transform doesn't converge uniformly, but it does converge in some other sense.

The signal $x(n)=\mathrm{sinc}(n)$ is not absolutely summable... but it is square-summable (i.e. it has finite energy). Therefore, the Fourier transform will converge but in the mean-square sense. In these cases, the transform converges to a discontinuous function (in the case of this specific signal, to a rectangular pulse).

Similarly, sequences that involve sinusoids like $x(n)=\sin(\omega_0n)$ are neither absolutely nor square summable, but a Fourier transform can be defined using distributions (Dirac delta functions).

In both cases the Fourier transforms are not continuous. They are not infinitely differentiable either.

As a general rule we can state that if some impulse response $h(n)\in l^1(\mathbb{Z})$, then its Fourier transform will be continuous. If the Fourier transform is not continuous, then you know that its inverse is not in the sequence space $l^1$.

You mentioned BIBO stability, so I think this may be of interest. There is a theorem that states:

  An LTI system with impulse response $h(n)$ is BIBO stable if and only if $h(n)\in l^1(\mathbb{Z})$.

So one can judge the stability of a system by looking at the continuity of its transfer function.


How to handle a logarithmic term in Kalman filter?

First of all, you're trying to evaluate the derivative of an exponential. If the base of the exponential is positive, the derivative exists. However, if the base is negative, the derivative does not exist. It would be pointless to use a Kalman filter when you have an exponential with a negative base.
Secondly, your model of the echo envelope is
$$
A(t) = A_0(\frac{t-\tau}{T})^{\alpha}e(\frac{t-\tau}{T})
$$
However, since it is an echo, it does not make sense for the echo envelope to be something else than 0 for $t < \tau$ where $\tau$ is the time of flight
You could add a unit step to the equation to make sure that A(t) = 0 for t < $\tau$
$$
A(t) = A_0(\frac{t-\tau}{T})^{\alpha} e(\frac{t-\tau}{T})  u(t-\tau)
$$


PSOLA how to calculate synthesized Marks

Let $t_a(n)$ be the time of the n-th analysis mark ("Pitch marker"), $P_a(n)$ the period at the n-th analysis mark (it can be defined as $t_a(n+1) - t_a(n)$, $t_s(n)$ be the time of the n-th synthesis mark, $n_s(k)$ is a number whose integer part is the index of analysis segment to be overlapped-added when processing the n-th synthesis mark (it can also represent a weighted interpolation between two consecutive analysis marks).

The synthesis marks will be defined by the recurrence:
$$n_s(0) = 0$$
$$t_s(0) = t_a(0)$$
$$t_s(k) = t_s(k-1) + \alpha(n) P_a([n_s(k)])$$
$$n_s(k) = n_s(k-1) + \alpha(n)$$

$\alpha(n)$ is a stretch factor. It is 0.5 if you want to double the frequency. If you want unvoiced frames to be processed naturally (no stretching for them), you can actually make it dependent on the voiced/unvoiced analysis. In this case: it is 1 if the analysis mark $[n_s(k)]$ is unvoiced, and 0.5 or whatever is your pitch-shift factor when it is voiced.


What's the difference between calculating the FFT of a constant sound for 5 seconds versus every 1 second and taking the average?

I would say that you answered part of your own question by noting the reduction of variance with averaging.

Most books tend to present signals that are either random or deterministic but  signals can be somewhat of both character to varying extent.

If your signal tends towards the deterministic, and high SNR, there isn’t a lot variation to reduce so the higher resolution of the longer FFT would probably be more preferable in most applications.

If you have a low SNR tone, you want to do as much averaging as other considerations would permit to reduce that noise.

One can also overresolve with long FFTs.  Pure steady tones are rare. Real tones have some bandwidth and you want your FFT bins to match that bandwidth.

Averaging has some other advantages like reducing the down stream data rates like if you are displaying real time STFT displays.

Depending on your application, there are other more specific tradeoffs between averaging and FFT length.


Maximum Magnitude Deviation between DFT and DTFT

okay, i'm gonna steal a diagram from this publication:

now this is an approximation that is based on the assumption that the discrete spectral points $|X[k]|$ are close enough to each other that only the discrete peak:

$$m_\textrm{d} = \max_k |X[k]| = |X[m]|$$

and the two adjacent neighbors $|X[m-1]$ and $|X[m+1]|$ need be considered. (these are shown as $k$ in the diagram above, not $m$, but i don't know if i should change the convention of discrete variables we're using and i don't wanna draw out another diagram.)

the assumption is, within this very tight range of $\omega$, that the peak is adequately approximated by a quadratic:

$$ |X(\omega)| = m_\textrm{c} - a(\omega - \omega_\textrm{m})^2 $$

where $a > 0$, and the "spectral peak":

$$\star \quad |X(\omega_\textrm{m})| = m_\textrm{c} = \max_\omega |X(\omega)|$$

so we have three unknowns: $|X(\omega_\textrm{m})| = m_\textrm{c}$, $a$, and $\omega_\textrm{m}$ and three equations by evaluating this quadratic at the discrete points $|X[m-1]|$, $|X[m]|$, $|X[m+1]|$:

$$\begin{align}
|X[m-1]| &= |X\left(2\pi \tfrac{m-1}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m-1}{N} - \omega_\textrm{m})^2 \\
|X[m]| &= |X\left(2\pi \tfrac{m}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m}{N} - \omega_\textrm{m})^2 \\
|X[m+1]| &= |X\left(2\pi \tfrac{m+1}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m+1}{N} - \omega_\textrm{m})^2 \\
\end{align}$$

so you can solve those three equations for the three unknowns.  we're not so interested in $a$, but we are interested in the frequency location of the peak $\omega_\textrm{m}$ and the height of the peak, $m_\textrm{c}$.  if you slug out solving this set of three equations and three unknowns, you get for those two parameters:

$$\begin{align}
\omega_\textrm{m} &= \frac{2 \pi}{N}\left( m + \tfrac12 \frac{|X[m+1]| - |X[m-1]|}{2|X[m]| - |X[m+1]| - |X[m-1]|} \right) \\
\\
|X(\omega_\textrm{m})|= m_\textrm{c} &= |X[m]| + \frac{\tfrac18 \big( |X[m+1]| - |X[m-1]| \big)^2}{2|X[m]| - |X[m+1]| - |X[m-1]|}   \\
\end{align}$$

Now this makes no assumptions about the nature of the spectrum except that the discrete spectral peak is at $|X[m]|$ and the neighboring spectral samples, $|X[m-1]|$ and $|X[m+1]|$ are no larger.  With no other information, this is about the best or only guess you can make.  It's not always perfect, but it should get you to within one-half FFT bin of the true DTFT peak.

It is perfect if you're looking at the peak of an isolated sinusoid, that was windowed with a gaussian window before going to the FFT, and you're looking at the dB or log magnitude.  Then this quadratic behavior in the log magnitude scale is mathematically the correct model and this quadratic peak interpolation is exactly correct.  So then you replace $|X[m]|$ and $|X[m-1]|$ and $|X[m+1]|$ with the logarithms (or dB values) of the same before using the above formulae to get the continuous-frequency peak.


Difference between convolving before/after discretizing LTI systems

Assuming you refer to the ZOH-discretization as shown in this figure from the mathworks site

the relationship between the continuous-time signals and the discrete-time signals can be derived as follows. The signal $u(t)$ is given by

$$u(t)=\sum_ku[k]g(t-kT)\tag{1}$$

where $g(t)$ is a rectangular impulse response (constant in the interval $t\in [0,T]$, zero anywhere else) representing the ZOH, and $T$ is the sampling period. Let $h(t)$ be the impulse response corresponding to the transfer function $H(s)$. Furthermore, let $f(t)$ be the convolution of $g(t)$ and $h(t)$:

$$f(t)=(g\star h)(t)\tag{2}$$

Then the continuous-time output signal $y(t)$ is given by

$$\begin{align}y(t)&=(u\star h)(t)\\
&=\int_{-\infty}^{\infty}u(\tau)h(t-\tau)d\tau\\&=\sum_ku[k]\int_{-\infty}^{\infty}g(\tau-kT)h(t-\tau)d\tau\\&=\sum_ku[k]f(t-kT)\tag{3}\end{align}$$

where I've used Equations $(1)$ and $(2)$. From $(3)$ the discrete-time output signal is easily obtained as

$$y[n]=y(nT)=\sum_ku[k]f((n-k)T)=\sum_ku[k]f[n-k]\tag{4}$$

with $f[n]=f(nT)$. So $y[n]$ is simply the discrete-time convolution of the input signal $u[n]$ and the total impulse response $f[n]$.

From this result it immediately follows that linearity must be satisfied, simply because convolution is a linear operation:

$$f(t)=\alpha_1(g\star h_1)(t)+\alpha_2(g\star h_2)(t)=(g\star (\alpha_1h_1+\alpha_2h_2))(t)\tag{5}$$

However, the (ZOH-)discretization of the concatenation of two transfer functions is not equivalent to the concatenation of the two discretized transfer functions. Referring to the figure above, the first case is equivalent to replacing $H(s)$ by the concatenation of $H_1(s)$ and $H_2(s)$. The other case involves concatenating two complete systems as shown in the figure, the first with $H(s)=H_1(s)$, the second with $H(s)=H_2(s)$. The difference between the two cases is that in the first case you only have one ZOH, whereas in the second one you get two ZOHs. The equivalent discrete-time impulse responses are

$$f[n]=(g\star h_1\star h_2)(nT)\tag{6}$$

and

$$f[n]=(g\star h_1)\star (g\star h_2)(nT)\tag{7}$$

which are generally not identical. The first case in Eq. $(6)$ is a ZOH discretization of the total transfer function $H(s)=H_1(s)H_2(s)$, whereas the second case in Eq. $(7)$ corresponds to a first-order hold (FOH) discretization of the total transfer function.


Does the output signal from the FIR filter have transient response?

Probably. Except for 1-tap filters, if your signal is of finite support, FIR filtering needs to assumes  values outside the support (before and after the values you were given)), and those values are often considered to be $0$. Hence, the extended sequence $(\ldots,0,0,x[0],x[1],x[2]\ldots)$ may produce some artifacts or transients in the output around time $[0]$, depending on the continuity and smoothness of $x$ close to $k=0$.


Fourier Transform Problem - absolute value, time-saving tricks, etc

As you suggested, you could simply solve the integral using $\cos(\omega_0t)=(e^{j\omega_0t}+e^{-j\omega_0t})/2$. But as you've also noted, with this type of problems there is usually some smart way using known transform pairs. And what you suggested appears to me a very sane approach: convolve the known transform of $\cos(\omega_0t)$ with the known transform of $e^{-at}u(t)$ and you're done. I would recommend to you to cross-check your result by solving the integral, which in this case is also quite straight-forward. But I can assure you that your result looks good.


why eigenvalues concerned in Harris corner detection?

Basic linear algebra. M is a symmetric matrix, so there's a decomposition:

$M=Q^{\mathsf{T}} \Lambda  Q$

where Q is a rotation and $\Lambda$ is a diagonal matrix.

We can rotate $[\text{$\Delta $x},\text{$\Delta $y}]$ by $Q^{\mathsf{T}}$

$Q^{\mathsf{T}} [\text{$\Delta $x},\text{$\Delta $y}]^{\mathsf{T}}=\left[\text{$\Delta $x}',\text{$\Delta $y}'\right]^{\mathsf{T}}$

giving

$S(x,y)=\left[\text{$\Delta $x}',\text{$\Delta $y}'\right]\Lambda \left[\text{$\Delta $x}',\text{$\Delta $y}'\right]^{\mathsf{T}}=\Lambda _{1,1} \left(\text{$\Delta $x}'\right)^2+\Lambda _{2,2} \left(\text{$\Delta $y}'\right)^2$

So $S(x,y)$ is large in both directions $\text{$\Delta $x}'$ and $\text{$\Delta $y}'$, if both eigenvalues are large.


How to reconstruct signal of its phase and magnitude functions?

You can look at the domains covered by the Fourier transform (FT) or the DTFT (Discrete Time Fourier Transform). The DTFT applies to discrete signals, and yields a (possibly $2\pi$) periodic transformation. Here, $\omega$ denotes the continuous normalized radian frequency variable. So:

$$ (s[k])_{k\in \mathbb{Z}} \stackrel{\mathrm{DTFT}}{\longrightarrow} S(\omega) = \sum_{k=-\infty}^{\infty} s[k] z^{-k}, \quad z=e^{\imath \omega}$$
while for the FT:

$$s(t) \stackrel{\mathrm{FT}}{\longrightarrow} S(\omega) = \int_{-\infty}^{\infty} s(t) e^{- \imath \omega  t}dt$$
with no specific periodicity in general.

So most likely, you are asked to recover a discrete signal. The complex quantity $S(\omega)$ can be recovered with  the amplitude and the phase you have. You can then recover the discrete signal using the inverse DTFT over one period:

$$s[k] =  \int_{-\pi}^{\pi} S(w)e^{\imath \omega k} d\omega \,.$$

You definitely need to check the  Discrete Time Fourier Transform (DTFT) and the associated pages.


How to find the output mean and autocorrelation of a non-linear system

Put your second equations into your first equation, express $Y(t)$ as a function of $S(t)$ and $N(t)$
Apply the definition for mean and autocorrelation.
Simplify and solve


Beginner question on aperiodic samples and DFT behavior

It sounds to me like you are experiencing time-domain aliasing. First, two basic facts:

Recall that if $x[n]$ is of length $N$ and $h[n]$ is of length $M$ then $x[n]\ast h[n]$ is of length $N+M-1$.

Fact 1 does not matter in what domain you perform the convolution/filtering. So, if you modify the DFT of $x[n]$ by "shaping" its spectrum, then there exists an  equivalent filter $h[n]$ in the time domain, which defines an implicit length $M$.

Taking an $N_1$-point IDFT of an $N_2$-point signal for $N_1 < N_2$ results in time-aliasing.

Fact 2 is the dual to frequency aliasing by undersampling. I'll illustrate this as follows: let $x_1[n]$ denote a signal of size $N_1$. Then, if we take an $N_1$-point DFT of $x_1$ and modify its spectrum so the size should be $N_2>N_1$, then taking an $N_1$-point IDFT is equivalent to obtaining the signal $x_2$ given by:
\begin{eqnarray}
x_2[n] = \sum_{p=0}^{\infty}x[n+pN_1]
\end{eqnarray}

This seems to be consistent with what you have mentioned. As you correctly noted, zero padding the original signal is exactly the solution! Essentially just zero-pad such that the terms in the above equation for $p>0$ are all zero-valued. Also, using the ideas above gives you a way to understand the exact amount of zero-padding necessary based on the implicit filter $h[n]$ your frequency domain shaping yields. Good luck!


Mathematical description of the ACF using fft2

When the signal is assumed to be ergodic, then its ACF can be computed using time averages which can also be computed using the following convolution of length $N$ sequence $x[n]$ by its conjugate symmetric version:

$$ \hat{r}_{xx}[m] = \frac{1}{N} ~~~x[m] ~~\star ~~ x^*[-m] ~~ $$

Taking the $2N-1$ point DFT of both sides yields the following:

$$ \hat{R}_{xx}[k] = \frac{1}{N} ~~~X[k] X^*[k] =  \frac{1}{N} |X[k]|^2  $$

In a MATLAB implementation, an fftshift is required to bring $r_{xx}[0]$ to the center position for convenience. Hence the ACF will be obtained by an $2N-1$ point inverse DFT/FFT as

rxk = (1/N)*real(fftshift( ifft( abs(fft(x,2N-1)).^2 ,2N-1) ) );

The 2D case follows very similarly to this 1D outline...


Why the bandwidth of the system with the minimum-order observer is higher than the full-order observer?

I think the key of that statement is the sentence

  [...] provided that the multiple observer poles are placed at the same place for both observers.

Suppose that you have a Bode plot, with poles at $\omega_1,\omega_2,...,\omega_n$. If you leave those poles as they are (i.e. you don't change their locations) and you add more poles to the system, then the bandwidth will decrease.

Why? Well, suppose we add a single pole at $\omega_{n+1}$. In that case, the slope of the magnitude as a function of the frequency will decrease for all the frequencies such that $\omega>\omega_{n+1}$. If the magnitude was increasing, it will increase slower or will become constant  (referring to the asymptotic Bode plot). If it was decreasing, it will decrease faster.

If you add even more poles and not just one, then, following the same reasoning, the bandwidth will decrease more abruptly.


Matlab is right, as usually. The DFT of a sine is imaginary. If you do the same with the sum of two cosines, you would get a real-valued result.


How to write e^(-1000|t|) in matlab title?

You can write e^(-1000*|t|) in Matlab in colored form if you want as

title('\ite^{-1000 |t|}','Color','b')
Here,"b" is for blue color ,"r" for red,"g" for greeen.

For subscript use _

For superscript use ^

For exponential symbol use \ite ,

For lambda symbol use \lambda

For alpha use \alpha.


We give a precise description of the pyramid algorithm in this paper: http://www.ipol.im/pub/art/2014/79/

It's a wavelet transform.

If you like to use the exact filters you can use Fourier transform but the filters are usually well approximated in spatial domain on a small support (8 x 8 px^2). Note that the direct convolution is faster than the Fourier implementation when the filters have few pixels (less than 16 x 16 px).


Calculating Affine Transformation matrix from 4 corner points

If you already have the pairs of corresponding two points from each coordinate,
you can use the function fitgeotrans
 to calculate the transformation matrix.

Once you have got the transformation matrix, the transformation can be done with imwarp.

Note that you are going to use Homography transformation rather than affine transformation in this case (because you have specified corresponding 4 points, which may cause trapezoidal distortion).

The code would be like this...

% assume that the input image is in the variable inputImage.
nPixelsX = 50; % assume that there is square image with 50 pixels for each axes after transformation.
nPixelsY = 50;

movingPoints = [73,86; ...
                117,88; ...
                70,127; ...
                115,128]; % set image's corresponding points
fixedPoint = [0,0; ...
              nPixelsY,0; ...
              0,nPixelsX; ...
              nPixelsY,nPixelsX]; % set ideal coordinate's corresponding points
transformationType = 'Projective';  % this should be homology transformation
tform = fitgeotrans(movingPoints,fixedPoints,transformationType); % make transformation matrix

resultReference = imref2d([nPixelsY, nPixelsX]);
desiredImage = imwarp(inputImage, tform, 'OutputView', resultReference); % transform


How to reconstruct the signal which is subjected to FIR filtering - Matlab

You can if the filter does not introduce any zeros in the frequency domain, which your filter unfortunately does because its frequency response is zero for $\omega=0$. What you normally have to do is zero-pad the impulse response $h(n)$ to the same length as the result of the convolution, i.e. $y(n)$. Then you compute the DFT (FFT) of the two zero-padded sequences:

$$H(k)=DFT\{h(n)\}$$

$$Y(k)=DFT\{y(n)\}$$

Now you have

$$Y(k)=H(k)X(k)\tag{1}$$

where $X(k)$ is the DFT (FFT) of the zero-padded input signal (such that it has the same length as $y(n)$). From (1) you can compute $X(k)$ for all indices $k$ for which $H(k)\neq 0$. In your case $H(0)=0$, so you can't compute $X(0)$. You can, however, compute all other values of $X(k)$ from (1).

For your example this means that you can determine the original time signal $x(n)$ from $X(k)$ up to its mean value. E.g., if you choose $X(0)=0$ and apply an inverse DFT you'll get the original signal but with a mean value of zero.

In Matlab it would look like this:

Nx = 100; x = randn(Nx,1);    % some input signal
h = [1,-2,1]'; Nh = 3;        % impulse response
y = conv(x,h);                % convolution
h0 = [h;zeros(Nx-1,1)];       % zero-padding
Y = fft(y); H = fft(h0);
X = Y(2:end)./H(2:end);       % undo filtering except for DC
X = [0;X];                    % add arbitrary DC value
x_ = real(ifft(X));           % reconstruct input signal (up to mean)
plot(x-x_(1:Nx));             % should be constant


What are the state of the art methods for person segmentation and person pose identification

Check out Yang and Ramanan's work with flexible mixture of parts.  It starts with a central part, but pairs of parts that are attached are detected jointly by estimating the probability that the image area being queried supports the co-occurence of the attached parts in the image.  Error is presumed in the attachment and is minimized using a spring-like model.  Each part is detected using HOG-style techniques, and all pairs of parts in the human 'tree' are iterated over to determine the detection.  In my experience it has a very good false alarm rate, but requires hardware acceleration for real-time applications.


Segmentation of half transparent material, e.g. glass

Why not just use a simple 2D FFT (guassian) high pass filer?

I did this real quick using MATLAB

Shard #1 using high pass FFT:

The same thing is done on #2.

Shard #2 using high pass FFT:

As you can see, the background and glass area is wiped out, and only the edges are traced.  I did not spend any time on it, but you can threshold the HP filtered output to have more crisp edges, or push the HP cut off higher.

Is this more the results you are looking to get?


Kalman filter for random observation matrix $G_t$

Very interesting question!

I do not know of any reason why you can't do this, but I have had a hard time finding an example.  I finally found this paper, which does an analysis of what happens when the measurements are randomly lost due to transport delays / errors in communication systems.

The main problem that they report is that the error covariance matrix iteration is now stochastic rather than deterministic (as in the usual Kalman filter formulation).


How to determine fundamental frequencies (beats/minute) of heartbeat? (matlab)

I would suggest processing the ECG signal through a standard QRS detection algorithm (e.g. Pan-Tompkins). Once you have isolated the QRS complexes in your ECG signal, you can simply count them, divide by the duration of your ECG signal, and convert to beats/min. Here is one implementation of Pan-Tompkins: https://www.mathworks.com/matlabcentral/fileexchange/45840-complete-pan-tompkins-implementation-ecg-qrs-detector. I am sure you can find others on physionet.


Filter coefficients for colored noise Voss Algorithm

Here, $a$ is a filter coefficient for a simple single-pole recursive low-pass filter. If you look at the documentation for MATLAB's filter() function, you'll see this corresponds to the difference equation $y(n) = x(n) - a y(n-1)$`. Let's analyze the filter in the $z$-domain to get its frequency response (just remembering that $z^{-1}$ corresponds to a single sample delay:

$$y(n) + a y(n-1) = x(n)$$
$$Y(z)*(1 + a z^{-1}) = X(z)$$
$$\frac{Y(z)}{X(z)} = \frac{1}{1 + a z^{-1}} = H(z)$$

That's the $z$-domain transfer function. We can substitute $z = e^{j\omega}$ where $\omega$ is the normalized frequency, where 0 is DC and $\pi$ corresponds to the Nyquist frequency of whatever sampling rate you're at. Now you should realize that if you use the same digital coefficients at different sampling rates, you will get filtered noise that varies depending on the sampling rate. If your sampling rate is $F_s$, then $\omega$ corresponds to the frequency $\dfrac{\omega F_s}{2 \pi}$. Let's find the frequency response:

$$G^2(\omega) = |H(e^{j\omega})|^2 = \left| \frac{1}{1 + a e^{-j\omega}} \right|^2$$

$$G^2(\omega) = \frac{1}{|1 + a(\cos(-\omega) + j \sin(-\omega))|^2}$$

$$G^2(\omega) = \frac{1}{|(1 + a \cos(\omega)) - j \sin(\omega))|^2}$$

$$G^2(\omega) = \frac{1}{(1 + a \cos(\omega))^2 + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + (1 - \cos^2(\omega))}$$

$$G^2(\omega) = \frac{1}{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}$$

$$G(\omega) = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

Let's plot the frequency response with $a = -0.9$: WolframAlpha plot. We can also take 20 Log10() to get the frequency response in decibels: WolframAlpha plot .

Now we can find out some things about this filter. For example, its gain at DC is greater than unity! $G(0) = \dfrac{1}{\sqrt{1 + 2a + a^2}}$. When $a = -0.9$, the filter increases the DC gain by a factor of 10, or 20 dB. Often filters are characterized in terms of their -3 dB point, which is where their magnitude decreases by $\sqrt{2}$ from their peak. Let's solve for this point:

$$\frac{1} { \sqrt{2} \sqrt{1 + 2a + a^2} } = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

$$\sqrt{2} \sqrt{1 + 2a + a^2} = \sqrt(2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega))$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

Now treat this as a quadratic in $\cos(\omega)$, let's substitute $x = \cos(\omega)$ and use the quadratic formula:

$$(a^2 - 1) x + 2ax + (-4a - 2 a^2) = 0$$
$$x = \frac{-2a \pm \sqrt( 4a^2 - 4(a^2-1)(-4a - 2 a^2)) }{2 (a^2 - 1)}$$

Now, we want the solution corresponding to positive frequencies, and it turns out that choosing the result with the minus sign gives us that. So, we can take the inverse cosine, and get our -3 dB frequency:

$$\omega = \arccos\left( \frac{-2a - \sqrt{4a^2 - 4(a^2-1)(-4a - 2 a^2)}}{2 (a^2 - 1)} \right)$$

Let's plug it in for $a = -0.9$ and verify:

$w = 0.0958385$

$H(z) = 7.07107 = \dfrac{10}{\sqrt{2}}$

It might be useful to go the other way, to design a filter by choosing a value of a to have some desired -3 dB point:

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + a^2 cos^2(\omega) - \cos^2(\omega)$$

$$(2 - \cos^2(\omega)) a^2 + (4 - 2 \cos(\omega)) a + \cos^2(w) = 0$$

$$\frac{ -(4 - 2 \cos(\omega)) +/- \sqrt{ (4 - 2 \cos^2(\omega))^2 - 4(2 - \cos^2(\omega))(\cos^2(\omega)) } }{2 (2 - \cos^2(\omega))}$$

And again we can verify by plugging in $\omega = 0.0958385$, and we'll see that we get $a = -0.9$.

So, for example, you could use these formulas to pick values of $a$ at different sampling rates, so that you have the same frequency -3 dB point. You can use the formula for gain at DC, for example, to normalize the filter so its gain at DC is unity. You can take the formula for gain, convert it to decibels, and see that it approaches roughly a -6 dB/oct slope, which is why it's being used to shape noise like this.


trainCascadeObjectDetector seems to be missing from MATLAB R2012b

Looking at the release notes for the computer vision toolbox, it looks like that functionality was not added until the R2013a release.


Convolution Theorem: Hamming Window on a Time Series and Fourier Domain

Fistly, you need to pad the data (y) and the filter (Hamming window). The minimum required length to avoid circular convolution = data length + filter length - 1. In this case, it is 100 + 7 - 1 = 106. For convenient, you may choose power of 2 length, in this case 128.

So, pad the data to 128, do fft on it.

Pad the filter to 128, do fft on it.

Multiply both ffts.

Do ifft.


Approximation of Hilbert Transform Using Very Short Hilbert Transform FIR / IIR Filter

This is more like an extended comment to chart the possible answers.

Hilbert transform is a frequency domain 90-degree phase shift of the signal. It has an antisymmetrical impulse response around time = 0. You specify that the approximation shall be causal (EDIT: this requirement has since been removed), so I think you need to reference the phase shift to that provided by another filter. If we denote the two filters F1 and F2, some usable alternatives are:

Case | Filter F1    | Filter F2           | Phase diff     | Magnitude frequency responses
-----+--------------+---------------------+----------------+------------------------------
  C1 | pure delay   | antisymmetrical FIR | exactly 90 deg | F2 has ripple
  C2 | pure delay   | all-pass IIR        | approx 90 deg  | No ripple
  C3 | all-pass IIR | all-pass IIR        | approx 90 deg  | No ripple

In all those cases you have two filters whose phase frequency responses have a phase difference of exactly or close to 90 degrees over the desired frequency band.

If the filter outputs are used "in quadrature" (summed, with the other multiplied by the imaginary unit) to create an analytic signal, any phase frequency response ripple will become magnitude frequency response ripple of the composite negative frequency removal filter.


Is it necessary to apply some window method, to obtain the FFT? Java

As you can see; When a sinusoidal signal of arbitrary frequency is input to the ‘N’ point DFT, the spectral response of the signal will not be a Delta function (impulse) because when we take into account of the missing end point to be the same as the beginning of the next period of the periodic extension of the sequence, the periodic extension of the signal is no more a sinusoid. One way to solve this problem is to ensure Coherent Sampling as shown in figure 1

But coherent sampling can be done only if the DFT engine receives pure monotone signals. So in a real world application, Coherent sampling is not feasible.

SO, just diminish the magnitude of the discontinuities by multiplying the whole signal with some function.

So, when you multiply by window function, it eliminates the jump discontinuities and removes the false frequency that is introduced when coherent sampling is not done (see Figure 3.

The length of the window should be the length of your signal.
Selection of your window could be depending on how much side-lobe height do you tolerate.

Image Source: http://electronicdesign.com/analog/choose-right-fft-window-function-when-evaluating-precision-adcs


Let $X(f)$ denote the Fourier transform of $x(t)$ where
$$\begin{align*}
X(f) &= \int_{-\infty}^{\infty} x(t) \exp(-j2\pi ft) \mathrm dt\\
x(t) &= \int_{-\infty}^{\infty} X(f) \exp(+j2\pi ft) \mathrm df
\end{align*}$$
which I will denote via $x(t) \leftrightarrow X(f)$.
The following transform pairs will be needed in what follows.
$$\begin{align*}
\delta(t) &\leftrightarrow 1\\
\delta(t-t_0) &\leftrightarrow \exp(-j2\pi f \,t_0)\\
\sum_{n=-\infty}^{\infty}\delta(t-nT) &\leftrightarrow
\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{T}\right)\\
\end{align*}$$

Given a signal $x(t)$, its sampled pulse train (at intervals of $T$ seconds) is
$$x(t) \sum_{n=-\infty}^{\infty}\delta(t-nT)
= \sum_{n=-\infty}^{\infty} x(t)\delta(t-nT)
= \sum_{n=-\infty}^{\infty} x(nT)\delta(t-nT).$$
Since multiplication in the time domain corresponds to convolution in the
frequency domain, we have
$$\begin{align*}
x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT) &\leftrightarrow
X(f)\circledast\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{T}\right)\\
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
X(f)\circledast\delta\left(f-\frac{k}{T}\right)\\
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
\int_{-\infty}^{\infty} X(f-w)\delta\left(w-\frac{k}{T}\right) \mathrm dw\\
x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT)
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
X\left(f-\frac{k}{T}\right) = \hat{X}(f)
\end{align*}$$
Thus, the Fourier transform of the impulse train formed by
sampling $x(t)$ at $T$ second intervals is $\hat{X}(f)$ which
is obtained by repeating $X(f)$ along the $f$ axis at intervals of
$T^{-1}$ Hz and summing the result.
Furthermore, $\hat{X}(f)$ is a periodic function of the frequency
variable $f$ with period $T^{-1}$ Hz.  That is, for all $f$,
$$\hat{X}\left(f + \frac{1}{T}\right) = \hat{X}(f).$$  Note that
all of this holds regardless of what $X(f)$ is: $X(f)$ could be
nonzero for all $f$, and the result would still be valid.

Now suppose that $X(f)$ is zero if $f < a$ or $f > b$.  Then,
$\hat{X}(f)$, obtained by repeating $X(f)$ periodically along
the $f$ axis,  is nonzero only in the intervals
$$\ldots, \left[a-\frac{2}{T}, b-\frac{2}{T}\right],
\left[a-\frac{1}{T}, b-\frac{1}{T}\right],
[a, b]
\left[a+\frac{1}{T}, b+\frac{1}{T}\right],
\left[a+\frac{2}{T}, b+\frac{2}{T}\right],
\ldots
$$
and so if
$$b -\frac{1}{T} < a \Rightarrow b-a < \frac{1}{T},$$
that is, the support $b-a$ of $X(f)$ is smaller than
the repetition interval $T^{-1}$ Hz, then the repetitions
of $X(f)$ do not overlap.  Indeed as the OP's figures show,
for a real-valued signal whose spectrum extends from $a = -f_0$
to $b = f_0$ (support of $X(f)$ is of length $2f_0$), sampling
$x(t)$ at intervals of $T = (3f_0)^{-1}$ (and thus repeating $X(f)$ at
intervals of $3f_0$ Hz on the frequency axis)
leads to no overlap while sampling $x(t)$ at intervals of
$T = (1.5f_0)^{-1}$ leads to overlap of the repetitions
of $X(f)$.

Finally, the OP asks about negative frequencies and the
interpretation of the pictures which show positive frequencies
only.  For a real-valued signal, $X(f)$ has conjugate symmetry,
meaning that $X(-f) = X^*(f)$, and so specifying $X(f)$ for
positive values of $f$ suffices.  In any case, the pictures
he is looking at are of $|X(f)|$ and $|\hat{X}(f)|$ which
are even functions of $f$, and so showing only the positive
axis saves space, though it does make the pictures look a
bit lopsided since only half the lobe is shown at low frequencies.
For complex-valued signals, $X(f)$ does not have conjugate
symmetry and $|X(f)|$ need not be an even function of $f$
and so the whole axis would need to be shown.  But, the general
development above is still applicable, and we still need
to sample at rate exceeding $b-a$ Hz, and it helps to keep in
mind that in this general case, each sample is actually two real
numbers, not one, since we are sampling a complex-valued signal.


Valid approach to implementing a band pass filter

then you'll be doing either "Overlap-add" (OLA) or "Overlap-save" (OLS) which is sometimes more clearly called "Overlap-scrap".  you need to read up about this.

then simply "apply[ing] gain in those bins" is multiplication in the frequency domain, which means convolution in the time domain.  if your frequency domain function is a rectangle function (some number for "those bins" and another number, like 1, for the other bins).  then you need to compute the impulse response that goes with it and make sure that FIR impulse response is not too long for your OLA or OLS.  do you have a good book that describes OLA or OLS?

but a bandpass filter in the time-domain (your option 1.) is simple and can be implemented simply.  maybe more difficult to design your IIR or FIR bandpass filter.