Calculating the output of a pole eigen signal in a difference equation

If the input is $x[n]=\left(-\frac12\right)^n$ for all $n$, then the output is infinite, because the convolution sum doesn't converge. Note that also the input grows without bounds for negative $n$.
If, on the other hand, the input starts at $n=0$, i.e., $x[n]=\left(-\frac12\right)^nu[n]$ then you can use the $\mathcal{Z}$-transform to compute the output:
$$Y(z)=X(z)H(z)\tag{1}$$
With $$X(z)=\frac{1}{1+\frac12 z^{-1}}\tag{2}$$
you obtain
$$\begin{align}Y(z)&=\frac{1+3z^{-1}}{\left(1+\frac12 z^{-1}\right)^2\left(1-\frac12 z^{-1}\right)}\\&=\frac{a}{\left(1+\frac12 z^{-1}\right)^2}+\frac{b
}{1+\frac12 z^{-1}}+\frac{c}{1-\frac12 z^{-1}}\tag{3}\end{align}$$
with some real-valued constants $a$, $b$, and $c$.
Inverse $\mathcal{Z}$-transform of $(3)$ gives
$$y[n]=-2a\;n\left(-\frac12\right)^n+b\left(-\frac12\right)^n+c\left(\frac12\right)^n\tag{4}$$

EDIT: First, to clarify, I use a basic example that shows that a double pole in the $\mathcal{Z}$-transform is no problem. Assume a system with impulse response $h[n]=(1/2)^nu[n]$ and an input signal $x[n]=(1/2)^nu[n]$. Computing the output in the time domain is straightforward:
$$\begin{align}y[n]&=u[n]\sum_{k=0}^n\left(\frac12\right)^k\left(\frac12\right)^{n-k}\\&=u[n]\left(\frac12\right)^n\sum_{k=0}^n1\\&=(n+1)\left(\frac12\right)^nu[n]\tag{5}\end{align}$$
The $\mathcal{Z}$-transform of the output sequence is simply given by the multiplication of the $\mathcal{Z}$-transforms $H(z)$ and $X(z)$:
$$Y(z)=\frac{1}{\left(1-\frac12z^{-1}\right)^2}\tag{6}$$
The inverse $\mathcal{Z}$-transform of $(6)$ equals the result $(5)$, which is very easy to confirm.

The probable cause of the OP's confusion could be the fact that for an input signal $x[n]=z_0^n$, $-\infty<n<\infty$, the output is given by
$$y[n]=H(z_0)z_0^n$$
if and only if $z_0$ is inside the ROC of $H(z)$, because otherwise the convolution sum doesn't converge and the expression $H(z_0)$ is meaningless. However, in the computation of the outputs $y[n]$ in the examples above we never evaluate $H(z)$ for any $z$ outside the ROC, so the situation is entirely different.


magnitude-only FFT can really only be done with the complex FFT followed by computing the magnitude.

the radix-2 Cooley-Tukey FFT (with the butterflies) will likely be the easiest to implement with a FPGA because of the self-similarity of the butterflies in each FFT pass.  (if you want some simple C code that you might translate into FPGA code, i can supply you with a very simple radix-2 FFT in C.  are you limited to fixed-point or can you do floating point in your implementation?)

if you're displaying spectrum magnitude in dB, then magnitude-squared (without the square root required for magnitude) is easier and just changes the scaling factor on the $\log(\cdot)$ function.

$$\begin{align} 20\log_{10}\left( \Big|X[k] \Big|\right) &= 10\log_{10}\left( \Big|X[k] \Big|^2\right) \\ \\ &= 10\log_{10}\left( \Re\{X[k]\}^2 + \Im\{X[k]\}^2 \right) \end{align} $$

so, for dB, don't bother with the square-root.

if your spectrum analyzer also wants log frequency along the frequency axis, that's another different problem to solve.


How to simplify the Fourier Series Using an Approximation?

Sounds like bandlimited synthesis.

  with the added ability to change the lowpass cutoff frequency of the signal dynamically.

That makes it even more complicated.

The simplest way to do this right is additive synthesis, as you know.   That way you can adjust each harmonic's frequency and phase and amplitude on the fly. Yes, there are a lot of calculations if the fundamental frequency is low, but you can use approximations of sin() if that speeds things up.  The 2-MAC resonant filter/oscillator (1, 2), for instance:

or various polynomial approximations with varying amounts of distortion:

If the waveform is fixed and a submultiple of the sample rate, you could use IFFT synthesis, but for anything else I imagine this would be difficult because the harmonics cover multiple FFT bins.  Maybe with overlapping and fading between chunks?

Otherwise the best solution depends on what exact waveform you're trying to synthesize.

If you only need a "sawtooth-like" waveform, for instance, you can use the DSF directly, which produces waveforms with harmonics that fall off linearly:

but if you want an actual square wave you use the DSF to generate a "bipolar" band-limited impulse train (BLIT):

and then leaky-integrate the sum. 2

These methods can't be used for arbitrary waveforms, however, and can't be used for swept-frequency chirps or swept-frequency LPFs (as you want?), since the harmonics suddenly disappearing from one sample to the next produces a pop or click which can't be filtered out after the fact.

You could probably use this method for the low harmonics and then combine with additive synthesis of the top harmonics and then fade those in and out to prevent the clicks.  In other words, as a harmonic disappears from one sample, you seamlessly replace it with an additive-synthesis harmonic at the same frequency and phase, which you can then fade out more smoothly.

Wavetable synthesis, fading between wavetables as harmonics disappear?
Then there are other methods.  :)

https://ccrma.stanford.edu/~stilti/papers/blit.pdf


Inverse DTFT of $H_1(\Omega)=\begin{cases} 10,& \frac{\pi}{3} \leq |\Omega| < \pi\\ 0,& 0 \leq |\Omega| < \frac{\pi}{3}\\ \end{cases}$

Let's use $\omega$ instead of $\Omega$ for the definition of the inverse DTFT and state the following:

$$h[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} H(\omega) e^{j \omega n} d\omega $$ Given that $H(\omega) = 1 $ , for $\pi/3 < |\omega| <\pi $ and zero otherwise, the integral becomes:

$$\begin{align}
h[n] &= \frac{1}{2\pi} \int_{-\pi}^{-\pi/3} e^{j \omega n} d\omega + \frac{1}{2\pi} \int_{\pi/3}^{\pi} e^{j \omega n} d\omega \\
&= \frac{1}{2\pi} \left(  \frac{ e^{j \omega n} }{j n} |_{-\pi}^{-\pi/3} + \frac{ e^{j \omega n} }{j n} |_{\pi}^{\pi/3}  \right)\\
&= \frac{1}{2\pi} \left(  \frac{ e^{-j (\pi/3) n} -e^{-j\pi n}  }{j n} + \frac{ e^{j \pi n} - e^{j \pi/3 n} }{j n}  \right)\\
&= \frac{1}{n\pi} \left(  \frac{ e^{-j (\pi/3) n}- e^{j (\pi/3) n} }{2j} + \frac{ e^{j \pi n} - e^{-j\pi n} }{2j}  \right)\\
&= \frac{1}{n\pi} \left(  -\sin(\frac{\pi}{3}n) + \sin(\pi n) \right)\\
&= \frac{\sin(\pi n) }{n\pi} - \frac{ \sin(\frac{\pi}{3}n)}{n\pi}  \\
\end{align}
$$

at this point recognize  that
$$\frac{\sin(\pi n) }{n\pi} = \begin{cases} 1 &,\text{ for } n=0 \\0 &, \text{ for } n \neq 0 \\ \end{cases} = \delta[n]$$

Hence the I-DFT becomes:
$$ \boxed{ h[n] =  \delta[n] - \frac{\sin( \frac{\pi}{3}n) }{\pi n} }$$

Note that your $H(\omega)$ had a gain of $10$ instead of $1$ so multiply this with $10$.


I do find this chapter quite well-written. The row or the column visions of a vector are mostly matters of convention and convenience. One or the other is perfectly fine.   In signal processing, when using algebra, the column version (a vertical vector) is used quite often. If you multiply a col-vector $x$ by matrix $A$, you get a col-vector $y$:

$$ y = Ax$$

You can get the row-vector version by transposing:

$$ y' = x'A'$$

If you decide to work with row-vectors, you could as well work with $v=uB$, but the matrix on the left seems more easy to read (at least with western languages).

When you deal with an image or image block $X$, you can apply a transform on both the columns, at $Y_1=AX$. $Y_1$ is a matrix whose columns are the transformation of each column of $X$. Then, you can work on the rows of $Y_1$, by applying the formula  $Y_1 A'$: every row is now affect.ed by the transform. The final result is thus $(AX)A'$, or $A((AX')') = A(XA')$


Measuring time distance between bass tones and send it as a signal

You don't exactly need the whole Fast Fourier Transform (FFT) "mechanism" to achieve what you are trying to do. Also, your biggest problem in this setup is going to be synchronisation. Here is why:

You are looking at something like this:

Source --> Low Pass Filter (Extracting the bass) --> Rectification --> Integration --> Threshold --> Motor

You can implement this setup with discrete convolution, which means that you can have this pipeline operating on a per sample basis. A new sample comes in, it is processed and as a result a new processed sample pops out of that pipeline.

At its most basic and simplest form, a digital filter is a set of coefficients that are multiplied with the different time instances of the incoming signal in a rolling window fashion. Therefore, the "problem" now is how do you calculate those coefficients. There is a very large number of methods to achieve this but the simplest way to design a low pass filter is a sinc windowed filter. For more information, please see this link.

So, now you have your filter coefficients. Let's say you have designed a low pass filter with 50 coefficients and a cut off frequency of 100Hz. The filter is returning a signal that contains ONLY the first 100Hz of the Source.

How do you use this to drive the motor?

If you are not doing Pulse Width Modulation (PWM) on the motor, to create different intensities of vibration, then all that you are after is an "on-off" signal.

This "on-off" signal can be produced by a threhsold function. In other words, if the signal becomes "louder" than some pre-set value, then turn the motor on, otherwise, keep the motor turned off.

To judge if the low pass filtered signal is loud enough, you need to first rectify it (take its absolute value) and then integrate it (apply yet another low pass filter to smooth the data). This is equivalent to classical rectification and smoothing out. The low pass filter for smoothing here can be a sinc filter as above, with a very low cut off frequency (so that it becomes an integrator). Alternatively, you could use a moving average filter. For a given length, say $N$, the moving average filter's coefficients are $h_i = \frac{1}{N}, i \in [0,1, \ldots, N-1] $.

At the end of this process, you will have a "pulse" that goes high when the region below 100Hz is "loud" and low otherwise.

Finally, apply a threshold to this pulse (for example $x>T$ where $T$ is some value in the range of the signal and you have your square pulse that drives the motor.

You now have to put this "on-off" signal in a tiny little data packet and send it via serial (or other way) to the Arduino which will be doing the hardware bit.

In fact, if you were to use fixed point arithmetic, you can handle part of this process entirely in the Arduino. Unfortunately, they don't go "audio quality" high in sampling so that you could carry out everything in the Arduino.

Your biggest problem here is going to be synchronisation. Due to the way that digital filters work, there is a delay between the input and the output. The longer your (FIR) filters are, the better filtering they will be doing but also, the longer the delay they will be introducing to your signal. On top of this, you are going to have to add the delay of the connection with the arduino and the "packet encoding and decoding" that has to take place anyway.

In other words, your "kick" will be occuring slightly later than the bass and this might be annoying, depending on your specifications. You would have a similar issue even if you were implementing this with analog filters.

If you want the most absolutely low delay between input and output then you are going to have to move towards more complex to design filters.

Hope this helps.

EDIT:

To drive the motor with PWM and produce different levels of vibration for different levels of "bass" in the signal, you need to substitute the stages of "recitication --> integration" with a simple "quantisation" stage.

Quantisation is simply a linear transformation from the range of possible values of your "bass" signal (say 0-255) to the range of possible PWM levels (say 0-8).

The fact that the signal can have practically infinite dynamic range (AS COMPARED TO THE VIBRATION LEVELS) does not mean that the PWM stage can produce infinite vibration levels with the differences between them being discriminable.

Just as it happens with joystick controllers, the levels of vibration are usually very small, for example 4 or 8. There is no point in driving the motor in one instance with 2.812721 ms pulses and in the next with 2.812728ms pulses. The difference is incredibly small to be felt.

Quantisation can be as simple as round(bassLevel/maxBassLevel*vibrationLevels) to produce a number between 0 and vibrationLevels which you can then map to PWM periods through a lookup table. For more information on quantisation, please see this link

So, now your pipeline is:
Source --> Low Pass Filter (Extracting the bass) --> Quantisation --> Motor

And this will be relatively faster in responding because you got rid of one low pass filter.


Kalman filtering for position using GPS,accelerometer and speed sensors

In my experience at the very least you need a nine or ten state Kalman (velocity, position, and angle -- that's nine states if you attempt to use Euler angles, ten if you use quaternions).  Practically, you need to use IMU offset, so add six more states.  If you use the IMU's compass function, add three more.  It may not be a bad idea to monitor the acceleration due to gravity -- one more.  If you use the vehicle's odometer function (from CAN), then add a state for proportionality error, and possibly offset error.
You need vehicle orientation.  That's because the ground is accelerating upwards at $9.81 \mathrm{m/s^2}$ (if Einstein is to be believed), or the accelerometer's proof mass is being pulled down with a force equal to $9.81 \mathrm{m/s^2/kg}$ (if Sir Isaac Newton is to be believed).  Either way, to get the actual acceleration with respect to the ground you need to null out the acceleration with respect to gravity, and you can only do that if you have the vehicle's orientation in 3D.
It's been over a decade since I worked on this last, but there has to be canned packages available open source -- have you done a search?


How is multiplexing on FM signals performed?

In U.S. broadcast FM, the various subcarriers (pilot tone, stereo L-R, RDS, etc.) are modulated into the upper sidebands (above the 15 kHz L+R/mono audio, but below 100 kHz, for instance 38 kHz DSB for the stereo subcarrier) of the baseband signal before that entire multiplex is used to wideband frequency  modulate a single, much higher frequency, FM carrier.

The FM carrier is modulated, but not shifted in center frequency.


Programming Wavelets for Audio Identification

Convolution can be used to find the note. Its best shown visually .

  Would this be a data structure holding e.g {sampleNumber, frequency} pairs?

Yes, the signal would be. The wavelet can be stored as a function but it can also be {sampleNumber, frequency} pairs.

  If a wavelet is an array of these pairs, how is it applied to the audio data?

Convolution:

$f*g = \int f(x)g(t-x)dx$

Convolution is commutative so it does not matter whether f or g is wavelet or signal or vice versa.

  How does this wavelet apply to the analysis when using an FFT?

In some applications in signal processing, convolution is calculated using the convolution theorem.

$\mathcal{F}(f*g) = \mathcal{F}(f) \mathcal{F}(g)$

so the convolution of f and g can alternatively be calculated by

$f*g = \mathcal F ^{-1}[ \mathcal{F}(f) \mathcal{F}(g)]$

  What is actually being compared to identify the signal?

The signal is being compared to the wavelet

  How exactly is a wavelet used digitally?

There are many applications, besides note identification, they can be used in lossy signal compression. See the Dirac codec.

Cross Correlation is a similar operation which can be used for identifying wavelets


Intensity scaling can be just as you say: multiplying all intensities by a scalar value. This will obviously change the histogram, but will do nothing like histogram equalization (it will just scale the x-axis of the histogram, barring some binning effects).

More generally, intensity scaling can change the range of intensities so that instead of going from $\tt x_{min}$ to $\tt x_{max}$ the pixel values range from $\tt x_{new\ min}$ to $\tt x_{new\ max}$ that might be better for display purposes (e.g. $\tt x_{new\ min} = 0$ to $\tt x_{new\ max} = 255$ or $\tt x_{new\ min} = 0$ to $\tt x_{new\ max} = 65535$).

See Gozales, Woods, and Eddis for example.


Anisotropic Diffusion Filter - Intuition Behind Parameters

This paper - The Structure of Images describes quite nicely to get to the diffusion equations from the point of view of a filter. It's about the isotropic case, but anisotropic is really just an extension of that. For anisotropic diffusion see Scale Space and Edge Detection Using Anisotropic Diffusion as Deniz suggests this one is quite good as it also includes discretised versions of the equations.

The basic idea is to solve the heat equation where $t$ is the scale of the image, i.e. the amount of filtering to apply. For anisotropic diffusion the diffusion coefficient $c$ is replaced by $g$ which is dependent on the image gradient and therefore filters preferentially depending on the gradient.

The divergence here really just arises from the definition of the heat equation.


Why does signal averaging reduces noise levels by more than $\sqrt{n}$?

Well, I would say the assumption that your noise is Gaussian is ill fitting.  If the noise is due to machine interference, it probably has some tonal characteristics.  Tones of the same frequency can reinforce or cancel each other out when added.

To get a better idea of what might be going on, you should:

1) Make a histogram of the noise

2) Take an FFT of the noise

The histogram should look like a bell curve and the FFT results should be flat.  If they aren't, your noise is not Gaussian.

Hope this helps.

Ced


Power Spectrum and Power Spectral Density

What is the link between them? Power should not be function of
  frequency (it is the integral on all frequencies of power spectral
  density)

As far as a Spectrum Analyzer is concerned, to display the power spectrum in dBm, what it does it to take the FFT and compute square of magnitude $P_n = 20  \times log_{10}(|X(f_n)|^2*1000/100)$ with $f_n$ being the discrete frequencies separated by resolution bandwidth (see it is set to 1kHz for this case). I multiplied by 1000 since unit is in dBm and I divided by 100 because to get to power, I have to pass voltage $X(f_n)$ through a resistance of $50\Omega$. ($\frac{V^2}{2R}$).

To derive power spectral density from power spectrum, I am sure it is just dividing power in each bin by RBW. $P_n/f_{rbw}$. Though I am not very sure whether it is the (dBm value divided by rbw) or (dBm value of linear power value divided by rbw).

  I think this representation has no sense: the ordinate should be the
  power spectral density, not its integral in a neighborhood around the
  frequency considered

Well, in this case it is the latter. I think that is how all commercial spectrum analyzers represent spectrum (if you choose units as dBm or Watts). The comment by MBaz correctly mentions this.

Note - The final plot will be slightly different because your FFT size output will be different from number of trace points used to display on the screen. In order to combine several FFT points into 1 trace point on screen, a detector logic is used.


Impulse Response from Frequency Response at Specific Values


How can I verify that a code belongs to a group?

Assuming these are really just sub-words that are common to all elements in a group, simply trying out all 8 possible (cyclic) bitshifts of the received sequence to see whether you can match any group representative (e.g. the first column in your GroupN table) would very likely be the most efficient approach – because 8 bit is really not that much, and bit shifts are very efficient on CPUs or on FPGAs (whatever you're using), and even more importantly (since "beacons" surely doesn't sound like high data throughput), the code is easy to understand, and premature optimization is the root of all evil.


Morlet wavelet time and frequency resolution

Consider shifted and scaled versions of a mother wavelet $\psi(t)$:

$$\psi_{a,b}(t)=\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right),\quad a>0,\;b\in\mathbb{R}\tag{1}$$

By the definition of the Fourier transform

$$\Psi(\omega)=\int_{-\infty}^{\infty}\psi(t)e^{-i\omega t}dt$$

it can be shown that the Fourier transform of $\psi_{a,b}(t)$ is

$$\Psi_{a,b}(\omega)=\sqrt{a}e^{-ib\omega}\Psi(a\omega)\tag{2}$$

Note that the scaling factor $a$ appears in the denominator of the argument of $\psi_{a,b}(t)$ and in the "numerator" of the argument of $\Psi_{a,b}(\omega)$. This implies that stretching in the time domain ($a>1$) implies compression in the frequency domain and vice versa. Consequently, if $\Delta t$ and $\Delta f$ denote the "spread" of $\psi(t)$ in time and frequency, it follows from (1) and (2) that the corresponding spreads of $\psi_{a,b}(t)$ - call them $\Delta t_{a}$ and $\Delta f_{a}$ - are given by

$$\Delta t_a=a\Delta t\quad\textrm{and}\quad\Delta f_{a}=\frac{\Delta f}{a}$$

As for the definition of $\Delta t$ and $\Delta f$, there is not way to "prove" the results that you stated. It is just a matter of defining the spread of a function. Since $|\psi(t)|$ is a Gaussian with standard deviation $\sigma=\sqrt{f_b/2}$ it seems natural to define $\Delta t$ as being proportional to the standard deviation, but there is no "correct" way to choose the proportionality constant. Once $\Delta t$ has been defined, the corresponding definition of $\Delta f$ can be obtained from the uncertainty principle.


How do I find the power of a specific frequency at a specific sample of a discrete time varying signal?

I realize this doesn't answer the detail of what you're asking, but I believe it answers the question in the title:

See Rick Lyons's writeup here for using the Goertzel algorithm to get the signal level at a specific frequency.

You can then follow up the equations:

with:

$$
\hat{p}(n) = | y(n) |^2
$$

and use $\hat{p}(n)$ as your "instantaneous power" at frequency $\omega =2\pi m/N$ --- which is sort of like  Lyons's eq. 13.83:

(though he's evaluating this one at the "end" of his DFT length).


How to find the difference equation directly from Direct Form II signal flow graph

Darkling, these things are quite clearly explained in standard signals & systems textbooks. But I assume you have little time to read (as most undergrad courseware are full of homeworks, preliminaries, projects that most of which I believe become harmful as they consume (more than necessary amount of) a student's learning time and just completing a homework or project without full understanding makes what sense..? anyway)

Consider the LCCDE associated with a causal LTI system (such a system is said to be at initial rest which assures that all the initial conditions ard zero; i.e., no output occurs before an input is applied) :

$$ \sum_{k=0}^{N} a_k ~y[n-k] = \sum_{k=0}^{M} b_k ~ x[n-k] $$ The order is $\max\{M,N\}$.

Now the LCCDE can be considered to be a cascade of two LTI systems; the RHS represents an IIR-LTI system represented by $A(z)$ and LHS represents an FIR-LTI system represented by  $B(z)$. So the system can be represented with the following block diagrams:

$$ x[n] \rightarrow \boxed{ B(z) } \rightarrow v[n] \rightarrow \boxed{ \frac{1}{A(z)} }\rightarrow y[n] \tag{1}$$

or equivalently
$$ x[n] \rightarrow \boxed{ \frac{1}{A(z)} } \rightarrow v[n] \rightarrow \boxed{ B(z) } \rightarrow y[n]  \tag{2}$$

This follows from the fact that both cascaded systems are LTI and order of LTI systems can be interchanged which follows from the fact that convolution operator is commutative and associative; i.e., $$h_1[n] \star h_2[n] \star x[n] = h_2[n] \star h_1[n] \star x[n] $$

Now, the block diagram of (1) can be implemented with the following signal flow graph which yields a Direct-Form-I structure :

And, the block diagram of (2) can be implemented with the following signal flow graph which yields a Direct-Form-II structure:

The intermediate memory (delays) are merged and the signal flow graph is simplified into the canonical direct-form II form finally as (assuming $N > M$)

From this last figure, it can be deduced that your LCCDE relating the input $x[n]$to the output $y[n]$ is :

$$ y[n] - \frac{3}{5} y[n-1] + \frac{38}{75} y[n-2] + \frac{2}{15} y[n-3] = x[n] - \frac{3}{10} x[n-1] + \frac{1}{3} x[n-2] $$


According to this IEEE article: https://github.com/liulinbo/slam/blob/master/A%20robust%20and%20easy%20to%20implement%20method%20for%20IMU%20calibration%20without%20external%20equipments.pdf

You can model errors this way:
$$ a = f*a' + g + b + \eta$$
where $a$ is the actual acceleration, $f$ is a 3x3 matrix to model scaling, misalignments, cross-axis and ... errors. $a'$ is sensor's data, $g$ is gravity, $b$ is 3x1 matrix to model bias, and $\eta$ is 3x1 matrix to model noise.

You can calibrate accelerometer by reading $a'$ from sensor data when the accelerometer is in static positions. In static positions, the only force effecting accelerometer is gravity. So by minimizing this summation, $f$ and $b$ can be computed(calibration):
$$ (|a_0|^2 - |g|^2)^2 + (|a_1|^2 - |g|^2)^2 + ... + (|a_N|^2 - |g|^2)^2 $$
where $N$ is number of static positions.


Where has this complex component come from?

At a quick glance, it looks like you're right. The impulse with an area of $1$ at zero frequency respresents a DC offset of 1. The $0.25 - j0.43$ complex value represents the amplitude and phase of the cosine at 2.5 Hz.

Specifically, the two frequency domain impulses at +/- 2.5 Hz correspond to time-domain complex exponentials:

$$
\left(0.25 - j0.43\right) \delta\left(f - 2.5\text{ Hz}\right) + \left(0.25 + j0.43\right) \delta\left(f + 2.5\text{ Hz}\right) \Leftrightarrow \\ \left(0.25 - j0.43\right)e^{j2\pi 2.5 t} + \left(0.25 + j0.43\right)e^{-j2\pi 2.5 t}
$$

(This can be shown from a table of common Fourier transform pairs, like this one)

Next, let's rewrite the $(0.25 \pm j0.43)$ factors in polar form:

$$
(0.25 \pm j0.43) = 0.50 e^{\pm j2\pi0.17}
$$

The right-hand side of the equation then becomes:

$$
0.50 e^{j2\pi0.17}e^{j2\pi 2.5 t} + 0.50 e^{-j2\pi0.17}e^{-j2\pi 2.5 t}
$$

$$
0.50 e^{j2\pi (2.5 t+0.17)} + 0.50 e^{-j2\pi (2.5 t+0.17)}
$$

We can then apply a form of Euler's identity:

$$
\cos(x) = \frac{e^{jx}+e^{-jx}}{2}
$$

$$
0.50 e^{j2\pi (2.5 t+0.17)} + 0.50 e^{-j2\pi (2.5 t+0.17)} = 1.0 \cos\left(2\pi (2.5t+0.17)\right)
$$

So, your signal should have a DC component of 1, plus a cosine at frequency 2.5 Hz with amplitude of approximately 1 and a phase offset of 0.17 cycles. From inspection, it appears that figure R1 is close to that form.


Does flat surface detection and orientation require special hardware to compute?

Modern augmented reality platforms such as Google's ARCore and Apple's ARKit seem to only operate on mobile devices,...

They only make sense on such devices.

  ... I'm guessing these libraries would never be able to work off of mobile platforms.

They do work outside of a mobile device... (and here).

  Flat surface detection inside an image or video

For anything to be detected by a computer vision algorithm it has to be in contrast to its environment. If not, then extraneous information would have to be used that does make contrast, e.g. obtain images in different wavelengths.

Provided that surfaces are in contrast to their surroundings it is possible to reconstruct 3D geometry (e.g. multiple surfaces coming together forming one object), by obtaining more than one images around the object of interest.

If a surface is flat, its orientation can also be be detected easily with a marker anyway.

  Orientation/angle of the camera view inside an image or video

This is called "Simultaneous Localisation and Mapping (SLAM)", it can be done with or without a marker (usually, when the camera is used as a sensor in robotics, it is done without a marker) and there are lots of libraries already on places like github that you can look for inspiration or further details of their inner workings.

Hope this helps.


Gradient of Total Variation (TV) Norm in Total Variation Denoising

I am by no means an expert on total variation, however I think you should check out this Wikipedia page. It doesn't directly answer your question, but I believe the lemma below illustrates the relationship between total variation and divergence.

There, it gives a lemma that follows from the Gauss-Ostrogradsky theorem and provides a proof for it,

$\int_{\Omega}f\,div(\phi) = -\int_{\Omega} \nabla f \cdot\phi$.

I think for some intution about divergence, it might be helpful to read Wikipedia's explanation.


The pre-envelope is also called analytic signal. Its Fourier transform is indeed given by the expression in your question:

$$X_+(f)=2X(f)u(f)\tag{1}$$

where $X(f)$ is the Fourier transform of the original signal, and $u(f)$ is the unit step function. Obviously, $X_+(f)$ has only positive frequency components. The analytic signal $x_+(t)$ with its Fourier transform given by $(1)$ is necessarily a complex-valued signal:

$$x_+(t)=x(t)+j\mathcal{H}\{x(t)\}\tag{2}$$

where $\mathcal{H}\{\cdot\}$ denotes the Hilbert transform. Note that the analytic signal is not the same as the complex envelope. For a band pass signal $x(t)$, the complex envelope is a low pass signal, whereas the analytic signal is a band pass signal.

From $(2)$, in order to compute the analytic signal you need to compute the Hilbert transform of $x(t)$:

$$\hat{x}(t)=\mathcal{H}\{x(t)\}=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\tag{3}$$

For the given signal you get

$$\hat{x}(t)=\frac{1}{\pi}\int_{-a}^a\frac{\cos(2\pi f_0\tau)}{t-\tau}d\tau\tag{4}$$

The result of $(4)$ can be written in terms of the cosine integral $\text{Ci}(x)$ and the sine integral $\text{Si}(x)$:

$$\hat{x}(t)=\frac{1}{\pi}\left\{\cos(2\pi f_0t)\left[\text{Ci}(2\pi f_0(t+a))-\text{Ci}(2\pi f_0(t-a))\right]+\\\sin(2\pi f_0)\left[\text{Si}(2\pi f_0(t+a))-\text{Si}(2\pi f_0(t-a))\right]\right\}\tag{5}$$

I don't think that anybody expected you to come up with that awful expression. Anyway, for large values of $f_0$, the first term in $(5)$ becomes very small, and the second term converges to $\sin(2\pi f_0)\Pi_a(t)$. So for large $f_0$ you get the expected result

$$x_+(t)\approx x(t)+j\sin(2\pi f_0)\Pi_a(t)=\Pi_a(t)e^{j2\pi f_0 t}\tag{6}$$

Obviously, the respresentation

$$x(t)=\text{Re}\left\{\Pi_a(t)e^{j2\pi f_0 t}\right\}\tag{7}$$

is always valid, but the complex-valued signal $\Pi_a(t)e^{j2\pi f_0 t}$ is no analytic signal, it's just a good approximation of the analytic signal for large values of $f_0$.

Note that for $x(t)=m(t)\cos(2\pi f_0t)$ with $m(t)$ a band-limited function, i.e. $M(f)=0$ for $|f|>B$, the complex-valued signal $m(t)e^{j2\pi f_0 t}$ is an analytic signal, as long as $f_0>B$. The problem with the function given in your question is that $\Pi_a(t)$ is not band-limited.


Relationship between order of filter and amplitude/phase response

Does a higher order simply mean a better approximation to the optimal filter with (cutoff, stopband) ?

That depends on how you define "optimal". In general a high order will give you a steeper rolloff

  Do higher orders mean different or less severe phase differences?

Higher order will result in more phase distortion

  Finally, I'd assume that it affects performance somewhat as well.

That depends on how you define "performance". The "best" filter always depends on the application.

In general higher IIR filter order will give you more steepness in the frequency but also more time smear and phase distortion. Any beginner class on digital filter should cover this. For example, take a look at https://ccrma.stanford.edu/~jos/filters/


What is the difference between Spectrogram and Scalogram?

Both convert a 1D signal into a 2D complex map, aiming at unveiling transient features of the data. Both use a functional microscope: the signal $s(t)$ is analyzed through variations of a single template function $\rho(t)$ with two parameters $p$ and $q$.

So the signal is analyzed through the form:

$$R_s(p,q) = \int s(t) \rho^*_{p,q}(t) d\cdot$$

which is 2D (the differential increment $d\cdot$ can be specific).

With the spectrogram, $\rho$ is a fixed (non-zero sum) window, and $\rho_{p,q}$ is the same window, shifted by $p$, and  modulated by a frequency proportional to $q$.
With the scalogram, $\rho$ is a fixed zero-mean function (a wavelet), and $\rho_{p,q}$ is the same wavelet, shifted by $p$, and  dilated by a factor proportional to $q$.

Why they have those properties is another story.


The phase term of your outermost cosine function is time-varying: $\phi(t):=t a (1+b \cos(ct))$. The instantaneous frequency is given by $\frac{\partial \phi}{\partial t} = a+ab \cos(ct) - abct \sin(ct).$ In particular, in a tiny neighborhood around $t=0$ we can argue that $f(t)$ has only one frequency component $a+ab$. At large values of $t$, the $abct\sin(ct)$ term takes control.

Notice that the frequency varies with time. This makes it difficult to intuitively interpret the Fourier transform which assumes that all frequencies exist for all times. For such signals, it may be better to use time-frequency analysis by looking at the spectra of windowed segments of $f(t)$ around different values of $t$.


You're not looking for edges (=borders between extended areas of high and low gray value), you're looking for ridges (thin lines darker or brighter than their neighborhood), so edge filters might not be ideal: An edge filter will give you two flanks (one on each side of the line) and a low response in the middle of the line:

ADD: If've been asked to explain the difference between an edge detector and a ridge detector more clearly. I apologize in advance if this answer is getting very long.

An edge detector is (usually) a first derivative operator: If you imagine the input image as a 3D landscape, an edge detector measures the steepness of the slope at each point of that landscape:

If you want to detect the border of an extended bright or dark region, this is just fine. But for the veins in the OP's image it will give you just the same: the outlines left and right of each vein:

That also explains the "double line pattern" in the Canny edge detector results:

So, how do you detect these thin lines (i.e. ridges), then? The idea is that the pixel values can be (locally) approximated by a 2nd order polynomial, i.e. if the image function is $g$, then for small values of $x$ and $y$:

$g(x,y)\approx \frac{1}{2} x^2 \frac{\partial ^2g}{\partial x^2}+x y \frac{\partial ^2g}{\partial x\, \partial y}+\frac{1}{2} y^2 \frac{\partial ^2g}{\partial y\, ^2}+x \frac{\partial g}{\partial x}+y \frac{\partial g}{\partial y}+g(0,0)$

or, in matrix form:

$g(x,y)\approx \frac{1}{2} \left(
\begin{array}{c}
 x & y
\end{array}
\right).\left(
\begin{array}{cc}
 \frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
 \frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right).\left(
\begin{array}{cc}
 x \\
 y
\end{array}
\right)+\left(
\begin{array}{cc}
 \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{array}
\right).\left(
\begin{array}{c}
 x \\
 y
\end{array}
\right)+g(0,0)$

The second order derivative matrix $\left(
\begin{array}{cc}
 \frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
 \frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right)$ is called the "Hessian matrix". It describes the 2nd order structure we're interested in.

The 2nd order part of this function can be transformed into the sum of two parabolas $\lambda _1 x^2 + \lambda _2 y^2$ rotated by some angle, by decomposing the Hessian matrix above to a rotation times a diagonal matrix of it's eigenvalues (Matrix decomposition). We don't care about the rotation (we want to detect ridges in any orientation), so we're only interested in $\lambda _1$ and $\lambda _2$

What kind of shapes can this function approximation have? Actually, not that many:

To detect ridges, we want to find areas in the image that look like the last of the plots above, so we're looking for areas where the major eigenvalue of the Hessian is large (compared to the minor eigenvalue). The simplest way to detect that is just to calculate the major eigenvalue at each pixel - and that's what the ridge filter below does.

A ridge filter will probably give better results. I've tried Mathematica's built in RidgeFilter (which calculates the major eigenvalue of the Hessian matrix at each pixel) on your image:

As you can see, there's only a single peak for every thin dark line. Binarizing and skeletonizing yields:

After pruning the skeleton and removing small components (noise) from the image, I get this final skeleton:

Full Mathematica code:

ridges = RidgeFilter[ColorNegate@src];
skeleton = SkeletonTransform[Binarize[ridges, 0.007]];
DeleteSmallComponents[Pruning[skeleton, 50], 50]

ADD:

I'm not a Matlab expert, I don't know if it has a built in ridge filter, but I can show you how to implement it "by hand" (again, using Matematica). As I said, the ridge filter is the major eigenvalue of the Hessian matrix. I can calculate that eigenvalue symbolically in Mathematica:

$\text{eigenvalue}=\text{Last}\left[\text{Eigenvalues}\left[\left(
\begin{array}{cc}
 H_{\text{xx}} & H_{\text{xy}} \\
 H_{\text{xy}} & H_{\text{yy}}
\end{array}
\right)\right]\right]$

=> $\frac{1}{2} \left(H_{\text{xx}}+H_{\text{yy}}+\sqrt{H_{\text{xx}}^2+4 H_{\text{xy}}^2-2 H_{\text{xx}} H_{\text{yy}}+H_{\text{yy}}^2}\right)$

So what you have to do is calculate the second derivatives $H_{\text{xx}}$, $H_{\text{xy}}$, $H_{\text{yy}}$ (using a sobel or derivative of gaussian filter) and insert them into the expression above, and you've got your ridge filter.


What is the optimal way to select by color?

The Wikipedia page Color difference gives several formulas, the most advanced being CIEDE2000. They are all based on Euclidean distance $\Delta E$ (straight line distance) in the approximately perceptually uniform CIE Lab color space, with various corrections.


Suggestion for a DSP board to implement a Kalman filter in real-time at a rate of around 4MHz or higher

The conclusion of my discussion with hops and Marcus was that getting a DSP to operate with an ADC and DAC at the latency required for my application (<=250ns for 4MHz) would be difficult and require writing some custom interface. We concluded that the best option is to continue to work with the FPGA and improve my VHDL implementation using pipelining, bit widening to reduce my usage of large numbers of bits and turning my functions into architecture blocks. An alternative solution was to attempt to use Vivado HLS to generate VHDL code from C code.


It seems to me based on your question, that you are still missing some important concepts related to the algorithm, but you seem to think that the problem is not the concept. Let me see if I can clear this up.

First of all, you seem to be referencing the book Multirate Systems and Filter Banks by Vaidyanathan (pages 53-55). This would have been helpful information that you should have included in your post (or whatever other reference you obtained this from). The exposition in the book is strictly for designing eigenfilters with a low pass response (ideally with unity gain in the passband and zero gain in the stopband). The tutorial paper On The Eigenfilter Design Method and Its Applications: A Tutorial and its references extend this method to the design of filters with an arbitrary response. You should be following those references for your problem. The key point is to use reference frequencies in the frequency response to cast the cost function in the desired quadratic form so that the Rayleigh Principle can be applied to find a minimum solution (namely, the eigenvector corresponding to minimum eigenvalue).

I'd be happy to help walk you through my understanding of the filter design procedure for your particular case, but I have limited free time. For this reason, I won't do so unless I hear that you are interested.

Okay, since this is a Q&A website, and I want to keep this short, I will be glossing over a lot of details that you should definitely read on your own. The tutorial paper referenced above indicates that we can obtain an arbitrary response for the eigenfilter by minimizing the cost function
$$ \xi = \mathbf{h}^\dagger \mathbf{P} \mathbf{h} $$
where I have adhered to the notation in the paper so $\dagger$ indicates the conjugate transpose (a.k.a. hermitian). Ignoring the constants and the weighting function from the paper, the matrix $\mathbf{P}$ can be computed as
$$ \mathbf{P} = \int_{\mathcal{R}} \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^* \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^T d\omega $$
where $\mathcal{R}$ is the region over which $D(\omega)$ is defined, $D(\omega)$ is your desired frequency response, $\omega_0$ is a fixed reference frequency, and $\mathbf{e}(z) = \left[\begin{array}{cccc} 1 & z^{-1} & \cdots & z^{-N} \end{array} \right]^T$. When I have more time, I will provide more details on how to go about computing this (and converting it from being complex-valued to real-valued and linear phase). Although, this is computable using trapz as you suggested in your original question, and so it may be an acceptable solution already.

I believe that your specific case for $D(\omega)$ can be defined as follows
$$ D(\omega) = \left\{\begin{array}{cc} 1 & |\omega| \leq \frac{\pi}{3}\\
0 & \frac{5\pi}{12} \leq |\omega| \leq \frac{7\pi}{12} \mbox{ or } \frac{9\pi}{12} \leq |\omega| \leq \pi\\
\frac{12}{\pi}\left(|\omega| - \frac{15 \pi}{24}\right)  & \frac{15\pi}{24} \leq |\omega| \leq \frac{17 \pi}{24} \\
\mbox{undefined} & \mbox{otherwise} \end{array}\right.$$
where $\omega$ can vary between $-\pi$ and $\pi$. Let me know if that is correct. This gives you a total of 7 nonzero regions over which you need to integrate to obtain the elements of the matrix $\mathbf{P}$.

A convenient reference frequency for you will be $\omega_0 = 0$ or DC. This leads to a special form where $\mathbf{e}(e^{j \omega_0}) = \mathbf{1}$ and $\mathbf{1}$ is a column vector of all ones with the same length. Also, $D(\omega_0) = 1$ so that simplifies things a bit. For now, if you choose a spectrum with conjugate symmetry, you should obtain a real response. It would be best to re-derive the formula in terms of $\mathbf{c}(\omega)$ (a vector of cosines analogous to the paper's definition for $\mathbf{e}(\omega)$ that exploits the symmetry in the coefficients). You might attempt to use this information together with the equation for the matrix $\mathbf{P}$ to find your filter coefficients. If you do and post the MATLAB that you generate, then I can do more to assist you with this.


The task is to filter x(t) when given y(t), where y(t) = x(t) + n(t).
  Great but first we need to build an appropriate filter. At this point:

No. The task is to filter $y[n]$ to achieve $x[n]$. Your filter acts on the available signal $y[n]$ to get a best estimate of $x[n]$ from a noisy observation of it. BTW. discrete-time uses $y[n]$ or at least $y(n)$ notation, and not $y(t)$. In DSP, $t$ is reserved for continuous-time.

  1) Are we assuming that we have access to some training samples of
  x(t) and/or n(t) in addition to y(t)?

trainig samples is a term related with adaptivity, pattern recognition, neural networks etc. In the classical Wiener context, you do not need training samples but you need data statistics. See answer 2 for how to obtain the statistics.

  2) Or are we just assuming that we have access to statistical
  properties(auto correlation, cross correlation) of x(t) and/or n(t) in
  addition to y(t)? If so since x(t),n(t) are not explicitly available,
  are these statistics computed from other signals available to us with
  similar statistical properties?

Yes. In the Wiener filter context, you either assume (for theoretical development) or estimate (for practical applcations) the correlations (or Power Spectral Densities) between all those signals. In practice you estimate them from avaiable data by several means, depending on the application. You need data for this; and that's probably what you call training samples.

  3) Or is it a combination of both, such as having access to samples of
  x(t), but having only access to second order statistics of n(t) (or
  vice versa)?

You do not have access to relevant noise-free $x[n]$ by any means; otherwise why filter $y[n]$ to get $x[n]$? But, in some cases, one can have noise free irrelevant blocks of $x[n]$ samples, in addition to the relevant noisy blocks. So you can use those noise free blocks to estimate certain parameters of the process $x[n]$, so as to filter out the noisy blocks.

  4) Is the noise n(t) always assumed to be uncorrelated to x(t) in
  Wiener filtering, or are the textbooks going for the simplified
  version? What is the default scenario?

Yes. The noise $n(t)$ (or better $v[n]$ in discrete-time context) is almost always modeled as additive, white, zero-mean Gaussian noise uncorrelated with $x[n]$. This simplifies all the mathematical development and in typical cases is a good model. But you can relax some of those assumption, if that would yield a better (but more complex to derive) filter...

  5) Why are all the examples coming up with all-zero filters but not
  using pole-zero filters? Is it just a matter of simplified textbook
  treatment or is the problem definition not amenable to using transfer
  function modelling approach? Yes the first approach is computationally
  easier, but these days we have significant computing power and I don't
  think it is a good justification anymore, if that's the only reason.

IIR filters are good for their computational efficiency (and in Wiener context, they also yield the minimum MSE compared to FIR) but they may suffer from unstability and causality problems, and their design and analysis is more complex in adaptive filtering context. For this reasons, it's simpler and safer to use FIR based filters for adaptive Wiener applications.

To learn more about the Wiener filters and adaptive Wiener filtering:

Simon Haykin - Adaptive Filter Theory.
Monson Hayes - Statistical Digital Signal Processing and Modeling.


What are some practical examples of complex signals?

Here is one aspect that often creates confusion: In many cases, complex signals are used solely as a mathematical convenience, not because the complex representation is required by anything physical.
A good example is circuit analysis in electronics: all physical quantities (voltages, currents, impedances) are real and you can do the analysis with these quantities and the associated physical equations directly.
However, it's a lot easier if you convert the real quantities to complex ones, run the analysis in the complex domain and convert the results back to real. Hence, that's the way everyone does it.
So you don't HAVE to use complex signals, it just gets a lot easier if you do.


How would I calculate the DC Gain for a Chebyshev Type 2 filter?

I suppose that you obtain a transfer function of the desired discrete-time system in the form

$$H(z)=\frac{b_0+b_1z^{-1}+\ldots +b_Nz^{-N}}{1+a_1z^{-1}+\ldots +a_Nz^{-N}}\tag{1}$$

From that transfer function you compute the coefficients of the second-order sections. Before doing that, you can make sure that the DC gain of $(1)$ equals $1$ (which of course only makes sense for a low pass or a band stop filter).

The DC gain of $(1)$ simply equals

$$H(1)=\frac{b_0+b_1+\ldots +b_N}{1+a_1+\dots +a_N}$$

Obviously, the DC gain is independent of the chosen filter type.


Instead of capturing with shutter on, I suggest you to capture a picture in low-light and set your exposure time as minimum as possible and the reduce the aperture. Capture a neutral surface(where R,G,B are equal,preferably neutral gray) so when averaging the entire frame you will know noise level change by particular color.Do your calculations on Bayer RAW instead of jpeg because RAW is directly from sensor without any extra processing.

  Do you think lens aperture will play a significant role?

Yes.Lens aperture will play a significant role. If you are reducing the aperture it means the amount of light reaching the sensor is less and it has to be compensated by either increasing exposure time or the gain(Increasing the gain increases the noise).


DCT and DFT can be thought of as a bank of bandpass filters.  So it's sort of a filter, but it doesn't remove any part of the signal.  The whole signal is still there, just in a different format.

Yes, you could do the DCT on the signal and then high-pass filter by multiplying by the DCT of a high-pass filter's impulse response.  I'd suggest using DFT/FFT first though; it's a little more intuitive, though for a final application DCT has better energy-compaction, and so can compress signals better.


which are saturation and intensity images in this figure?

When the text refers to Hue component image, it is to represent what hue value these colors hold, not an actual picture in the text. When the text says H is shown in figure b, it means the numerical representation of these three colors (or hue) are illustrated more clearly when we take away saturation.

That being said, the saturation for all four figures (a to d) is being held constant as well. Figures a and d having high saturation, while figures b and c have no saturation. So I wouldn't say any of these figures are a Saturation image, by the definition you might be looking for. To make that picture you would need a square of a single color (Hue) and then change the saturation of each of the four corner squares.

Likewise, all of these figures have constant intensity in each corner. So I wouldn't say any of these are Intensity images either.

All of these figures have "component images" in them: H, S, and I.

Figure a and b illustrate the effects of changing Hue on an image with static saturation and intensity, while figure c and d show what happens when you "blur" (taking color averages over small portions of the image) the first two figures.


No you don't need to zero-pad each chirp. MATLAB's cspd function divide the signal into segments of length nfft. Would you want is to zero-pad each segment. Luckily for you, the cspd function can do it for you. Take a look here:

  If $\tt nfft$ is greater than the signal length, the data is zero-padded. If $\tt nfft$ is less than the signal length, the segment is wrapped so that the length is equal to $\tt nfft$.

So what you should do is increase the value of nfft from $512$ to $2048$.


What's the difference between the Gabor and Morlet wavelets?

The Gabor wavelet is basically the same thing.  It's apparently another name for the Modified Morlet wavelet. Quoting from Wavelets and Signal Processing:

  [The Modified Morlet wavelet] does not satisfy the admissibility
  condition but is nonetheless commonly used.  Sometimes this wavelet is
  called the "Gabor wavelet," but that term is improper because Gabor
  had nothing to do with wavelets.  He was ... one of the founders of
  time-frequency analysis.

That book is a collection of papers, and that paper ("The Wavelet Transform and Time-Frequency Analysis") is by Leon Cohen (of time-frequency distribution "Cohen class" fame), so I think it's reasonably authoritative.

At the very least, it sounds like the confusion is just a naming disagreement.  According to A Friendly Guide to Wavelets (pg. 114), Gabor was the first person to propose the use of Gaussian windows for time-frequency localizations, so his name tends to get attached anytime they're involved.


Frequency response creates non-causal impulse response

I am familiar with similar phenomena. During measurements with usage of sweep-sine, when you convolve your signal with inverse filter you are obtaining your IR. Although if you consider perfect case of convolving your original sweep $x(t)$ with the inverse filter $f(t)$, then you should get perfect Kronecker delta. Unfortunately what you will get is an impulse with pre and post-ringing:

One way to minimise this effect is to use so called Kirkeby Algorithm, described here by A.Farina (my guru of acoustic measurements): Advancements in Impulse Response Measurements by Sine Sweeps. I am sure will find it quite fascinating.

One last tittle-tattle:

  Commonsense suggests to me that it should have a very sharp attack as the main wave hits followed by diminishing ripples corresponding to echoes.

You are definitely right, although in real world measurements it is possible to have preceding spike before main impulse in some cases. For example if you place your source close to the wall, then wave from speaker will excite mechanical wave in a wall and wave will propagate in the material with speed of sound much higher than in case of air. Nextly wall will emit acoustic wave that will be captured by your microphone.


Does the Shannon theorem not apply when the amplitude of a wave is changed faster than half the time period of the wave?

Once you start changing the amplitude you are increasing the bandwidth of the signal. That's called "amplitude modulation" and the highest frequency is now the sum of the original frequency and the highest frequency in the modulation signal.
The sampling theorem still holds. You still only need twice the bandwidth but the bandwidth has increased due to the modulation


How to minimize peak-to-peak amplitude of periodic waveform through FFT phase shifts?

There are algorithms for determining low-peak-factor signals, for example:

  "Synthesis of low-peak-factor signals and binary sequences with low
  autocorrelation (Corresp.)", M. Schroeder, 1970

I have used this for flat-spectra signals to minimise peak to peak amplitude.

However I think this may be an XY problem. Do you have a limit on how much distortion you can introduce to the signals? Have you tried using a dynamic range compression algorithm? There are limiting algorithms which can increase the RMS amplitude significantly without introducing significant audio artifacts.


z-Transform Methods: Definition vs. Integration Rule

Just to avoid a misunderstanding: the $\mathcal{Z}$-transform is a transform defined for sequences, comparable to the Laplace transform for continuous functions. What you are talking about is not the $\mathcal{Z}$-transform, but methods for converting analog to digital (actually, discrete-time) systems. [And it doesn't help that one of those conversion methods is called the matched Z-transform, which is just an unfortunately chosen name.]

There are several methods to achieve such a conversion, and all of them have certain properties. There is no ideal or "correct" mapping, and some methods may be more well-suited for certain applications than others.

The methods mentioned in your question are direct mapping methods, where the variable $s$ in the continuous-time transfer function is simply replaced by a (rational) function of $z$. Apart from the two methods mentioned in your question (1. Tustin's method - also called bilinear transform, and 2. forward Euler), there is also the backward Euler method, which replaces $s$ by $(1-z^{-1})/T$. Note that Tustin's method (bilinear transform) is equivalent to the trapezoidal integration rule. The important advantage of the bilinear transform for the conversion of analog systems is that the frequency axis ($j\omega$-axis) of the $s$-plane is transformed to the unit circle of the $z$-plane, and, consequently, the frequency response of the resulting discrete-time system is just a compressed (warped) version of the analog frequency response. One consequence of this is that analog frequency selective filters with optimal magnitudes (Chebysev, Butterworth, elliptic) are mapped to discrete-time filters satisfying the same optimality criteria. Furthermore, since the left half-plane is mapped to the inside of the unit circle in the $z$-plane, stability (and the minimum-phase property) are preserved by the transformation.

This is different for the forward and backward Euler methods. With those methods, the $j\omega$-axis is not transformed to the unit circle of the $z$-plane, and the left half-plane is not mapped to the region $|z|<1$. In case of the forward Euler method, the $j\omega$-axis is mapped to the line $\textrm{Re}\{z\}=1$, and the left half-plane is mapped to the half-plane $\textrm{Re}\{z\}<1$. This means that a stable analog prototype filter might be transformed to an unstable discrete-time filter. The backward Euler method transforms the $j\omega$-axis to the circle $|z-\frac12|=\frac12$, and the left half-plane is mapped to the region inside that circle: $|z-\frac12|<\frac12$. Consequently, stability and the minimum phase property are preserved, but the analog frequency axis is not mapped to the unit circle of the complex plane, so the resulting frequency response is not just a warped version of the original analog frequency response.

Note that there are also other very common methods for converting analog to digital systems, which are not based on direct mappings from the $s$-plane to the $z$-plane. These methods are based on time domain criteria. The two most common of these methods are the impulse-invariant method and the step-invariant method. They preserve the shape of the impulse (or step) response of the analog system by defining the desired digital system as the one whose impulse (step) response is a sampled version of the corresponding analog impulse (step) response. These methods introduce aliasing in the frequency domain, but - unlike the bilinear transform - they do not introduce frequency warping.


Amplitude in frequency domain of the same signal measured at different bandwidth (f_max - f_min) and different FFT lines

Are you using the correct FFT form for spectral analysis?  You need to divide by N when mapping time->frequency domain with an N-point FFT for a power spectrum. If you want power spectral density you have to factor in the bin bandwidth and window effective noise bandwidths.
The amplitude of a spectral line represents the power of a single signal frequency, power spectral density reflects power per unit bandwidth, which is different.


SoapBoxLabs have very good technology for children, including pronunciation scoring


CDMA is not a protocol. It's a multi-access scheme.
A BER can't be generally given for a multi-access scheme. You'll need to define your access characteristics, SNR scenario, look at the specific multi-user detector for your specific modulation.

Hence, BlackMath is right: you'll need to go down the long route: Learn about CDMA's basics; do a lot of estimation theory, and apply it to the problem that you're looking at.

For many types of CDMA, detection probabilities are relatively easy to derive; things get a bit complicated because you'll need to define how the different users access the channel, in terms of occupation probability, and synchronity, as well as probability distribution of receiver powers.

Then you'd need to model how your multi-user detector deals with that. As you can see, this requires stochastic modelling of your scenario.


This is how you can try to design a short and simple FIR filter by hand. Note that this method is just supposed to be enlightening, a really useful filter can be designed using some software. If you want a zero at 60 Hz you indeed need to place it at an angle of $0.24\pi$ on the unit circle of the $z$-plane. The general relation between angle and frequency is

$$\phi=2\pi\frac{f}{f_s}$$

where $f$ is the desired frequency, and $f_s$ is the sampling frequency. Note that if the filter coefficients are real-valued, zeros (and poles) always occur in complex conjugate pairs. So if you have a zero at $z=z_0$ you must also have a zero at $z=z_0^*$. This just reflects the fact that the spectrum of a real-valued filter is symmetric: $H(\omega)=H^*(-\omega)$. Since you want to design a causal FIR filter, all poles are at the origin of the $z$-plane, and you only have control over the zeros.

So let's design an FIR filter of order $3$ (i.e. length $4$) with a zero at $z_0=e^{j0.24\pi}$, and - because it's a high pass filter - a zero at $z=0$. The zero at $z=0$ contributes a factor

$$H_1(z)=1-z^{-1}$$

to the total transfer function, and the zero at $z_0=e^{j0.24\pi}$ (and its complex conjugate counterpart at $z_0^*=e^{-j0.24\pi}$) contributes a factor

$$H_2(z)=(1-z_0z^{-1})(1-z_0^*z^{-1})=1-(z_0+z_0^*)z^{-1}+|z_0|^2z^{-2}=\\
=1-2\Re\{z_0\}z^{-1}+z^{-2}=1-2\cos(0.24\pi)z^{-1}+z^{-2}$$

because $z_0+z_0^*=2\Re\{z_0\}=2\cos(0.24\pi)$, and $|z_0|=1$ (the zero is on the unit circle). The total transfer function is

$$H(z)=H_1(z)H_2(z)=(1-z^{-1})(1-2\cos(0.24\pi)z^{-1}+z^{-2})=\\
=1-[1+\cos(0.24\pi)]z^{-1}+[1+\cos(0.24\pi)]z^{-2}-z^{-3}\tag{1}$$

From (1) you can read the filter coefficients:

$$h_0=1,\quad h_1=-[1+\cos(0.24\pi)]=-2.4579,\quad h_2=-h_1,\quad h_3=-h_0\tag{2}$$

Note that we just imposed the locations of the zeros but we didn't apply any normalization. A useful normalization for a highpass filter is to require $H(-1)=1$, i.e. the response at the Nyquist frequency (half the sampling frequency) is $1$. Since

$$H(-1)=\sum_n(-1)^nh[n]=h[0]-h[1]+h[2]-h[3]=2h[0]-2h[1]$$

we can define new normalized filter coefficients by

$$\hat{h}[n]=\frac{h[n]}{2(h[0]-h[1])},\quad n=0,1,2,3$$

The magnitude of the resulting frequency response of this filter looks like this (the phase is linear due to the symmetry of the coefficients):

You see that the three constraints (zeros at DC and at 60 Hz, and unity response at Nyquist) are satisfied. The transition from stopband to passband is of course not very steep but that's as good as it gets with an FIR filter of order 3.


yes - a video signal would be considered 3rd order tensor.

scalars are 0-order tensors, vectors are 1st order tensors, matrices are 2nd order tensor,  data volumes are 3rd order tensors and so on for multi-dimensional arrays.


determine two signals with a scale factor

A necessary (but not sufficient) conditions for $f_2$ to be a temporally scaled version of $f_1$ is that a spectral representation with a logarithmic frequency scale (such as the constant-Q transform) of $f_1$ is a translation of a log-frequency spectral representation of $f_2$.

Practically, given two signals, you can perform the test and evaluate $a$ by computing the CQT of $f_1$ and $f_2$, cross-correlating them and looking at the location of the peak. The strength of the peak might give you an idea of the spectral similarity of the two signals irrespectively of their temporal scale; and the position of the peak will give you the temporal scaling factor.

This type of representation robust to temporal scaling is useful in music signal modeling, where the different notes produced by a music instrument are - in a very rough approximation - temporally scaled versions of themselves.


Identify fundamental frequencies in short-period waveform

Finding the sum (and/or difference) of source frequencies in an FFT result usually means that there is an inter-modulation distortion in the source mixing (mixer).  Likely a multiplication instead of a strictly linear summing of source signals.

Regarding your last bit of thinking, trying to do blind source separation using just a short FFT is trying to do something unusual.  Any DSP solution is extremely unlikely to be simpler.


Say a discrete signal $x(n)$ and its Z-transform $X(z)$, discrete time Fourier transform $X(e^{j\omega})$, and discrete Fourier transform $X(k)$. The relationships between DFT and ZT / DTFT are respectively:
$$
X(k) = X(z)|_{z=e^{j2\pi k/N}} \tag{1.1}
$$
$$
X(k) = X(e^{j\omega})|_{\omega=2\pi k/N} \tag{1.2}
$$
Since the original signal $x(n)$ can be perfectly reconstructed by IDFT, yes you can certainly obtain DTFT from DFT. However, you can't get DTFT by padding zeros as long as possible but by a technique called frequency domain interpolation. First is reconstruction of ZT from DFT:
$$
X(z) = \sum_{n=0}^{N-1}x(n)z^{-n}=\sum_{n=0}^{N-1}\Big[\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_N^{-nk} \Big]z^{-n} \\ =  \frac{1}{N}\sum_{k=0}^{N-1}X(k)\Big[ \sum_{n=0}^{N-1}W_N^{-nk}z^{-n} \Big] \\ = \frac{1}{N}\sum_{k=0}^{N-1}X(k)\frac{1-W_N^{-Nk}z^{-N}}{1-W_N^{-k}z^{-1}} = \frac{1-z^{-N}}{N} \sum_{k=0}^{N-1}\frac{X(k)}{1-W_N^{-k}z^{-1}} \tag{2}
$$
where $W_N^{nk} = e^{j\frac{2\pi}{N}nk}$. Define the interpolation function
$$\varPhi_k(z) = \frac{1}{N}\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}} = \frac{z^N-1}{Nz^{N-1}(z-W_N^{-k})} \tag{3}$$
and $X(z)$ can be written as
$$
X(z) = \sum_{k=0}^{N-1}X(k)\varPhi_k(z) \tag{4}
$$
Now let $z = e^{j\omega}$ we can derive that
$$
X(e^{j\omega}) = \sum_{k=0}^{N-1}X(k)\varPhi_k(e^{j\omega}) \tag{5}
$$
and the interpolation function
$$
\varPhi_k(e^{j\omega}) = \frac{1-e^{-j\omega N}}{N(1-W_N^{-k}e^{-j\omega})} = \frac{1-e^{-j\omega N}}{N(1-e^{-j(\omega-2\pi k/N)})}\\ = \frac{1}{N}\frac{\sin(N\omega/2)}{\sin\big[(\omega-2\pi k/N)/2\big]} e^{-j[\omega(N-1)/2 + k\pi/N]} \\
=\frac{1}{N}\frac{\sin(N\omega/2)}{\sin\big[(\omega-2\pi k/N)/2\big]} e^{jk\pi(N-1)/N} e^{-j(N-1)\omega/2} \tag{6}
$$
Let's write it in a more convinient form as
$$
\varPhi_k(e^{j\omega}) = \varPhi(\omega - 2\pi k/N) \tag{7}
$$
where $\varPhi(\omega)$ is the Fourier transform of a $N$-point rectangle window divided by $N$
$$
\varPhi(\omega) = \frac{1}{N} \frac{\sin(\omega N/2)}{\sin(\omega/2)} e^{-j(N-1)\omega/2} \tag{8}
$$
Finally, we get the equation of reconstructing DTFT by interpolating DFT
$$
X(e^{j\omega}) = \sum_{k=0}^{N-1}X(k)\varPhi(\omega-2\pi k/N) \tag{9}
$$


IIR Cascaded Filter Response, Unexpected Results

Looks like a simple sign error in the denominator. If you use a = [1 -1.866685 .91058], you get something that looks notch-like at around 380 Hz.

  The Form is Direct Form II

This is numerically the worst. Shouldn't matter here, but it's good practice to either stick with Direct Form I or Transposed Form II.


Pitch is psychoacoustic.  Thus, a harmonic frequency can be stronger than the fundamental, without changing the pitch heard.  This is common in some stringed instruments.  And a note's pitch frequency could even be missing from a notes frequency spectrum.  For instance: An old telephone circuit can cut off frequencies below 200 Hz, but people still have no problem recognizing a low male voice pitched at around 100 Hz.

For many timbres, humans will hear the harmonic spacing as the note pitch.  In your plot, the harmonic spacing appears to be around 130 Hz, so that will be the note pitch perceived (if played at the same sample rate at the plot axis suggest).  Not 260.


Why dft of moving average impulse response has 0 angle

You are looking at the frequency response at multiples of $2 \pi /5$, where the response is either 1 or 0.

If you add several 0 values to the end of irs, you'll get the frequency response at more intermediate frequencies, where you can see that phase is not zero.


Generating a signature to identify RF modulation types

Yes, this is quite a studied problem in certain fields. I co-own a patent on it and we used one out of many approaches. We broke the problem into a feature extraction step and then a classification step. Some features can be the bandwidth, the continuity of the phase, the amount of phase clusters for digital modulation, etc. Then we used a Mahalanobis distance comparison to known modulations. I've also seen linear discriminants and neural networks used to process the features.

Certain modulations also exhibit certain shapes once you square them, cube them, etc. and then compute the power spectral density. Often if you do an FM or AM demodulation you can see spectral lines for digital modulations. Cyclostationary computations can also be helpful in digital modulation detection. The post above is right in regards to the search topic: "modulation classification"


You did not specify what sampling frequency you used to compute Fourier Transform. Software like MATLAB/GNU octave are only simulating an analog signal with high enough sampling rate. As a very general rule of thumb, if you know that time domain signal has a high periodic component of 1KHz, a 10x sampling frequency will give comprehensive results.
The pseudo-random waveform you showed on the left is a real and periodic signal and so we see real part of the FT is even-symmetric and imaginary part is odd-symmetric essentially meaning you can drop the left half plane from FT and focus on frequency values greater than zero. If you choose a sampling frequency correctly, magnitude response of the FT should show you peaks on frequencies which have the major share in your time domain signal e.g. 1KHz and 3KHz cosine signals may makeup (of course along with many smaller contributions from other cosine signals) the waveform. The real and imaginary parts are useful to compute the phase plot from FT essentially telling you how much delay should be added to cosine signals in order to make an approximation of your time-domain waveform.


How the Gain Term $ K \left( n \right) $ Is Derived? Why Is It Called Gain?

In all adaptive signal processing schemes, be it a Least Mean Squares (LMS), Recursive Least Squares (RLS) or a Kalman Filter, The fundamental concept is the update of some parameter: such as the vector filter coefficients of a Transversal Tapped Delay Line (an FIR indeed) structure of LMS or a state of a dynamical system as in the Kalman filter.

The generic structure of the update mechanisms resembles the following:

 parameter(n+1) = parameter(n) + K(n)*errorfunction(n)

in which the thing called parameter is updated from its previous value at index n to its current value at (n+1), based on a correction term that is some meaningful function of error (denoted by errorfunction(n)) and weighted by a specific multiplier K(n) commonly referred to as the GAIN.

As you can see, the gain multiplier, determines the amount of correction that will be applied to the current value of the parameter so as to update it to its next value. That's why its called as the gain. The larger the gain K(n) the more correction there will be hence more responsive will the filter be untill it becomes unstable due to too much gain etc... The gain K(n) shall be made variable based on the convergence properties of the filter.

The specific form and derivation of an optimal gain term is completely dependent on the filter structure and the necessary mathematical manipulations of the signals inside the processing blocks, so as to get the required update equations...


Where did we get the DC term of the Accumulator from DTFT?

Note that the equation
$$x[n]=y[n]-y[n-1]\tag{1}$$
doesn't uniquely determine the sequence $y[n]$. If some $y[n]$ satisfies $(1)$, so does $y[n]+c$ with some constant $c$. Consequently, Eq. $(1)$ determines $y[n]$ only up to a constant, which corresponds to a DC term in the frequency domain.
Hence, the equation
$$Y(e^{j\omega})=\frac{X(e^{j\omega})}{1-e^{-j\omega}}\tag{2}$$
is only valid for $\omega\neq 0$, or, more generally, for $\omega\neq 2\pi k$, $k\in\mathbb{Z}$.
The most straightforward way to derive the DTFT of
$$y[n]=\sum_{m=-\infty}^nx[m]\tag{3}$$
is to realize that $y[n]$ is the convolution of $x[n]$ with the unit step sequence $u[n]$, and, consequently, the DTFT of $y[n]$ is given by
$$Y(e^{j\omega})=X(e^{j\omega})U(e^{j\omega})\tag{4}$$
where $U(e^{j\omega})$ is the DTFT of the unit step $u[n]$:
$$U(e^{j\omega})=\pi\delta(\omega)+\frac{1}{1-e^{-j\omega}}\tag{5}$$
Several ways to derive the result $(5)$ are discussed in the answers to this and this question.
Note that in $(5)$ I use $\delta(\omega)$ instead of $\sum_k\delta(\omega-2\pi k)$ to avoid cluttered notation. The fact that the DTFT is always $2\pi$-periodic is understood.


This depends a bit on what you want to get out of this and how much effort/work  you are willing to put in.

Doing a room EQ that actually works and makes it sound consistently better is quite complicated. There are commercial systems available but they tend to be complicated and expensive or tied to a specific product.

A really good freeware option to play around with is Room EQ Wizard https://www.roomeqwizard.com/

If you want to do this yourself, I suggest you break it down in a few separate steps that can be tackled individually

Figure out how to make a good measurement in a room (noise, reverb, microphone, ADC, SNR estimation, etc.)
Figure out what to measure (using the tools from step 1) that's representative enough of the entire room/speaker/listener situation
Extract perceptually relevant parameters from these measurements
Turn these parameters into some sort of correction that would move the parameters towards some perceptual target
Design a filter and/or algorithm based on this correction
Design the process/HW/system to actually implement and apply the filter or correction

All of these steps have been studied in detail and there is a fair bit of literature available on any of them.


An FFT approach to filtering/resampling should work fine, but you don't want to just truncate or pad with zeros.  That will add effects to the image.  See here for more information about FFT resampling.

You could also just do a FIR filter.  If you upsample by a ratio of $\frac{p}{q}$ then you would want your anti-aliasing filter to have a normalized bandwidth of $\frac{q}{p}$.  You could make the bandwidth a little smaller to get rid of more alias energy, as desired.  If you downsample by a ratio of $\frac{p}{q}$ then you would want the filter bandwidth to also be about $\frac{p}{q}$.  Again, you could make it a little smaller to get rid of more alias energy.

Once you have your filter you just need to do a two-dimensional convolution and trim the extra resulting samples on each edge.  Many signal processing libraries have built-in two-dimensional filtering routines that can do this for you.


Alamouti code after or before IFFT block in OFDM

The basic concept of OFDM is to divide a high-bitrate datastream into $N$ low-bitrate datastreams and to multiplex these low-bitrate datastreams in frequency. That is, every datastream is assigned to a distinct frequency band (so-called subcarrier) that does not interfere with the other frequency bands. Orthogonality allows the frequency bands to be packed close together which increases the spectral efficiency of OFDM. As it turns out the frequency multiplex can be implemented very efficiently using the IFFT algorithm that transforms a frequency-domain signal into a time-domain signal. Every $N$ complex QAM symbols are one input vector to the IFFT which produces a vector of $N$ time-domain samples that are then transmitted sequentially over the channel. The construction of the transmit signal in frequency domain has several advantages:

Fast fading channels (often the case for multipath mobile channels) might significantly attenuate some frequencies. OFDM enables adaption to such a channel by switching off some subcarriers or decreasing the modulation order on the respective subcarriers.
Usage of cyclic prefix as guard interval makes channel equalization easy to implement.
Dynamic resource allocation is possible (see LTE uplink)

Your suggestion to take $N$ QAM symbols, apply the FFT and then transmit the FFT output corresponds to transmitting the spectral analysis of some time-domain signal one bin per symbol interval. This is totally different from applying the IFFT and I don't see any advantage in doing so.

  What is the reason for alamouti code to be implemented before IFFT? Can it also be implemented the way I thought after cyclic prefixer block?

Let's assume we apply the Alamouti space time block code (STBC) after the cyclic prefix. Let $x_k$ be the discrete-time output of the cylcic prefix block. Then the coded signal vector is given by
$$
\mathbf{x} =
\left(
\begin{matrix}
  x_{k+1} & -x_k^*\\
  x_k     & x_{k+1}^*
\end{matrix}
\right)
$$
If we assume a perfect channel with channel matrix $\mathbf H = (1\quad1)$ and if we ignore noise, the received signal is
$$
\mathbf r = \mathbf H \mathbf x =
\left(
\begin{matrix}
  x_{k+1} + x_k & -x_k^* + x_{k+1}^*
\end{matrix}
\right)
$$
It is clear that we first have to decode $\mathbf r$ before we can demodulate it with the FFT. For doing so, the channel matrix has to be known. In general, the channel matrix is unknown at the receiver and thus has to be estimated. This channel estimation is done in frequency domain in OFDM systems, i. e. after the FFT block. As a consequence the STBC decoder has to be placed after the FFT as can be seen from the following figure:

(Figure is taken from Kim, J et al.: Receiver Designs for Alamouti Coded OFDM
Systems in Fast Fading Channels)

This in turn means that the STBC encoder has to be placed before the IFFT.


The first thing you should do is examine your signal collection hardware and do everything humanly possible to improve the signal-to-noise ratio of the analog signal applied to your A/D converter. I realize this is easier said than done, but do your very best to reduce the noise in your analog signal.
Next, you might consider the following process:

Based on the following figure, examine your noisy signal samples and identify the $T_n$ time instants of your QRS complex spikes.
Create multiple signal segments, $S_n$, where each segment comprises $N$ samples just prior to and including the QRS complex spike sample.
Average the $S_n$ signal segments.
Examine the "averaged" $N$-sample segment to see if your can detect an averaged p-wave.


Why should one choose the maximum as reference for power-to-decibel conversion?

dB is a power ratio, so when we see units in dB we are seeing units on a relative scale. The reason that a 0 dB reference is so common is because this is simply normalizing the number scale to 1. $10Log(1) = 0$ dB

This is similar to using 100% (1) instead of absolute numbers.


Is it not possible to perform non-linear filtering in frequency domain?

Linear filters in the spatial domain have a direct equivalent in the frequency domain, so you can transform your data and filter between spatial<->frequency domains and get equivalent results. However the transforms that we use for converting between spatial<->frequency domains are only valid for linear systems - a non-linear filter such as a median filter therefore has no frequency domain equivalent, so we can't perform e.g. median filtering in the frequency domain.

Note that this does not mean that there are no non-linear operations that we can perform in the frequency domain, e.g. we could reduce all frequency components below a certain threshold value to zero, it's just that there are no counterparts between the two domains for non-linear operations, so the frequency domain threshold operation described would have no equivalent operation in the spatial domain.


What are good estimates for phonemes length?

NLTK has an interface/downloader for a sample of the TIMIT phoneme corpus (http://nltk.sourceforge.net/doc/api/nltk.corpus.reader.timit-module.html), and it should be straightforward to get the lengths of the phoneme recordings in that corpus.


Evaluate phase of $X(\omega)$ without computing $X(\omega)$

hint: x(n) has symmetric coefficients which as a filter is a type 1 linear phase filter as the OP has stated. Such a filter will have a phase slope that is dependent on the number of taps only.


Why Wavelet developed when we already had Short-time Fourier transform

The short-time Fourier transform doesn’t offer better analysis of data than the discrete Fourier transform, it offers a different kind of analysis.  The DFT offers an exact decomposition of data to a frequency representation. The STFT offers an approximate decomposition to a time/frequency representation.  Which is better depends on what you are after.  The same holds true of the Wavelet transform.  Wavelet transforms can be thought of as decomposition to a time/frequency representation, but wavelet transforms generalize the concept of decomposition.  Different wavelet functions have been created so you can choose a decomposition that suits your needs.


FMCW Radar: Need for Second stage FFT for velocity determination?

Actually the second round of FFT is exactly doing that, as by definition differentiation of phase is frequency. In realistic case if target is moving, the target position will not change but its phase will change along the column (also called slow-time, where the row in your case is called fast-time), the second FFT shows this changing phase as Doppler frequency. Furthermore, what second FFT also does it something similar to pulse integration. More the number of rows larger is your integration gain.


Formula for 2D Image Filter for Scratch Detection

I agree with Peter K.  It looks like the authors of the Multidirectional Scratch Detection and Restoration in Digitized Old Images (E. Ardizzone, H. Dindo, and G. Mazzola, EURASIP Journal on Image and Video Processing, 2010) paper tried to reformulate the equation (shown by Peter K.) from the Analyzing Oriented Patterns paper by Kass and Witkin (1987; Computer Vision, Graphics and Image Processing, Vol. 37. No 3, pages 362 to 385).

Comparing the equations posted by the OP with the equations of a Gabor Filter, it feels like the author's are trying to make a novel contribution with their paper by reformulating the work of Kass and Witkin.  Although neither paper includes references to Gabor filters directly, Ardizzone et al. justify the use of their reformulated "steerable" filter by stating,

  The filter order n can control the slope of the subband, so that we chose n=4 to concentrate the filtering in a precise zone of the spectrum.  This explains the use of this bandpass filter instead of a Gabor one.

My experience with image filters is limited, but I believe a Gabor filter can realize the same apparent functionality as the equation given by the OP including bandwidth selection (see this entry, for example).

I think "subband" refers to one of the directionally dependent filters and the "center" is the origin of the Gabor kernel relative to the image on which the kernel is being applied.

I hope this helps.


What is the purpose of a window function?

I've read that a window function is a mathematical function that is zero-valued outside of some chosen interval.
Can someone please explain what this really means?

It means exactly what it says, but I assume you really want to know why do we care?
The purpose of window functions in signal processing is to remove (hopefully) extraneous information, so that subsequent analysis produces more meaningful results. Checkout wikipedia (and other sources) for specifics, but windows functions do pre-processioning in the time domain to help improve results of later processing in the frequency domain.
This hopefully helps explain why there are so many different window functions.  Basically, different window functions are designed to help with different problems present in the time domain data.  Often there is no right/wrong window functions, there are just more/less useful functions for solving the problem at hand.


Modelling the amplitude and phase through an amplifier

I apologize, the following is a bit rough, but I've put together an example in MATLAB illustrating all the different elements I think you need. That is:

a fist-order model for an amplifier (a common approximation for most op. amps.)
how to discretize it to find an IIR filter
one way to implement a realization of the IIR filter (straight forward state-space implementation)
how to generate a transfer-function estimate from measured data
one way (not the most robust, but perhaps the most obvious) to fit the transfer-function coefficients to the transfer-function estimate

To read up on the theory behind all this I would recommend the following text-books:

Chi-Tsong Chen, Linear System Theory and Design (3rd edition). Oxford University Press, 1999.

Petre Stoica and Randolph Moses, Spectral Analysis of Signals. Prentice Hall, 2004.

Lennart Ljung, System Identification: Theory for the User (2nd edition). Prentice Hall, 1999.

The last book is somewhat heavy, but concentrate on "output error" methods, and you should be very well covered for systems where you know the order of the model, or "subspace" methods if you don't. These methods are available in the MATLAB system identification toolbox.

%% Amplifier model
wb = 1e4; % Amplifier bandwidth
ga = 1; % Amplifier gain

% Modeling amplifier as having fist-order dynamics
% If using Octave run these in the command window:
%pkg install -forge control
%pkg load control
H = tf(ga,[1/wb 1]); % Continuous-time transfer function

%% Finding a realization (convert to state-space formulation)
% V/U = ga/(s/wb + 1)
% V*(s/wb + 1) = ga*U
% s*V = -wb*V + ga*wb*U
% dv/dt = -wb*v + ga*wb*u
%
% Standard form:
% dx/dt = A*x + B*u
% y = C*x + D*u
%
% with x = v, y = x
A = -wb; B = ga*wb; C = 1; D = 0;

%% Finding a discretization:
Ts = 1/1e5; % sampling time
% Using Euler integration
Ad1 = 1 + Ts*A; Bd1 = Ts*B; Cd1 = C; Dd1 = D;
% Using exact discretization (use when system has analog-digital converters)
Ad2 = exp(A*Ts); Bd2 = (Ad2 - 1)/A*B; Cd2 = C; Dd2 = D;

%% Computing a step-response to demonstrate one way to implement an IIR filter
u = [zeros(1,100) ones(1,200)]; % input
x = zeros(2,numel(u)+1);
y = zeros(2,numel(u));

t = 0:Ts:(numel(u)-1)*Ts;

% Comupting filter outputs
for k = 1:numel(u)
    x(1,k+1) = Ad1*x(1,k) + Bd1*u(k); % response of Euler integration scheme
    y(1,k) = Cd1*x(1,k) + Dd1*u(k);

    x(2,k+1) = Ad1*x(2,k) + Bd1*u(k); % response of exact discretization
    y(2,k) = Cd1*x(2,k) + Dd1*u(k);
end

figure(1)
plot(t,y(1,:),t,y(2,:))

%% Computing the response to noise for identification of parameters
% The frequency response of a system is in general most reliably found when
% exciting the system with white noise. A pseudo-random binary sequence is
% theoretically optimal when you don't know anything about your system, but
% can in practice cause problems due to the switched nature. A Gaussian
% noise source is usually a good substitute.
un = randn(1,1e4); % input
xn = zeros(1,numel(un)+1);
yn = zeros(1,numel(un));

% Comupting filter outputs
for k = 1:numel(un)
    xn(k+1) = Ad2*xn(k) + Bd2*un(k);
    yn(k) = Cd2*xn(k) + Dd2*un(k);
end
% Remove transient as you want to look at stationary response
un_ = un(201:end);
yn_ = yn(201:end);

% Compute power spectral densitites and cross-power spectral densitites
% using a simplified version of the Welch method
N = length(yn_); % length of original sequence
L = 5; % number of PSD averages
M = floor(N/L); % length of sequence segments
W = 0:(2*pi)/M:2*pi-(2*pi)/M; % normalized PSD frequencies

WIN = hanning(M); WIN = WIN(:).'; % optional windowing (if the segments are short)
Pwin = 1/M*sum(abs(WIN).^2); % window "power" correction (for Welch method)

% PSD for input sequence
Puu = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';
    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puu = Puu + (Uft.*conj(Uft));
end
Puu = 1/(L*Pwin)*Puu;

% Cross-PSD between input and output sequence
Puy = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';

    yn_segment = yn_((k-1)*M+1:k*M);
    yn_segment = yn_segment(:).';

    yn_windowed = yn_segment.*WIN;
    Yft = 1/sqrt(2*pi*M)*fft(yn_windowed);

    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puy = Puy + (Yft.*conj(Uft));
end
Puy = 1/(L*Pwin)*Puy;

% Transfer-function estimate
Huy = Puy./Puu;

figure(2)
subplot(2,1,1)
semilogx(W,abs(Huy))
subplot(2,1,2)
semilogx(W,unwrap(atan2(imag(Huy),real(Huy))))

%% Fitting a model to the output data

% finding the z-transform of discrete-time model (Euler integration)
% using the standard expression
z = tf('z',Ts);
Hd1 = Cd1*(z - Ad1)\Bd1 + Dd1;

% finding the z-transform of discrete-time model (exact discretization)
% using an alternative method
b0 = Bd2*Cd2;
a0 = -Ad2;
a1 = 1;
Hd2 = tf(b0,[a1 a0],Ts);

% compare the continuous and discrete-time models
figure(3)
bode(H,Hd1,Hd2)

zv = exp(1i*W); % starred transform (between z- and s-domains)

% H = b0/(a1*z + a0)
% H(W) = b0/(a1*exp(i*W) + a0)
% a1 = 1; % in order to have unique coefficients for a0 and b0
% H(W) = b0/(exp(i*W) + a0)
%
% (exp(i*W) + a0)*H(W) = b0
% exp(i*W)*H(W) = -a0*H(W) + b0
% exp(i*Ts*w)*H(W) = [-H(W); 1]*[a0, b0]
% zeta = phi*theta;

zeta = zv(:).*Huy(:);
phi = [-Huy(:), ones(size(Huy(:)))];

theta_true = [a0; b0]
theta_est = phi\zeta


Eureka! You can get a zero centered waveform of an audio clip by applying a window, i used hann (cosine) window and then inverting the left and right sides.


I assume you set the first momentum of your true measure to $\mu=0$, so your  question actually reads

  gyroscope with sensor drift does that mean that $E[f(x)] \ne \mu$?

No. It might still be the case that the sensor value, in the long term, drifts centred around the true value.

The thing is that $f(x)$ is a random variable with non-zero variance and a time dependency. That means: If there's drift, that random variable is not a stationary process.


Find the $\mathcal Z$-transform of this function?

Doing some calculations is appreciated. However, it is preferred that you write them in the question because a picture can vanish more easily. Do not hesitate to edit the math.

To help you with an hint in checking your answer, remember that:

$$ a^{-n} = \left(\frac{1}{a}\right)^n$$

a standard trick to play with known tables.


Filter pixels with respect to specific wavelength frequencies

Aliasing irreversibly transforms frequencies above the Nyquist limit to appear as frequencies below the limit. Almost no mathematical algorithms are effective compared with a physical anti aliasing filter in front of the sensor.

However, if the sampling frequency of the sensor itself is higher than the numerical aperture of the lens, it may not need the low pass anti aliasing filter any more since lens itself would serve as a low pass filter.

Regarding the color filter, I think it depends on your application. If you are try to image fluorescence, it is possible without a color filter if the emission spectrum is distinguishable with excitation spectrum. Background subtraction operations will recover the pixel information on the fluorescence. Even if several emission spectra are multiplexed, there is still possibility to differentiate them by solving a linear matrix equation under the assumption that the sensor system is a linear system with different wavelengths.

But color correction filters in digital photography may still needed to get more uniform exposure in the desired channels. Otherwise, the denoising process, though maybe possible through some image processing algorithms, is still challenging in complete noise removal without the prior knowledge of the noise distribution.


$ \frac{1}{f^2} $ Weighting for Least Squares FIR

Well the page you link to says:

  'differentiator' for type III and type IV filters, using a special weighting technique. For nonzero amplitude bands, the integrated squared error has a weight of (1/f)2 so that the error at low frequencies is much smaller than at high frequencies. For FIR differentiators, which have an amplitude characteristic proportional to frequency, the filters minimize the relative integrated squared error (the integral of the square of the ratio of the error to the desired amplitude).

Bolding mine.

So, following the paper:
$$
\epsilon_2 = \int_0^\pi W(\omega)(A(\omega)−D(\omega))^2d\omega
$$
for the differentiator, we want:
$$
W(\omega) = \frac{1}{\omega^2}\\
D(\omega) = \omega
$$
and $A(\omega)$ has the same form as before.

So
$$
\epsilon_2 = \int_0^\pi \frac{1}{\omega^2}(A(\omega) - \omega)^2 d\omega\\
= \int_0^\pi \left[\frac{A^2(\omega)}{\omega^2} - 2\frac{A(\omega)}{\omega} + 1 \right] d\omega
$$
Just substitute the expression for $A(\omega)$ and you're done.


that's pretty simple: in digital domain, real-world frequencies have no inherent "meaning".

That is, a period of a real-world sine of frequency 10 Hz might be 2, 200, 1233 or whatever samples long, depending on the sampling rate.

Thus, periods in DSP can only be measured in samples; for example, a sine of period 1 second has a period of 5.4 samples if the sampling rate was 5.4 S/s.

Now, logically, that also means that frequencies can only be related to the sampling rate. Hence, your $\frac{f_{LO}}{f_\text{sample}}$ (which is a unitless thing!) oscillator corresponds to an analog oscillator of frequency $f_{LO}$ (which has Hertz as a unit).

With your question in mind, it strikes me as odd that you're looking at a figure describing a receiver that mixes in the digital domain – in software defined radio, most receivers try to minimize the rate at which the ADC needs to run, simply for cost, efficiency, complexity, ADC dynamic range, and noise bandwidth reasons, and thus have the IQ mixer on the analog side – I'd call your design a bit of a "speciality". So, if this is you learning about direct conversion receivers, you might have picked an unlucky reference. If this is you learning about different types of digital receivers: That is extremely close to the theoretical perfect Software Radio, that does all the signal processing in digital domain.

Also note, in the figure, $f_s \overset !\ge 2\cdot f_{max}$ of the low-pass filter, not only one time that frequency (I assume "GHz" just means "order of magnitude is Gigahertzes"). That's why direct sampling is really an expensive technology for high-frequency signals (and why undersampling with filters or an IF makes a lot more sense, most of the time, since you don't care about everything happening between 0 Hz and the maximum frequency of your signal of interest, usually, but only about some limited bandwidth around a center frequency).