Smallest FFT buffer size given zero-padding

Zero padding does not increase frequency resolution, but interpolates the spectrum you have by adding more samples in between as an interpolation. The frequency resolution is well known to be set by the size of your buffer before zero padding with the relationship:

$$f_{res} = 1/T$$

Where T is the length of your buffer in time (notice that this relationship is independent of the sampling rate!). However this classic relationship is the width of each bin by their 3.0 dB bandwidths or half power point (viewed as a filter-- it is helpful and accurate to view the FFT as a bank of filters), while to completely resolve two frequencies of equal magnitude, due to their coherent addition in the FFT, the frequencies need to be separated by greater than the half magnitude width of each bin (or the 6.0 dB bandwidth).  The 6.0 dB bandwidth for a rectangular window (no window) is 1.21, so to resolve two tones of equal magnitude and any arbitrary phase the length N must be greater than:

Minimum number of samples, two tones same magnitude (no Window):

$$N > 1.21 f_s/f_{\Delta}$$

Where $N$ is the number of samples, $f_s$ is the sampling rate and $f{\Delta}$ is the frequency separation both in Hz, using the relationship:

$$T = \frac{N}{f_s}$$

Note that using any window will decrease this resolution further (the best resolution is obtained with the rectangular window, which is no window) at the needed benefit of dynamic range. Thus choosing a window is always a trade between frequency resolution and dynamic range, and therefore the window is chosen based on the anticipated maximum power difference between closest frequencies and number of samples needed.  For example, the 6.0 dB BW of the Hamming window is 1.81 which modifies the above formula accordingly:

Minimum number of samples, two tones same magnitude (Hamming Window):

$$N > 1.81 f_s/f_{\Delta}$$

Please refer to this paper by fred harris On the Use of Windows for Harmonic Analysis that provides detailed considerations for window selection and resolution bandwidth for resolving tones, particularly the section G on p 178 "Minimum Resolution Bandwidth" with reference to his figure copied below, along with Table 1 on page 176 that lists detailed metrics for all common windows.

In the case of multiple tones at different power levels and in the presence of background noise, it is helpful to look at the complete magnitude response of the kernel (frequency response of the "filter" for each bin). Note again, since the tones in different bins add coherently they will sum in magnitude based on their relative phases (and subtract!), while noise components will sum in power. This is shown below by comparing the normalized magnitude plots of the rectangular and Hamming kernel, with a horizontal line showing the -6 dB half magnitude point to relate to the relationships derived above. The width of the main lobe of the rectangular window is 1.21 and for the Hamming window it is 1.81, and with reference to the figure above we see how this is the minimum separation for two tones of equal magnitude, being the boundary condition where a single maxima can exist.  But we also see the details related to dynamic range considerations.

Specifically consider two tones separated by 20 dB, the rectangular window will require greater than 5.36 bins ($f_{\Delta} > 5.36 f_s/N$) while the Hamming window will be able to resolve with separations greater than 3.32 bins in this condition:

Minimum number of samples, two tones 20 dB difference (no Window):

$$N > 5.36 f_s/f_{\Delta}$$

Minimum number of samples, two tones 20 dB difference (Hamming Window):

$$N > 3.32 f_s/f_{\Delta}$$

The paper referenced above details this further with consideration to multiple tones and in the presence of noise, where due to the effect of the side-lobe roll-off the window chosen has a significant impact.

For consideration of noise, margin is added to your maximum separation in power between tones; theoretically for AWGN you can statistically always encounter a condition where the noise added to one tone in magnitude compared to the noise subtracted from the next is beyond your designed separation, therefore it becomes a probability of false detection problem, where you design a margin based on the probability of not separating two tones in your designed resolution. I suspect this would quickly be insignificant if you are working with high SNR's in a bin (and remember the FFT has a processing gain of 20LogN decreasing the noise in one bin compared to your overall noise). This is best demonstrated with an example: Consider the case of adjacent tones of equal magnitude with a noise limit we bound to 1.5 (which could also then reduce the adjacent tone by 0.5), this would result in a worst case separation in the tones of 20Log(1.5/5)= 9.5 dB.  If we therefore added 9.5 dB to our total maximum designed separation in magnitude between adjacent tones (and then use the process we used above to compute number of FFT samples needed) we will ensure that we will always be able to detect the adjacent frequencies for all cases that the noise itself does not cause a 3x magnitude separation between bins. (Placeholder here to show awesome graph of SNR vs false alarm rate in this case: This would be an interesting joint probability problem to solve at some point if I or someone else has the time and energy to work through it assuming not readily available somewhere).

So to answer your question, your smallest buffer size is set by the actual frequency resolution you need, with consideration to total number of tones in close proximity and dynamic range needed (max distance in magnitude between tones). Zero padding is done if you want to interpolate more samples (without adding any more resolution; it is just a smoothing), which you may do for convenience (to get the FFT out to the closest power of 2 for example, or to fill in more samples on a plot.

In the case of detecting the frequency of a single tone, the zero padded FFT will more readily point out the precise frequency location (but this is not different from what you can do more efficiently through interpolation using the information derived from the closest bins if you did not zero pad), but if you had multiple frequencies together within the span of the frequency resolution defined above, you would not be able to resolve them (they would appear as one frequency).

As @OlliNiemitalo has correctly pointed out that zero padding WILL allow you to visibly observe the obtained frequency resolution described above (which is what I have done to plot the responses in between bins in the plot above).

In this post is an example plot showing the difference between zero padding to increase the number of samples versus adding more actual samples:

What happens when N increases in N-point DFT


Pre-Warping is not working on bilinear transform

is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz?

Yes it is possible... .
I use MIM (Magnitude Invariance Method) for this. Here's plot showing 2nd order LPF responses for fc=20kHz-25kHz (1kHz step) @ 44.kHz sampling:

MIM method works well for even lower sample rate implementations but, there are also issues not sorted out yet.
Here's Matlab implementation of the MIM and PIM (Phase Invariance Method) method.


How do ADC/DACs deal with out of sync signals near Nyquist frequency?

How does it faithfully recreate the waveform if the samples are taken a points that are 0?

Not at all. Things at $\frac{f_\text{sample}}2$ might just as well be constants. You can't tell [0 0 0 0] (sampled 24 kHz signal at zero-crossings) from [0 0 0 0] (sampled zero constant signal).

Nor can you tell [0.1 -0.1 0.1 -0.1] (sampled 24 kHz signal with amplitude 1 just a bit after each zero-crossing) from [0.1 -0.1 0.1 -0.1] (sampled 24 kHz signal of amplitude 0.1, sampled at its extrema!).

You need the sampling rate to be higher (not only just as high) than two times the highest signal frequency component for this to be unambiguous.


Ok, get the solution, the sigmoid-funtion should be updated during the iterations:

x0 = 0
v = 0
# now starting with newton-method iterations
for i=1:10
    I_t = σ(x0 + v) - g(x0) # temporal derivative !!! sigmoid-function should be updatetd here
    I_x = dσ(x0)  # spatial derivative does not change
    v = v - (I_t/I_x)
    print("$(v) \n")
end


What is the algorithm to do a Discrete Hilbert Transform?

The Discrete Hilbert Transform is an "ideal" (implying that this is not the practical implementation) linear time-invariant filter with input $x[n]$ having output

$$ \mathscr{H}\Big\{ x[n] \Big\} \triangleq \hat{x}[n] = \sum\limits_{i=-\infty}^{\infty} h[i] \, x[n-i] $$

Where the impulse response of a discrete-time Hilbert transformer is:

$$ h[n] = \begin{cases}
 \frac{\big(1 - (-1)^n\big)}{\pi n} \quad & n \ne 0 \\
\\
0 & n = 0
\end{cases}$$

Because $1 - (-1)^n = 0$ for even $n$, this can be restated as

$$ h[n] = \begin{cases}
 \frac{2}{\pi n} \quad & n \text{ odd} \\
\\
0 &  n \text{ even}
\end{cases}$$

This is not a causal impulse response, nor is it finite in length.  To make it finite in length, you would need to window it with a decent window:

$$ h[n] = \begin{cases}
 \frac{2}{\pi n} w[n] \quad & n \text{ odd} \\
\\
0 &  n \text{ even}
\end{cases}$$

where $w[n]$ is some window function of width $L+1$ samples (and $L$ is an even positive integer).  If it were a Hamming Window, it would be:

$$ w[n] = \begin{cases}
  0.54 \ + \ 0.46 \cdot \cos\left(\pi \frac{n}{L/2} \right) \qquad & |n| \le L/2 \\
  0 & |n| > L/2 \\
\end{cases}$$

If it were a Kaiser window it would be

$$ w[n] = \begin{cases}
  \frac{1}{I_0(\beta)} \, I_0\left(\beta \sqrt{1 -  \left(\frac{n}{L/2}\right)^2 } \right) \qquad & |n| \le L/2 \\
  0 & |n| > L/2 \\
\end{cases}$$

where

$$ I_0(u) \triangleq \sum\limits_{k=0}^{\infty} \frac{(-1)^k \big( \tfrac{u}{2} \big)^{2k}}{(k!)^2} $$

$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $L+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $L$) and, in my centered and symmetrical case, must be even. $\beta$ is a "shape parameter", maybe around 5 or 6, i dunno.

Now, to make this causal, your impulse response has to be delayed by $\frac{L}2$ samples to be $h[n-\frac{L}2]$, but then should also all other signals that this Hilbert output is compared to, they should also be delayed by $\frac{L}2$ samples to keep the phase relationship correct.

With a windowed finite-length impulse response (which is what we call an "FIR"), the Hilbert transform output is

$$ \hat{x}[n] = \sum\limits_{i=-L/2}^{L/2} h[i] \, x[n-i] $$

and delaying the output so that the filter is causal you get

$$ \hat{x}[n-\tfrac{L}2] = \sum\limits_{i=0}^{L} h[i-\tfrac{L}2] \, x[n-i] $$

but with this delayed output $\hat{x}[n-\tfrac{L}2]$, you must compare that only to the like delayed input $x[n-\tfrac{L}2]$ in order for the two signals to have their 90° phase relationship (which is fundamentally what the Hilbert Transform is about).  Note that every even-indexed sample of $h[n]$ is zero, so the number of taps is not really $L+1$ but is $\frac{L}2 - 1$ taps with non-zero coefficients.


Given a Horizontal Derivative Kernel What Would Be the Output Image?

The solution provided is correct. Here is why.

To compute the derivative at index(3,5)

d(3,5) = value_at(4,5) - value_at(2,5)
d(3,5) = 150 - 0
d(3,5) = 150

Hence the solution provided there is correct. Hope this helps :)


Design Filter with Arbitrary Magnitude and Phase Response in Matlab

for FIR, the simplest way is to inverse FFT you two-sided complex frequency response:

$$ h[n] = \operatorname{iFFT} \left\{ H[k] \right\} $$

(use a Kaiser window to window the FIR to a length that is acceptable to implement.)

i presume you have access to firpm( ) and/or firls( ).

so first divide your complex frequency response into real and imaginary parts:

$$ H_r[k] = |H[k]| \cos\left(\arg\{H[k]\} \right) $$
$$ H_i[k] = |H[k]| \sin\left(\arg\{H[k]\} \right) $$

send both $H_r[k]$ and $H_i[k]$ to the Parks-McClellan or Least Squares optimal FIR designer, but set ftype to 'hilbert' for $H_i[k]$.

sum the two resulting real-valued FIRs to get your result.

lastly, if you want to design an IIR to an arbitrary magnitude and phase spec, then that's a much bigger problem and consideration.  maybe use Prony to design to your impulse response $h[n]$ or use Matt's dissertation that i first learned of now, or use Greg Berchin's FDLS.  i dunno.  no quick solution for you there.


What does it mean to rotate an image by 0.1 degrees?

From the article:

  First the image is rotated by an infinitesimal amount $\alpha$. This is done by computing the rotated coordinates of each pixel and interpolating the gray level values at the new coordinates.

Interpolation is key to very small rotations. It allows to sample the image between the pixel center points as if it was a smooth 2-dimensional function. You can choose the interpolation method in skimage.transform.rotate using the parameter order. It controls the order of the interpolation method which (I think) is an interpolating cardinal B-spline. Higher orders give better results but are slower to compute.

If the tangent vector calculated by rotation, subtraction and division is increasingly non-zero as you go further away from the center, then that is an expected result and you can conclude that some rotation took place.


Your approach will not model the actual channel but can be used to create the received waveform as having passed through a Rayleigh fading channel. I describe considerations to actually create the received waveform using your approach in either flat and frequency selective conditions below, which should give insight into why it would not be a channel simulation.

To simulate a channel

The Rayleigh channel can be simulated with an FIR filter with a significant number of taps that match the sample rate of the waveform and enough taps that the span of the filter exceeds the delay spread of the channel.  As for "correct ways" of simulating a Rayleigh fading channel itself as a FIR filter, please refer to "Jakes' Channel Model" which uses a "sum of sinusoids" approach to modeling the Rayleigh channel. More information on Jakes' model and other similar approaches is available here.

Details on how to proceed with the "Dictionary Approach"

The dictionary would represent the effect of a "flat fading" Rayleigh channel as a "1 tap filter" as you described if your dictionary was generated at your symbol rate and you multiply (not filter!) each transmitted symbol with a sample from your dictionary. In this fashion each symbol would have an independent amplitude and phase from the Rayleigh channel and that fade would be constant over a symbol duration consistent with flat fading. However the fade would be completely independent from one symbol to the next, so in this case would poorly represent an actual channel.

One solution to that problem is to perform a moving average on the generated dictionary to create a new dictionary with some coherence from symbol to symbol. Much better would be to create an upsampled dictionary by $M$ and then do a moving average over $M$ samples. The higher $M$ is, the more accurately you will represent a true Rayleigh channel.  In this approach the transmitted symbols would also need to be interpolated to match the rate in the dictionary. Your transmitted symbols would be interpolated for pulse shaping and other purposes in your modelling, so the final dictionary using this last approach described can be decimated at the output of the moving average filter to match the number of samples required in your waveform for all other reasons to accurately represent your model.

This approach will also give you the "knob" you are looking for to change between frequency selective and flat fading; what is described above by creating a symbol length moving average over $M$ samples, with $M$ samples per symbol, is at the threshold position between flat fading and frequency selective fading.  If you decrease the span of the moving average to be less than the symbol rate, it will become frequency selective fading, as you will have uncorrelated fading samples within your symbol period which is frequency selective fading.

So in summary to create a received waveform that has passed through a Raylegih channel do the following:

1) Create your transmitted waveform (from $N$ complex baseband symbols) that is sufficiently upsampled to model your channel bandwidth of interest (usually for pulse shaping and spectral regrowth or whatever else you may be modelling; this will typically be between 2 and 10 samples per symbol). I will use $K$ to be the total number of samples, and lower case $m$ to be the upsample factor, so then the total number of samples will be $K = Nm$.

2) Create a vector of $KM$ samples of independent complex Gaussian values using the process you described, where $M$ is an additional smoothing factor.

3) Do a moving average over $L$ samples on the vector created in step 2 above.  Choose $L$ based on modeling either flat fading or frequency selective fading where the threshold boundary between the two cases is when the moving average span is equal to the symbol duration ($L \ge mM$ is flat fading and $L < mM$ is frequency selective fading).

4) Decimate the $KM$ samples a the output of the moving average smoothing filter to match the sample rate of your $K$ interpolated transmitted waveform to create the final dictionary. (Decimate by $M$).

5) Finally multiply each of the $K$ samples in the transmitted waveform by each of the $K$ samples in the final dictionary as a sample by sample multiplication to create the received waveform as having gone through a Rayleigh fading channel.

The higher $M$ is, the more accurate the channel result would be but I would think a factor between 6 to 10 would both be easy to manage and be sufficient for creating the coherence from sample to sample. For more frequency selectivity, an even higher $M$ may be needed using this approach. Ultimately with high frequency selectivity each sample is decorrelated and the use of the moving average is unnecessary (the degree of frequency selectivity would be controlled by the interpolation factor of your waveform). I do like the use of the moving average as it gives you the "knob" to control the degree of fading between flat and frequency selective.

To be most accurate (although a lot more complicated), you could optionally vary L over the course of the moving average process within the range consistent that is consistent with your desired fading environment (flat or frequency selective)...this changes the coherence itself over time which would more accurately represent channel conditions, but that may be a second order effect; meaning not critical to your model).


A random process is a collection of random variables, one for each time instant.  Here, the $t$-th random variable is denoted by $X(t)$ where the assumption is that $t \in (-\infty,\infty)$ and what you call the PDF is actually a whole collection of PDFs, one for each of the (uncountably infinitely many) random variables $X(t)$.  Now, we are told for any given (fixed) value of $t$, say $t=0.313$, the random variable $X(t)$ (specifically $X(0.313)$ in our chosen example) is a function (namely a rect function here) of the random variable $T_0$.
$$X(t) = \operatorname{rect}\left(\frac{t-T_0}{T}\right); \qquad
X(0.313) = \operatorname{rect}\left(\frac{0.313-T_0}{T}\right)$$  Now, the rect function can only take on values $0$ or $1$ and so we immediately can say that regardless of the numerical value of $t$ (doesn't have to be
$0.313$) or whatever value the
random variable $T_0$ might have assumed, $X(t)$ can only take on values $0$ or $1$, that is,

  for each value of $t \in (-\infty,\infty)$, $X(t)$ is a Bernoulli random variable.

As such, $X(t)$ does not have a PDF but has a PMF (probability mass function), unless of course you like to use impulses and define what are called discrete density functions to incorporate both notions into one.

The only question remaining is what is the parameter $p$ (the probability that $X(t)$ equals $1$) of this Bernoulli random variable,
and the value of $p$ depends on what value we assume for $t$. That is, the answer is not something like $p = 0.3$ but rather a multi-line equation like
$$p(t) = \begin{cases}\cdots & \text{if}~ t\in \cdots,\\
\cdots & \text{if}~t\in \cdots,\\
0, &\text{otherwise.}
\end{cases}$$
To figure out the details, note that $\operatorname{rect}(\cdot)$ has value $1$ exactly when its argument has value $\leq \frac 12$. Thus, as a function of $t$, $\operatorname{rect}\left(\frac{t-T_0}{T}\right)$ has value $1$ exactly when $$T_0-\frac 12T \leq t \leq T_0+\frac 12T.$$
"Wait a minute", you splutter, "Didn't you say that $t$ had a fixed value
such as $0.313$? What's all this about as a function of $t$??? $t$ is not a variable; it is a fixed value" Well, we need to turn the thing around a bit. Notice that for a fixed value of $t$ (happy now?) $X(t)$ has value $1$ if and only if $T_0$ is such that
$$t - \frac 12T \leq T_0 \leq t+\frac12T,$$
that is, for the chosen fixed value of $t$, and keeping in mind that
$T_0 \in [-T, T]$
\begin{align}
P\{X(t) = 1\} = p(t) &= P\left(t - \frac 12T \leq T_0 \leq t+\frac12T\right)\\
&= P\left(T_0 \in \left(\left[t- \frac 12T, t+ \frac 12T\right]\cap
\big[-T,T\big]\right)\right)
\end{align}
That intersection on the right depends on what we have chosen $t$ to be,
and so we consider various cases. For brevity, denote the interval
$\left[t- \frac 12T, t+ \frac 12T\right]$ as $\mathcal A$ and the interval
$[-T.T]$ as $\mathcal B$

If $\color{blue}{t < -\frac{3T}{2}}$, then $t+\frac 12T < -T$ and so the right end of the interval $\mathcal A$ is smaller than the left end of $\mathcal B$
and so $\mathcal A \cap \mathcal B = \emptyset$. We get that $p(t) = 0$
for $t \leq -\frac{3T}{2}$. Similarly, for $t > \frac{3T}{2}$, the
left end $t-\frac 12T$ exceeds $T$, the right end of $\mathcal B$. So

  $p(t) = 0$ if $|t| > \frac{3T}{2}$.

If $\color{blue}{-\frac{3T}{2} \leq t < -\frac T2}$, then the right end $t+\frac T2$
is in $[-T,0]$ but the left end $t-\frac T2$ is smaller than $-T$, that
is, $\mathcal A \cap \mathcal B =  \left[-T, t+\frac T2\right]$.  Since
$T_0 \sim U[-T,T]$, we get that
$p(t) = \frac{1}{2T}\times \left(t+\frac T2 - (-T)\right)
= \frac{1}{2T}\times \left(t+\frac {3T}{2}\right)$. Note that
$p(t)$ increases from $0$ at $t=\frac{3T}{2}$ to $\frac 12$ as
$t$ approaches $-\frac T2$.  Similarly,
for $\frac T2 < t \leq \frac {3T}{2}$, $p(t)$ decreases from
$\frac 12$ to $0$ linearly.  We have that

  $p(t) = \frac{1}{2T}\times \left(\frac {3T}{2} - |t|\right)$ if
  $\frac T2 < t \leq \frac{3T}{2}$

If $\color{blue}{|t| \leq \frac{T}{2}}$, then both the left end
$t-\frac T2$ and the right end $t+\frac T2$ of $\mathcal A$ are in
$\mathcal B$, and so

  $p(t) = \frac 12$ if $|t| \leq \frac{T}{2}$.

The process in question is not what is called a stationary process and it does not have a PDF (or PMF) but rather a whole bunch of different PDFs (or PMFs), one for each time instant $t$. Read the first part of
this answer to learn a little more about these ideas.


Is there a name for the result of summing the bins of an FFT?

It has the name "first element of the time domain vector times the length of the vector".

reason:

With the DFT of a $N$ long sequence $y[n]$ being

$$Y[k] = \sum_{n=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}\text,$$

your sum of FFT bins is

\begin{align}
\sum_{k=0}^{N -1}Y[k] &= \sum_{k=0}^{N -1}\sum_{n=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}&\text{finite sums can be exchanged in order}\\
&=\sum_{n=0}^{N -1}\sum_{k=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}
&\text{$y[n]$ is constant over all $k$}\\
&=\sum_{n=0}^{N -1}y[n]\underbrace{\sum_{k=0}^{N -1}e^{-j2\pi k\frac{n}N}}_{=\begin{cases}N\cdot1 &n=0\\0&\text{else}\end{cases}}&
\text{sum-orthonogonality of the DFT}\\
&=y[0]N
\end{align}


How to Obtain Mathematically High Frequency and Low Frequency Component Separately Using Bilateral Filter?

One you have convolved your original image with the gaussian kernel (of appropriately chosen variance), then simply subtract the original image from the result of the convolution. This will give you the image with only the high frequency components of your original.


Joint behavior cannot always be deduced from individual behavior. For example, if $X$ and $Y$ are (nondegenerate) random variables with
finite means, then $P\{X<E[X]\}$ and $P\{Y<E[Y]\}$ both are nonzero. This
is because there must be probability mass to the left of the mean (and to the right of the mean too) except in the degenerate case when $X$ has value $E[X]$ with probability $1$. However, $P\{X<E[X],Y<E[Y]\}$ can be zero. Don't believe me? Consider $X\sim$ Bernoulli$(\frac 12)$ and
$Y=1−X$ which is also Bernoulli$(\frac 12)$. Both have means $\frac 12$
but $P\{X < \frac 12, Y < \frac 12\} = 0$.

Turning to ergodicity, one characteristic of an ergodic process is that
over long periods of time (say duration $T$), a typical sample path has value $\alpha$ or less for approximately $F_X(\alpha)\cdot T$ time units where $F_X(\cdot)$ is the first-order CDF of all the random variables
comprising the process. That is, we can get good approximations to the
value of the first-order CDF at $\alpha$
by observing what fraction of time a sample path spends
below level $\alpha$. So, going back to your question, if $\{X(t)\}$
and $\{Y(t)\}$ are ergodic processes, we cannot deduce that the
joint process is also ergodic, for reasons similar to the example
I gave in the first paragraph above. Now, if you impose conditions
such as $\{X(t)\}$ and $\{Y(t)\}$ are independent processes
in addition to being ergodic processes, then Yes, the joint
process is also ergodic.


How can extract the cosine transform formula used for 2D by scipy.fft.dct

I'd take bets that the 2D-DCT is just separable into identical row-wise and column-wise DCTs. Anything else would be surprising.
Anyways, you know the 2D-DCT is a reversible linear operation: It maintains dimensionality, and has only a trivial null space. Therefore,  you can analyze what it does directly from applying it to an arbitrarily chosen basis of your matrix space.


image reconstruction using a set of kernels

I think what you would like to do is something similar to an Auto-Encoder, in which one would append an Encoder to reconstruct back the image. The 'kernel', or more appropriately the intermediary representation is used to synthesize an image, and an appropriate reconstruction objective is minimized to recover the input. On the other hand, you could always train a second upsampling network, which reconstructs the input image, given a 1D/2D feature.


'imread' function in MATLAB gives image data.what this image data is?is that intensity of pixels?and how MATLAB computes this data?

'imread' function just read the RGB(in case of color image) or black/white(in case of black and white image) value of each pixel in image and shows it in a matrix format.
it is not 'calculating' these values its just reading them from the image/path you provide to imread function.


Is there any processing needed to playback audio?

This is unlikely to be a signal processing question. This is about the internals of the dsPic's firmware libraries and would be more suitable on another stackexchange site.

Things that could go wrong:

ADC input and PWM output might have a different data format (for example 12-bit right justified for the ADC, 8-bit for the PWM). You might need to scale / offset the data.
ADC_BUFFER_SIZE and OCPWM_FRAME_SIZE might be of different sizes, requiring you to use an intermediate buffer and counters to keep track of data pouring in/out of this buffer. Or, if the implementation allows it, you might instead read/write partial buffers of a smaller size - identical for both calls.


Short lowpass filter that also extrapolates lines

The most straight-forward solution is a length $N$ FIR filter with impulse response

$$h[n]=\begin{cases}\frac{N}{N-1},&\quad n=0\\
-\frac{1}{N-1},&\quad n=N-1\\0,&\quad\text{otherwise}\end{cases}$$

For $N=3$, the output sample computed from three input samples x = [1 3 5] is $\frac32\cdot 5 - \frac12\cdot 1 =7$, as it should be. Note that this filter is very cheap, it has only two non-zero taps regardless of the value $N$.


Yes, we can explicitely construct such an example: The DFT matrix of dimension $N$ is given by
$$F=\{\exp(-j2\pi ik/N\}_{i=0,\ldots,N-1,k=0,\ldots,N-1}$$

Now, if $T,\Omega\subset\{0,1,\ldots,N-1\}$ the corresponding DFT matrix is given by

$$F_{T\rightarrow\Omega}=\{\exp(-j2\pi ik/N\}_{i\in\Omega,k\in T}$$

We can directly implement this in e.g. Python:

N = 17  # the order/dimension (P in your notation)
L = 8   # the dimension of T, Omega
# choose some random subsets of 0,...,N-1 for T and Omega
T = np.random.choice(np.arange(N), L, replace=False)
Omega = np.random.choice(np.arange(N), L, replace=False)

# Calculate the DFT matrix
t, o = np.meshgrid(T, Omega)
F = np.exp(2j*np.pi*t*o/N)

import numpy.linalg as lg
# Print its rank. It should be always equal
# to L to make F a bijection (i.e. invertible)
print (lg.matrix_rank(F))

When running above code with N=17, the rank is always 8. I've tried setting N=16, and it turns out that it can happen that the rank goes to 7 for some random configuration, i.e. F would not be a bijection.


Understanding IIR , Low Pass, High Pass, Chebyshev, Elliptic Filter?

You did not posed any particular problems and asking for animations? I think that books are way better. We don't know what are your topics, but knowing you did all the studying and have two weeks so I suggest to:

Do the ELEC301x – Discrete Time Signals and Systems course from the Rice University. They do have a good part about design of digital filters, especially about differences between different topologies (including placement of zeros and poles), just watch.
Design few filters in MATLAB, do some filtering and try to understand their characteristics. fvtool is great, especially if you want to visualize anything you need.
Do the exams from past and if you feel confident then just go for a pint!
Read all the questions and answers on DSP SE regarding digital filters ;)

Yep, and after typing into google "digital filter applet" for example you can get:
link1, link2, link3.


You can use interpolation (or equivalent zero-padding) to get your desired sub-divided 512 frequency points from your 220 FFT bins.  Use a Sinc interpolation kernel for a rectangular window, or the transform of the window for other windows.  If you want more real frequency resolution (as in clear separation of multiple peaks present in the frequency spectrum), you will have to use (overlap or "smear") more than 10 mS of data (perhaps over 50 mS).


Can state covariance matrix in Unscented Kalman filter contain negative values?

actually a covariance matrix can have negative values off the diagonal.  variables can be negatively correlated, but the diagonal terms are variances which must be positive.

The update you show can be a problem which is why updates of $P^{-1}$ are preferred.


How do I rotate a warping grid without changing the warped image's orientation?

After three and a half days of cursing, a couple of beers with a friend helped solve the problem. My mistake was that the reordering of points needed to be not exactly transpose, but close. So:

To rotate a grid 90 degrees counter-clockwise:

transpose the grid and reflect it along the horizontal central line
apply the transformation

x = y'

y = -x',

where primes denote old coordinates.

To rotate a grid 180 degrees

reflect the grid around its horizontal and around it's vertical axis
apply the transformation

x = - x'

y = - y'

To rotate a grid 270 degrees:

transpose the grid and reflect it along the vertical central line
apply the transformation

x = - y'

y = x'


Determine poles and zeros of a specific filter design

There's a good reason why band pass filters usually have an even filter order. As you may know, real-valued filters have either real-valued poles and zeros, or their poles and zeros occur in complex-conjugate pairs. If the order is odd, there must always be at least one real-valued pole and zero. In the case of a band pass filter, a real-valued pole is not very useful, because poles create peaks in the frequency response, but a real-valued pole can only create a peak at DC (frequency zero), or at Nyquist (the maximum frequency given the sampling rate). However, a band pass must attenuate the signal at DC and at Nyquist, so peaks at those frequencies are undesirable.

So a reasonable choice would be to place a complex conjugate pole pair at (or close to) the center frequency of the band pass filter, and leave the remaining real-valued pole at $z=0$ (where it doesn't affect the magnitude response). The three zeros should all be placed on the unit circle to create frequency response zeros. A reasonable choice would be to have one zero at DC and a complex pair of zeros close to Nyquist.

The figure below shows the magnitude response and the pole zero diagram of a possible choice of poles and zeros for a third order band pass filter with the given pass band.


When using a constant tone audio beacon, beware of room echoes causing multi-path interference and distortion, especially around the leading and trailing portions of your received waveforms.

Try using a frequency sweep instead of a constant tone for your transmit waveform.  This might provide you with a sharper correlation peak that is less likely to have duplicate peaks.


imfindcircles can take in other "Name, Value" pairs as parameters. In your case, a really helpful one is the Sensitivity parameter. I would change your first line to include this (also note that 'bright' is the default value for 'ObjectPolarity', so I took it out). The value `0.95' here can be tuned to get more or fewer circles identified.

[centers, radii] = imfindcircles(W, [30 65], 'Sensitivity', 0.95);

Since you mentioned circles may vary between pictures, it would be beneficial use some other parameters to get better results across the set of images. Check out the MATLAB documentation online for specifics.


A-GPS will give you ephemeris data. That will be used to update your initial guess for the frequency of the downlink signal, indeed. The relevant remaining uncertainty  is in the local oscillator of the GPS receiver itself. And that can be large compared to Dopplers!

So, yes, assuming the frequency probability density to be symmetrical convex around the nominal frequency, that would put an ML estimator at the doppler frequency.

But: your GPS receiver typically doesn't run the first time. It already has acquired estimates about its own frequency error in the past, and together with parameters like current temperature and the time that has passed since the last successful fix, your initial frequency error estimate doesn't have to be 0; for example, if I remember correctly (it's really a bit of a fuzzy memory, please check the source code) RTKlib uses an extended Kalman predictor for such purposes.

Also, the highest likelihood estimate for an initial frequency error is not necessarily one that makes for the best starting point in a search algorithm!

For example: Due to the ambiguities in multiples of (chip rate)/(chip sequence length), it might be cleverer to first test a frequency hypothesis that's intentionally a bit off, because it can give a drastically improved expected convergence rate when looking for a frequency error optimum – in the end, this is a relatively complex optimization problem that's not globally convex.


Help Understanding the Robert Cross Operator

What I do not get is Wikipedia. Why does it go on finding the gradient?

You subtract two nearby pixel values. That is a gradient, or as close as you can get on a discrete image.

  Also that looks like the Pythagorean Theorem not a gradient equation.

As the book text say, "it is often better to consider the two templates as providing the components of an edge vector". If you want the length of a vector, you use the Pythagorean Theorem.

The better question is: why use the max or the sum of the absolute values, or: what's the difference?

To show that, it's best to try the filter on a rotationally symmetric template:

that way, we can see how it acts on different edge angles:

As you can see, adding the absolute values of the gradient vector (this is sometimes called the "L1-Norm") gives larger edge strengths for horizontal and vertical edges than diagonal edges. Taking the max of the two absolute values of the gradient vector (also known as the "L-infinity Norm") gives slightly larger edge values for diagonal edges. Using the ordinary Euclidean norm (i.e. the Pythagorean Theorem) gives exactly the same edge strength for any angle, which is of course the right result.

The only downside (I'm aware of) is that this requires two multiplications and a square root for every pixel, which used to be very expensive operations when the Roberts operator was invented. They're still expensive on mobile or embedded platforms today, but on a PC, you might not even notice the difference.


Determinant of Hessian approximation (SURF)

I think $w$ is a factor that makes the ratio of  $D_{xy}$ to $D_{yy}$ the same as $L_{xy}$ to $L_{yy}$. This is so the value of the determinant for the simplified kernel roughly matches that of the continuous version.
You could write it that way. But if you leave it as $wD_{xy}$ then if used in other expressions you don't have to muck around with adding square roots. e.g. $(wD_{xx})^3$ expressed your suggested way would be $\sqrt{w}^3D_{xx}^3$


Equation for impulse train as sum of complex exponentials

The summation on the left side of your equation represents a time domain discrete time periodic signal $x[n]$ whose period is N.
$$ x[n] = \sum_{k=-\infty}^{\infty} \delta[n-kN] $$

And the summation on the right side of your equation is just the same (periodic) signal $x[n]$ represented via its Discrete Fourier Series DFS synthesis form which is an alternative representation for periodic signals
$$ x[n] = \sum_{k=0}^{N-1} a_k e^{j \frac{2 \pi}{N} k n} $$

where the DFS coefficients $a_k$ are all $1/N$ found via DFS analysis equation:
$$ a_k = \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-j \frac{2 \pi}{N} n k} $$

I hope you can see why all $a_k$ are $1/N$ (it can be easily seen from the properties of impulse inside the summation...)


Computational burden of EMD/Huang-Hilbert vs wavelet

There are relationships between fixed wavelets and EMD, for certain classes of signals, cf. Empirical Mode Decompositions as data-driven wavelet-like wavelet-like expansions. A related relation between a fixed basis and an adaptive one can be found between Fourier transforms and Karhunen-Loève/PCA type decompositions.

Of course, when the vectors used for representations are learnt instead of being fixed beforehand, they are likely to require more computations. And are likely to be more difficult to implement online. Plus, DWT are often critically sampled, while EMD are redundant (hence greedier). Yet, some papers (see below) evaluate it proportionally to the burden of the FFT. Roughly to that of a redundant stationary wavelet transform.

Then, for non constrained versions, EMD can be quite sensitive to border effect, local spike instability, etc. As a result, from one signal to a similar other, a different number of IMFs can be necessary, which is not the case of a fixed-level DWT (but you can find that with adaptive wavelet packets).

Further, each can be sensitive to the features you can extract.

So yes, it is not so surprising to me, with experience on seismic data and combustion signals with sporadic instabilities. Yet, it can work on some data.

Finally, some online implementations exist:

On-line empirical mode decomposition biomedical microprocessor for Hilbert Huang transform, 2011

  On-chip implementation of Hilbert-Huang transform (HHT) has great
  impact to analyze the non-linear and non-stationary biomedical signals
  on wearable or implantable sensors for the real-time applications.
  Empirical mode decomposition (EMD) is the key component for the HHT
  processor. In tradition, EMD is usually performed after the collection
  of a large window of signals, and the long latency may not be feasible
  for the real-time applications. In this work, the architecture of
  on-line EMD for biomedical signals is proposed. The on-line
  interpolation method with data reuse as well as component and
  iteration loop decomposition is applied to obtain low latency and low
  hardware cost. The first chip of EMD processor is fabricated in UMC
  90nm LL process and consumes 57.3µW.

Online Empirical Mode Decomposition, 2017

  The success of Empirical Mode Decomposition (EMD) resides in its
  practical approach to dissect non-stationary data. EMD repetitively
  goes through the entire data span to iteratively extract Intrinsic
  Mode Functions (IMFs). This approach, however, is not suitable for
  data stream as the entire data set has to be reconsidered every time a
  new point is added. To overcome this, we propose Online EMD, an
  algorithm that extracts IMFs on the fly. The two key elements of
  Online EMD are a sliding window to compute local IMFs, and a stitching
  procedure to gradually append local IMFs to the final result. Using
  synthetic data we show that the decomposition quality of Online EMD is
  similar to classical EMD. We also present results obtained with a real
  data set to expose the practical advantages of Online EMD when dealing
  with data stream or large data set.

On the computational complexity of the empirical mode decomposition algorithm, 2014

  It has been claimed that the empirical mode decomposition (EMD) and
  its improved version the ensemble EMD (EEMD) are computation
  intensive. In this study we will prove that the time complexity of the
  EMD/EEMD, which has never been analyzed before, is actually equivalent
  to that of the Fourier Transform. Numerical examples are presented to
  verify that EMD/EEMD is, in fact, a computationally efficient method.


Identifying the correct source of packets

This is a very broad topic, and a full answer requires a few textbook chapters. A few pointers:

This problem is not solved in the physical layer; in other words, don't think in terms of frequencies, bits and packets.
The physical layer can handle some problems such as multipath, error correction and detection, and synchronization. The problem you're referring to is solved in the upper layers in the stack, and the solution depends on the particular network implementation.
For example, in WiFi there is an authentication protocol that ensures (with high probability) that every received packet comes from the actual source. A man-in-the-middle wouldn't be able to generate packets that would be recognized as valid by the receiver; again, this would happen somewhere in the network stack, not at the physical layer.


transfer function and 'causal' signal - evaluate transfer function or use z-transform of input?

For a discrete-time LTI system with transfer function $H(z)$, the response to $x[n]=z_1^n$ equals $y[n]=H(z_1)z_1^n$ if $z_1$ is inside the region of convergence of $H(z)$. This relationship holds whether or not the system is causal or stable.

From this relationship it follows that if the system is real-valued and stable, i.e., if its impulse response is real-valued and if $H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}$ exists, the response to a sinusoidal input $x[n]=\cos(n\omega_0)$ is given by $y[n]=|H(e^{j\omega_0})|\cos[n\omega_0+\phi(\omega_0)]$.

I think your confusion comes from assuming something you call a "causal context". It doesn't matter if a system is causal or not for above relationships to hold. You can apply a sinusoidal input signal to a causal system, and as long as the system is real-valued and stable, the output will be given by the relationship mentioned above. If the input is a sinusoid switched on at a finite $n=n_0$, the system's response can be computed by solving the convolution sum or, equivalently, by using the $\mathcal{Z}$-transform.


Quantifying goodness of amplitude estimation

Your understanding of the FT as several extremely narrow bandpass filters is fine. However, you mix up the terms SNR and noise variance. SNR is defined by

$$SNR=\frac{E(s(t)^2)}{E(n(t)^2)}$$

where $E[n(t)^2]=\int_\mathbb{R}S(f)df$ due to the Wiener-Khintschin Theorem ($S(f)$ is the PSD of the noise). So, if you have a measured signal $x(t)=s(t)+n(t)$ with $s(t)=A\sin(2\pi f_0t)$ and $n(t)$ is bandlimited AWGN with bandwidth $B$ and density $N_0$, the SNR decreases when the noise bandwidth increases. So, in order to have a good SNR, you should keep the bandwidth as small as possible (as narrow, as your $s(t)$ permits).

However, your problem is to estimate the amplitude $A$ of $s(t)=A\sin(2\pi f_0t)$ where $f_0$ is known to you. Optimally, you take the Fourier transform of your signal and see what is the value of $X(f)$ at the frequency $f_0$. What will you see there? Does this measurement depend on SNR? Does this measurement actually depend on the noise bandwidth? No, the measurement $X(f_0)$ is independent of the noise bandwidth (as long as the noise has contribution at $f_0$). So, by performing bandpass filtering, increasing/decreasing bandwidth you will not get any better estimates.

The only thing which you can do is to increase the duration of the measurement window. If you have a finite length measurement of duration $T$, your frequency domain essentially has resolution $1/T$. Hence, on each frequency bin (assuming the DFT) the noise variance is $N_0/T$ (assuming AWGN). So, the longer the measurement, the better the estimate.

You can also understand this from a different point: What you are interested in, is the value

$$\rho=\frac{\int_{f_0-\Delta_f}^{f_0+\Delta_f}S(f)df}{\int_{f_0-\Delta_f}^{f_0+\Delta_f}N(f)df}$$

where $2\Delta_f$ describes the frequency domain resolution, $S(f)$ is the signal spectrum and $N(f)$ is the noise spectrum. The numerator is always the same (independent of $\Delta_f$) since $S(f)$ is just a (weighted) Dirac at $f_0$. But, the denominator depends on $\Delta_f$.


3D reconstruction of an image based on depth of intensity levels?

This is an extremely difficult task, one which is a very active line of research. I've managed to find a semi-recent paper on the subject. I won't go into the details, but here's a few things that this paper found that can be used.

Textures: If you see a given texture, they can help you figure out how far something is away. The wood texture here would be quite ideal for this.
Fixed object sizes: Identify common objects, see how large they are, and use them to help with depth analysis. The phone book, soda bottle, etc should work well with this.
Focus: If you know how the picture was take, focus can help you identify ranges.
Shadows- If you know where the light was coming from, shadows can help you identify how tall objects are, and thus how far they are away.


Map time/scene differences between two similar videos

This was an interesting read and it looks like you're not too far off! I think that what you have left to do is divide one of the signals into intervals (let's call it signal 1) and cross-correlate each segment of signal 1 with the entirety of signal 2. The peak in the cross-correlation will indicate where there is the strongest match for that segment in signal 2. You then calculate the difference in frames between where the segment was in signal 1 and where its strongest match is in signal 2.
The third step in the post you mentioned appears to be Kernel Density Estimation. The author placed a Gaussian at the center of each segment from signal 1. The standard deviation of the Gaussian, I assume is one-sixth the segment length, and the amplitude of the Gaussian is the frame difference. The Gaussian should extend from time 0 to the length of signal 2. He then summed all the Gaussians to obtained a continuous function.
I think you can achieve the same thing if you do the following: once you have the frame difference as a function of time, you convolve that function with a Gaussian whose standard deviation is one-sixth the segment length you used to divide signal 1.
I hope this helped!


Why is the first axis of STFT frequency bin

It is really a matter of convention of the communities that use STFT. Medical devices tend to have time scrolling left to right or over writing left to right. For passive SONAR, most displays that Im familiar with have time water falling down and frequency increasing from left to right. There are exceptions as well. I've also seen amplitudes as marching mountains, color coded, and intensity coded as well. Much depends on the amount of time history that should be displayed and the time duration of signal events that are of interest.  There are also perceptual consideration based on human vision that come into play and effect how STFT data is displayed. There are also places where an STFT display is the input to image processing algorithm.


What would the target matrix to train Neural Network?

Example with the following numbers (I use random numbers):

There are 6 possible classes (face expressions in your case).
You have 2000 examples (lets say 2000 photos of faces from which you know the correct face expression).
There are 30 features. From each example you have extracted the 30 features and know the correct class.

With these numbers, the input matrix would be a matrix of $30 \times 2000$, where each column $j$ is an example, and it has the 30 features of it. Hence, the position $(i,j)$ contains the feature $i$ of the example $j$.

The target matrix would be a matrix of $6\times 2000$. Where each column $j$ represents the correct class of the example $j$. The column is formed by zeros, and it has a $1$ at the row that indicates the class. If we have 6 classes, and the first example corresponds to class 3, the first column of the target matrix will be $[0 \ 0 \ 1\  0\  0\  0]$. If the second example corresponds to class 5, the second column of the target matrix will be: $[0 \ 0 \ 0 \ 0 \ 1\  0]$, and so on.

More info here.


Convolution effects width of the signal?

If $y(t)$ is the signal resulting from the convolution of $x_1(t)$ with $x_2(t)$ then it will have the same bandwidth as $x_1(t)$. $x_1(t)$ has (presumably one-sided) frequency support in $[0, 500]$ Hz and $x_2(t)$ has frequency support in $[0, 1000]$ Hz. So, when performing pointwise multiplication in the frequency domain, the higher frequency components in $x_2(t)$ between $500$ Hz and $1000$ Hz will be multiplied by zero (the value of $X_1(f)$ at those frequencies). The Nyquist frequency required to capture $y(t)$ without distortion should be any frequency greater than $1000$ Hz (twice the highest frequency in the signal).


How to calibrate three orthogonal cameras？

Use an asymmetric calibration plate. Such as the new circle grid in OpenCV. Check here. You do not need to find the correspondences in the gui, but rather run the calibration directly. Also, do not forget to initialize the optimization problem using the median of the obtained poses.

To calibrate multiple view setups, I would also recommend Multi Camera Self Calibration tool. It uses a simple laser dot and is much easier to calibrate, especially when the number of cameras increase. Of course, pay attention to the synchronization.


Peak to Average Power Ratio of an exponentially increasing signal

This feels a lot like there's quite a bit of homework you should do about this, but let's give you a start:

Write down your signal explicitly as

$$s(t) = p(t) \cdot e^{\lambda t}\text,$$

with $p$ being your periodic oscillation. You sample that, and divide it into vectors. Each vector has $N=128$ samples. Let's write it this:

\begin{align}
s_k[n] &= s\left(\frac{n+N\cdot k}{f_\text{sample}}\right),& n = 0,\ldots, N-1\\
&= p\left(\frac{n+N\cdot k}{f_\text{sample}}\right) \cdot e^{\lambda\frac{n+N\cdot k}{f_\text{sample}}}
\end{align}

where $k$ gives us the frame number.

Your question is how things develop when you increase the frame number!

So, consider the quotient of the next frame and the current frame:

\begin{align}
\frac{s_{k+1}[n]}{s_k[n]}
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot \frac{e^{\lambda\frac{n+N\cdot (k+1)}{f_\text{sample}}}}{e^{\lambda\frac{n+N\cdot k}{f_\text{sample}}}}\\
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot e^{\lambda\frac{n+N\cdot (k+1)}{f_\text{sample}}-\lambda\frac{n+N\cdot k}{f_\text{sample}}}\\
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot e^{\lambda\frac{N}{f_\text{sample}}}\\
\end{align}

Now, it'd be a nice simplification if we can just restrict the oscillations' period to something that is an divisor of $N$. That especially means that $p()$ becomes $N$-periodic (under the given sampling rate).

\begin{align}
&= 1 \cdot e^{\lambda\frac{N}{f_\text{sample}}}\\
&=\text{const.}
\end{align}

In other words, the next frame is identical to the previous frame, aside from a constant factor.

Since scaling with a constant factor doesn't change the PAPR at all, this was easy to answer for oscillations that have a frequency that divides the sampling rate, only. The rest isn't really unusual (it's leakage, essentially), but the math is slightly more annoying, so I'll leave it to you.


I'm not sure why the problem occurs for you with the x22 example, but I believe the problem is because the deconvolution is not being performed over the right length.

Deconvolution of an FIR (inverting an FIR) requires, in general, convolution using an infinite duration impulse response.  So just using the default lengths in performing the deconvolution using an FFT might not generate the correct result.

I've implemented the same problem in R, but the last (additional) row of plots shows the results obtained if I do the deconvolution over a longer signal than the default. This generates the right results.

As you can see, my implementation generates bizarre results for the x21 example rather than the x22 example in yours.

R Code Below

x1 <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0)
x21 <- c(0,1,1,1,0)
x22 <- c(0,1,1,1,0,0)
x23 <- c(0,1,1,1,0,0,0)

my_conv <- function(x,y)
{
  out_length <- length(x)+length(y)-1
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) * fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

my_deconv <- function(x,y, out_length=0)
{
  if (out_length == 0)
  {
    out_length <- length(x)+length(y)-1
  }
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) / fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

par(mfrow=c(5,3))
plot(x1)
plot(x1)
plot(x1)
plot(x21)
plot(x22)
plot(x23)
y1 <- my_conv(x1,x21)
plot(y1)
y2 <- my_conv(x1,x22)
plot(y2)
y3 <- my_conv(x1,x23)
plot(y3)
plot(my_deconv(y1,x1))
plot(my_deconv(y2,x1))
plot(my_deconv(y3,x1))
plot(my_deconv(y1,x1,50))
plot(my_deconv(y2,x1,50))
plot(my_deconv(y3,x1,50))


Where is the flaw in this derivation of the DTFT of the unit step sequence $u[n]$?

There are infinitely many signals that make the following equality hold: $$y[n]-y[n-1]=\delta[n] \qquad (1)$$
The only thing that matters is that $y[0]-y[-1]=1$, and then the rest of the coefficients of $y$ can be determined under the restriction that Eq. $(1)$ states (i.e. the substraction of consecutive samples must be $0$ for $n\neq 0$).
In other words, Eq. $(1)$ will be achieved by any signal $y[n]$ such that $$y[0]=y[-1]+1 \land y[n]=y[n-1] \ \forall n\neq0$$
Another way to see this is that any function that is basically $u[n]$ with an offset (a constant value added) will satisfy $(1)$. This explains the statement made by robert bristow-johnson in his answer: differentiators destroy this information (such as taking a derivative in continuous time destroys evidence of any constant value in the original function).

To sum up, I believe that the proof is flawed because the procedure followed could use any function of the form $u[n]+C$ with $C\in\mathbb{R}$, and this would lead to many functions having the same Fourier transform, which is indeed wrong as the Fourier transform is a bijection. Maybe the author deliberately decided to ignore anything related to DC values, conscious that in order to show that $F(\omega)$ is the DTFT of $f[n]$ he would need the accumulation property (whose most popular proof is derived from the DTFT of the unit step - ergo, a pretty circular proof). The proof is not strictly wrong, as everything it states (the formulae for $F(\omega)$ and $U(\omega)$, the decomposition of the unit step, the difference equation) is true, but it would require the accumulation property to show why $F(\omega)$ doesn't have any Dirac deltas.


Given the input and the output, how to determine the impulse response?

$\mathcal Z$-transforming the input and the output and then dividing them, we obtain the transfer function

$$\begin{array}{rl} H (z) &= \dfrac{1 - z^{-1} - 4 z^{-2} + 4 z^{-3}}{1 - 3 z^{-1} + 2 z^{-2}}\\\\ &= \dfrac{z^3 - z^2 - 4 z + 4}{z (z^2 - 3 z + 2)}\\\\ &= \dfrac{(z-1)  (z-2) (z+2)}{z (z-1) (z-2)}\\\\ &= \dfrac{z+2}{z} = \color{blue}{1 + 2 z^{-1}}\end{array}$$

Computing the inverse $\mathcal Z$-transform of $H (z)$, we finally obtain the impulse response.


Getting non smooth spectrum after performing fft on numerically obtained acceleration data

Most of the time, a non-parametric PSD estimation based on FFT alone (called a Periodogram) will provide a random looking spectrum of the numerical data.

To reduce those random variations, or to get a smoother looking estimate, you can do one of the two things:

Use Periodogram averaging such as Welch's method.
Use a model based parametric PSD estimate instead.

Note that the first approach will trade-off spectral resolution for reduced variance, whereas the second method will be high resolution but requires a consistent knowledge of the associated physical model of the process.


TL;DR:  Subspaces are low-dimensional, linear portions of the entire signal space that are expected to contain (or be close to) a large part of the observable and useful signals  or transformations thereof, with additional tools that allow us to compute interesting things on the data
We are  given a set of data. To manipulate them more easily, it is common to embed them, or represented them in a well-adapted mathematical structure (from the plenty of structures  we have in algebra or geometry), to perform operations, prove things, develop algorithms, etc. For instance in channel coding, group or ring structures can be better adapted. In a domain called mathematical morphology, one uses lattices.
Here, for standard signals or images, we often suppose a linear structure: signals can be weighted, added: $\alpha x+ \beta y$. This is the base for linear systems, like traditional windowing, filtering (convolution), differentiating, etc.
So, a mathematical  structure of choice lies in vector spaces. Vector spaces equipped with tools: a dot product (that can be used to compare data), a norm (to mesure distances). These tools help us compute. Indeed, energy minimization and linearity are strongly related.
Then, a data of $N$ samples naturally lives in the classical linear space of $N$ dimension. It is quite big (think of million-pixel images). It contains an awful lot of other "uninteresting" data: any $N$ dimensional "random" vector. Most of them are and will  never be observed, have no meaning, etc.
The reasonable quantity of  signals that you can record, up to variations, is very small relatively to the big space.
Even more, we are often interested in structured information. So if you subtract noise effects, unimportant variations, the proportion of useful signals is very very tiny within the whole potential signal space.
One very useful hypothesis (heuristic, to help discover) is that those interesting signals live close together, or at least along regions of the space that "make sense". An example: suppose that some extraterrestrial intelligence has no other detection system than  a very precise dog detector. They will get, across the Solar system, almost nothing, except many points located on something vaguely looking like a sphere, with large empty spaces (oceans), and  sometimes very concentrated (urban areas). And the point cloud moves around a center, with a constant periodicity, and rotating on itself. Those aliens have discovered something!
Anyway, the partial-sphere looking point cloud is interpretable... maybe a planet?
So, our dog point cloud could have been fully 3D, but they are concentrated on a 2D surface (lower dimension), that seems relatively regular (in altitude) and smooth: most dogs live at intermediate altitudes.
These smooth low-dimensional parts of space are sometimes called smooth manifolds or varieties. Their structure and operators allow to compute things. For instance: distances, distributions, etc. Inter-dog distances make more sense when computed along the Earth surface (in spherical 2D coordinates) than directly through the planet with the standard 3D norm! But this  can still be complicated to deal with. Let us simplify this a bit more.
Looking a little closer, the dog points are almost located on close-to-flat surfaces: countries, even continents. Those flat surfaces are portions of linear (or affine) subspaces. Still, you can now compute inter-dog distance, more easily, and design an algorithm for dog matching that will make you rich.
The story continues a bit. Sometimes, natural data does not assemble around a clear structure, directly. Unveiling this inherent structure is at the core of DSP. To help us in this direction, we can resort to data transformations to concentrate them better (Fourier, time-frequency, wavelets), filtering.
And if we find a suitable subspace, most algorithms become simpler, more tractable, and so on: adaptive filtering, denoising, matching.
[ADDITION] A typical use is the following: a signal can be better concentrated with a well-chosen orthogonal transform. In the meantime, a zero-mean random Gaussian noise remains Gaussian under an orthogonal transformation. Typically, the covariance matrix can be diagonalized. If you sort the eigenvalues in decreasing ordre, the smallest ones tend to flatten (they correspond to noise), and the highest more or less correspond to the signal. Hence, by thresholding the eigenvalues, it because possible to remove the noise.


what is the difference between $X(j\omega)$ and $X(\omega)$ notation

Both notations are common and correct. As pointed out by Yuri Nenakhov, the advantage of the argument $j\omega$ is that it coincides with the complex (Laplace transform) variable $s$ when its real-part is zero. Note that in the complex $s$-plane the frequency axis is the imaginary axis. So $j\omega$ has nothing to do with complex frequency (which makes no sense).

So if the Laplace transform $X(s)$ of a signal $x(t)$ exists, and if the imaginary axis is inside its region of convergence, then the Fourier transform is obtained by setting $s=j\omega$.

Note that this does not work in general! In general you can't get the Fourier transform by replacing $s$ with $j\omega$ and vice versa. Two conditions must be satisfied in order for this to lead to a correct result:

Both transforms must exist (in the sense that the corresponding signal $x(t)$ has a Laplace transform and a Fourier transform).
The imaginary axis $s=j\omega$ must be inside the region of convergence of the Laplace transform.

An example where replacing $s$ by $j\omega$ doesn't work, even though both transforms exist, is the step function:

$$\begin{align}&x(t)=u(t)\\\text{Laplace transform: }&X(s)=\frac{1}{s}\\
\text{Fourier transform: }&\hat{X}(j\omega)=\pi\delta(\omega)+\frac{1}{j\omega}\neq X(s)|_{s=j\omega}\end{align}$$


Motion can be described using displacement with time, velocity or acceleration. To get these quantities from the accelerometer data the following link may help.
http://www.chrobotics.com/library/accel-position-velocity


Signal segmentation into states ( Stop/mov)

It seems like you already got all the right components. Your goal is to detect segments where speed!=0 in the data. Since speed is the derivative of position, you could transform the data into the first difference series - simply subtract the preceding point from each data point, i.e. $x_t - x_{t-1}$, and you'll have a crude estimate of the derivative.
The series of differences then can be segmented with any standard method for peak or changepoint detection. In particular, since you know the background value (zero), and the speed appears roughly constant during each movement episode, it seems that epidemic detection of changes in mean should be able to solve this nicely. I am aware of a couple packages for this in R, but it's also relatively easy to implement it yourself if it's not available for your programming language (just search for "epidemic changepoint detection"). Multivariate detectors in the same framework are also available, if you were interested in incorporating several sensors to detect movement, for example.


Closed form of $\mathcal Z$-transform : decomposition signal $x(n)$

Or you can proceed this way: \begin{align}\mathcal Z\left\{x(n)\right\}= X(z) &= \sum_{n=-\infty}^\infty x(n)z^{-n}\\
& = 0 + \sum_{n=0}^{2N} \left|n - N\right|z^{-n} + 0\\
& = \sum_{n=0}^{N} \left(-n + N\right)z^{-n} + \sum_{n=N}^{2N} \left(n - N\right)z^{-n}\\
& = N\sum_{n=0}^{N}z^{-n} - N\sum_{n=N}^{2N}z^{-n} - \sum_{n=0}^{N} nz^{-n} + \sum_{n=N}^{2N} nz^{-n}
\end{align}

And changing the boundaries for the second and fourth term, putting $k = n - N$, you get :

$$X(z) = N\sum_{n=0}^{N}z^{-n} - Nz^{-N}\sum_{k=0}^{N}z^{-k} - \sum_{n=0}^{N} nz^{-n} + z^{-N}\sum_{k=0}^{N} kz^{-k} $$

You can continue from there.


Benefits of oversampling in a realtime noise-cancellation system

The system has zero time delay (sampling, calculations, etc. are instantaneous)

This is an invalid assumption. Almost every digital system has at least one sample delay. DAC and ADC or clocked off the same master clock. So at sample time n, the ADC is acquiring sample n while at the same time the DAC is outputting sample n-1. The only way around this is using a much higher master clock, donwsampling that clock and clocking the A/D and D/A on different phases. However that's equivalent to running at a higher sample rate in the first place.

So your ability to cancel is given by how much the amplitude of the sine wave changes during 1 sample. The residual signal is simply given by $$e(t) = \sin(f_0 \cdot t) - \sin(f_0 \cdot (t-T_s)) = 2 \cdot \cos(f_0 \cdot (t- T_s/2)) \cdot \sin(f_0 \cdot T_s/2)$$
 where $T_s = 1/f_s$ is the sample time. The error is simply a function of the ratio of signal frequency and sample rate. For $f_0 << f_s$ the error is very small.

The only way to get around this is to crank up the sample rate or use some a priori knowledge of the signal to do some prediction and speculative cancellation.


In case of Complex DFT spectrum, why the x axis range from N/2 to N point mean a negative frequency?

it's because the DFT is periodic with period $N$.

$$  X[k+N] = X[k] \qquad \forall k \in \mathbb{Z} $$

this periodic extension is also true about the discrete-time signal $x[n]$

$$  x[n+N] = x[n] \qquad \forall n \in \mathbb{Z} $$

doesn't matter if either time-domain or frequency-domain representations are complex or not.

that means for $-\frac{N}{2}\le k < 0 $ (the negative frequencies), then $X[k] = X[k+N]$.  the explicit negative frequency component $X[k]$ has a copy $X[k+N]$, but that is in latter half of the DFT result.

If  $-\frac{N}{2}\le k < 0 $, then  $\frac{N}{2}\le k+N < N $.


Is there a problem with using the WAV format for non-audio signals?

I think WAV will be a good format for your data as it supports:

Integer Hz sampling frequencies (like your 1000 Hz).
Integer sample formats (like 16-bit, 24-bit, 32-bit).
A floating point sample format (32-bit IEEE float).
Multiple channels.
Seeking (quickly reading the file at a specific time point in the audio).

Other nice things about WAV:

You can use FLAC to compress WAV files for storage and transfer.
WAV files can also contain non-PCM data (such as markers and metadata).


They are system norms, a metric that you can compare two different systems in terms of their generalized gain and spread. You can look these up no need for attaching physical motivation.
You don't. They are for assesing and using in the optimization programs.
Yes, but only if you know what you are doing.
Can be. But again, they are for assessment and optimization. If you hit the minimum with a PID controller then it won't make any difference. But hitting the minimum is pretty much unlikely.
Unfortunately no.
see 5 (you can pull some tricks by bounding the nonlinearity with a linear system etc. or for very limited nonlinear operators you can apply them)
Not applicable. They are different things.
They are different norms. But both of them are norms. Read about them.
See the documentation of hinfsyn of matlab
Yes


About speech preprocessing in mobile phones

Yes, the cellular phones use various forms of compression to convert the captured analog audio (speech) into the digital bitstream for transmission through 2G/3G/.../.

The specific method used depends on the GSM version which might dictate its own bandwidth constraint and backward or forward compatibility issues into the audio encoding stage.

Most typically a 2G/3G system utilizes a variant of voice code based on a linear predictive coding (LPC) scheme which converts the analog sound captured by the microphone into a set of coefficients (parameters) according to the vocal tract model of the speaker. Unlike PCM, mp3, ogg, ra, etc., this is not a waveform coder but a parametric voice coder therefore.

The speech is reproduced at the recevier side by using a similar code excited LPC (CELP) which uses various excitation signals for the vocal tract model to reproduce the input signal as natural as possible.

You can find more details in the literature by the keywords vocoder, LPC, CELP, GSM speech coding etc...


Is it possible to generate Hamming code with minimum distance $d$ which can correct more than 1 bit error?

You cannot detect 2-bit errors (reliably) with a Hamming code distance of 3 if you are correcting 1-bit errors.

At any rate, you can design Hamming codes correcting two errors.  Let's assume 32bit of data and n Hamming code bits.  The states we then need to be able to encode is
0 bit error = 1 state
1 bit error = 32+n states
2 bit error = (32+n)(31+n)/2 states
For n=10, this amounts to 1+32+10+21*41 = 904 states that can be encoded with 10 bits of redundancy.  Depending on how complex you want your correction machinery to get, you might want to use more bits in order to make multiple error correction less complex.


Help required with BM3D Image denoising algorithm implementation

Step 3 is not correct. The first part is ok (each denoised result is weighted by some likelihood of correctness, that was obtained previously from the variance of the data), but you don't have to add them up. Each slice from the 3D array is a 2D patch that is back-projected to the resulting image.

The sum happens when 2 back-projected patches do overlap in the final image domain.

About back-projection

Back-projection is the process of taking a 2D slice from the denoised stack and put it back in image space.
This is required because unlike most variants of NL-means BM3D denoises patches and not the central pixel of a patch.

So, you have 2 images when implementing BM3D:

the original (noisy) one;
the result (denoised) image.

The spatial location of a patch should not change during the denoising process, otherwise the image content would be destroyed. So, each slice from the denoised patch stack goes back to its original place.

Now, one issue remains: what to do when back-projected patches overlap? What value should we keep for a single pixel?
The solution chosen by BM3D is that each pixel contains a weighted mean of the patches that overlap in each pix.

To implement this, you can simply pre-multiply each slice of the denoised stack by its weight, add a third image to store the sum of the weights inside each pixel, sum for each pixel its value from all the denoised patches that contain it, and in a final step (once all the patch stacks haves been processed) divide the result image by the weight image.

Image processing main loop (outline)

There is no "unleft" pixel in the final result. The outline of the main loop is as follows:

for every pixel of the image, form a reference patch;
for a given reference patch, find n <= N (B=16 in your case) patches whose distance to the reference patch is less than some threshold and make a stack;
denoise the stack;
back project the patches in the denoised stack;
back to 1) for the next pixel).

Thus, you always have one stack for a given patch (degenerate case: only 1 patch in the stack), and a pixel is always involved in several patches.


mne raw plot looks like vertical lines instead of the usual signal

You need change the scalings parameter. In this case, you need increase it.


okay, there is a slightly newer paper by Ivo Mateljan titled Signal selection for the room acoustics measurement.  i ended up reviewing the paper for the 1999 IEEE Mohonk conference, and when it turned up that Ivo was not going to make it to the conference to present the paper, they ask a reviewer to sorta stumble through the issues of the paper.  i ended up talking about what's wrong with MLS more than about Ivo's technique.

if you cannot get a copy of Ivo's paper (it costs money from IEEE), email me and i will send a pdf to you.  in that paper he has a section titled The Low-Crest Multisine which i paste in the text below.  perhaps you can modify it to fit with your selected frequencies.  i dunno.

i still haven't figured out SE's formatting to make the code look like code.  someone else can edit this answer to do that.

Various algorithms were proposed for lowering the multisine crest factor. The following algorithm is slight modification of algorithms proposed in [9] and [10]. It has fast convergence and resulting multisine has good phase randomisation property.

Algorithm: Low-crest multisine (LCMS) generation.

Generate uniform random sequence u of length N.
Transform the sequence to random sign sequence s:

    for i = 1 : N

        if (u(i) > 0.5)

            s(i) = 1

        else

            s(i) = -1

    end

Generate the multisine using the phase values of DFT(s).
If the number of iterations is less then 10, clip the time series to 85% of the maximum value, else clip to 90% of the maximum value.
Apply DFT to the clipped time series.
Generate the new multisine using the phase values obtained in the step 5.
Calculate the crest factor. If it is higher than objective value go to the step 4.

Although there is no proof for the convergence of this
algorithm, it has never failed in a more than thousand runs with
different starting seed factor. The number of iterations grows
exponentially when decreasing crest factor. For a multisine
sequence of length 32768 it needs only 4-5 iterations to halve the
crest factor from 4 to 2, but to get the crest factor < 1.5 it needs
several hundreds or more than thousand iterations.

Phase randomisation properties of this signal can be
estimated by testing how it distributes nonlinear distortion in the
estimated PIR. The following numerical experiment is done.
The input sequence of length 4096 drives 9th order bandpass
Butheworth filter (0.1-0.3 fs). The output signal is distorted with
2nd and 3rd order nonlinearities, which are equivalent to 1% of 2nd
harmonic distortion and 0.7% of the 3rd harmonic distortion.

Impulse response is estimated using crosscorrelation and scaled
to maximum value 1. Tables I and II show the rms and peak level
of the impulse response distortion (in percentage of the maximum
value). First 50 values of the impulse response are excluded from distortion calculation to better estimate peak/rms tail distortion. Also, the crest factor of the generated discrete sequence and a crest factor at the filter output are shown.

These results show that low-crest multisine generates the same level and a distribution of distortions as a true random phase multisine, regardless the crest factor value. MLS sequence generates slightly lower distortion rms level (1-1.5dB), but peak/rms ratio is much higher, especially with 2nd order distortion (10dB higher than with multisine excitation). This high peak/rms ratio disqualifies MLS as a reliable excitation signal for echo detection.

For all type of an excitation signal the crest factor at the filter output is almost the same. This is a general conclusion for narrow band systems. The benefit of using signals with lower crest factor is significant only for wideband systems, and when it is necessary to use all of the available measurement dynamic range.

RPMS is the only signal, which do not exhibit a significant crest factor transformation. This fact, constant amplitude and true random phase property fully qualify it as a periodic white noise signal.

[9] Schoukens, J., Pinelton, R., Ven der Ouderaa, E., and Renneboog, J., "Survey of Excitation Signals for FFT Based Signal Analysers", IEEE Trans. Instrumentation and Measurement, vol. 37, September 1988.
[10] Ven der Ouderaa, E., Schoukens, J., and Renneboog, J., "Peak Factor Minimisation of Input and Output Signals of Linear Systems", IEEE Trans. Instrumentation and Measurement, vol. 37, June 1988.


Failed to implement Goertzel algorithm in Python

Just had a look, and from the wikipedia definition of the Goertzel algorithm, the frequency in the cosine weight should be a normalized frequency (as for the DFT, by the way). If you modify your code as below, you should get the right output (note also that your computation of the power led to negative powers -sic-, removing a redundant factor 2 solved that issue as well).

The modified code:

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size) # JLD: in Hz
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    f_step_normalized = 1./window_size # JLD: step in normalized frequency
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step_normalized # JLD: here you need the normalized frequency
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - w_real * d1 * d2) # removed factor 2: already in w_real!
        )
        freqs.append(f * sample_rate)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    ##t = np.linspace(0, 1, 44100)
    # JLD: if you do this, the sampling rate is not exactly 44100 anymore:
    #    linspace includes the boundaries 0 and 1, and there are therefore
    #    44100 samples for a bit more than 1 second...
    sample_rate = 44100 # in Hz
    t = np.arange(sample_rate) / np.double(sample_rate) # in seconds
    # JLD: with 1000 samples, a sine at 441Hz leads to a DFT with only one
    # non-nul component:
    sine_wave = np.sin(2*np.pi*441*t)[:1000]
    freqs, results = goertzel(sine_wave, sample_rate, 0, 22049)
    print np.array(results)
    pylab.figure()
    pylab.clf()
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               marker='x',
               label='computed power')
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               linestyle='--',
               marker='o',
               label='returned power')
    pylab.xlim([0,1000])
    pylab.ylabel('Power')
    pylab.xlabel('Frequency (Hz)')
    pylab.legend()
    pylab.show()

I did not thoroughly test that code, and I can't guarantee it to be bug-free, but for your simple example, it seems to work well. I obtain the following figure:

Hope that helps!

Regards!

Jean-Louis!


The problem lies in step 1. Since you need the color information, you should not convert it to grayscale. Convert the image to HSV color space and work on the HUE(H) to filter the objects with color you require.


Matlab and R functions for fft(): N-point FFT and dimension issue

R allows you to loop over rows or columns through the use of the apply family of functions. Note that apply is not necessarily faster than looping but more in line with the functional programming nature of R. There are other apply functions such as lapply, vapply, sapply etc but for your purposes, apply should suffice.  Also, note that you may have to transpose the result using t(result). 1 means go over rows and 2 means go over columns. So, if X is the matrix,

apply(X, 1, fft) does the fft of the rows.

apply(X, 2, fft) does the fft of the columns.

Note that the use of .row or .col is my convention. I could have used x, y , or z because that's just a dummy variable from apply's point of view.

Also, you can type ? apply in R in order to obtain the help for apply. Sometimes they can be kind of obscure but atleast they tell you what arguments are available.


Finite length truncated exponential sequence $\mathcal Z$-transform, zeros over circle explanation

The $z$-transform is

$$\frac{1}{z^{N-1}}\cdot\frac{z^N-a^{N}}{z-a}$$
The zeros will therefore be given by
$$\begin{align}
z^N&=a^{N}\cdot1\\
&=\left(a^{N}\right)e^{j2k\pi}\\
&=\left(ae^{\frac{j2k\pi}{N}}\right)^{N},\,k=1,2,\cdots,N-1
\end{align}$$
Hence,
$$\boxed{z_k=ae^\frac{j2k\pi}{N},\,k=1,2,\cdots,N-1}$$


Why NRZ-L and NRZ-I line coding have DC component?

Assuming bi-polar encoding, and assuming that the number of transmitted 1's is equal to the number of 0's, then over very long time intervals the positive pulses "cancel out" the negative pulses and the average voltage is zero.

However, over smaller time intervals, there's likely to be more 1's than 0's, or vice versa. This is what introduces a DC component to the signal.

Intuitively, you may think about it this way: if you filter the line-encoded signal with a very narrow low-pass filter, you'll see a non-zero, varying output. The DC value is the expected value of the filter's output.

In order to reduce or eliminate DC you need a line encoding that forces the number of positive and negative pulses to be the same, even over very short time spans. One example is the Alternate Mark Inversion (AMI) encoding.


Applying Window Function to Low Frequency Signals

A 1s Hann window will have (in the frequency domain) a main lobe width of 4Hz. This will severely mask your signal components.

Even a rectangular window will be ~2Hz main lobe.

You should use a greater windowing time depending on the accuracy you need for frequency detection.


What is the optimal Wiener-filter for this purpose?

The noise which you defined will only shift your signal and only has zero frequency. i think you intended to define the noise's power. Also in addition to the noise power you need to know your signal's power spectrum to obtain Wiener filter.

Wiener filter will be (1/H)*(S/(S+N)). in this equation H is frequency response of your sensor, S is your signal's power spectrum and N is your noise power spectrum.

to obtain H you have to take derivative of your sensor unit step response to obtain unit impulse response, then take Fourier transform of unit impulse response.


Manual Extraction of Filter Phase Response with MATLAB using sweep sine

You missing one, very important step. In order to get the impulse response of your filter while using the sweep sine, you must get the output of your system $h[n]$ to excitation with sweep-sine signal $x[n]$:

$y[n]=h[n]\star x[n]$

In your case that's the filtering with filter function. Now in order to get the impulse response you must convolve your signal with the inverse filter $f[n]$. Generally that's the time-inverted sweep signal (although special kind of pre-conditioning, amplitude scaling, etc. might be done as well):

$h[n]=y[n]\star f[n]$

Now you have your impulse response that you can use for further analysis. Please take a look at this answer tackling very similar problem.

For more theory please see:

  A. Farina - Advancements in Impulse Response Measurements by Sine
  Sweeps

  A. Farina - Simultaneous Measurement of Impulse Response and Distortion with a
  Swept-Sine
  Technique


RMS in frequency domain using Plancherel theorem

The array $[\sqrt{x_{1}^2 + y_{1}^2 + z_{1}^2}, \sqrt{x_{2}^2 + y_{2}^2 + z_{2}^2},...,\sqrt{x_{3000}^2 + y_{3000}^2 + z_{3000}^2}]$ is a vector of the Euclidean lengths of 3000 3D-space vectors.
The designation RMS (root mean square) is usually used for the square root of the arithmetic mean of the squares of serial values of the same kind. It is not used to designate the value $\sqrt{x^2 + y^2 + z^2}$, which is called the space vector length. The values used for RMS calculation can be received from a series of measurements of the same observable. The components of a 3D vector are not three instances of measurement of an observable; in fact, these vector components are viewed as three observables: the x component, the y component, and the z component.
You have 3000 samples from a series of measurements of acceleration, each instance of measurement produces three vector components of a 3D acceleration vector. When calculating RMS, you iterate over individual items which you select from a vector of measurement samples -- selecting one sample of a total of 3000 samples, squaring the selected item value, summing up these squares in the course of iteration -- then you divide the total sum by the number of samples (3000 here), and finally calculate the square root. You calculate three RMS values: an RMS for the x component of the selected acceleration vector, an RMS for the y component of the same vector, and an RMS for the z component of the same selected vector, like this:
$$
RMS_x = \sqrt{{1\over N}\sum_{n=0}^{N-1}x_{n}^2}; \,
RMS_y = \sqrt{{1\over N}\sum_{n=0}^{N-1}y_{n}^2}; \,
RMS_z = \sqrt{{1\over N}\sum_{n=0}^{N-1}z_{n}^2}
$$
These $RMS_{x,y,z}$ components are components of a root-mean-squared acceleration vector, $\overrightarrow {RMS} = [RMS_x, RMS_y, RMS_z]$
FFT and the Plancherel theorem may or may not be required for processing of your data, dependent on what processing you are going to do on the measurement data.
Notice the usage of the word vector for two different entities: the vector as an element of generalized vector space (as a set of subsequent samples from a waveform) versus the Euclidean vector (a directed segment). See https://en.wikipedia.org/wiki/Vector_(mathematics_and_physics).
It may happen that you need to compute the RMS of acceleration vector lengths, irrespective to acceleration direction in space. If this is the case, the formula is
$$
RMS_{scalar} = \sqrt{{1\over N}\sum_{n=0}^{N-1}(x_{n}^2 + y_{n}^2 + z_{n}^2)}
$$
subscript scalar means that, unlike $\overrightarrow {RMS}$, $RMS_{scalar}$ it is not a vector. Notice also that $RMS_{scalar} = \sqrt{RMS_x^2 + RMS_y^2 + RMS_z^2}$


How to find the output signal of a filter using state space matrices?

We don't know anything about your system, but general rule is:

$$\begin{array}
&y[n] &= \mathbf{C}x[n] + \mathbf{D}u[n]\\
x[n+1] &=\mathbf{A}x[n] + \mathbf{B}u[n] \\
\end{array}$$

Where $u[n]$ is the input signal value, $x[n]$ is state vector value at a given point in time, and $y[n]$ is obviously the output.

What you need to do is to:

Initialize the initial state $x[0]$ to some arbitrary value (usually vector of zeros).
Run the calculation by consequently calculating values of $y[n]$ and $x[n+1]$ to be used in the next step.

Speaking Matlab'ish (not tested, cause I am in travel, but given all informations its extremely easy to figure out what is going on):

b = [0 2 3];
a = [1 0.4 1];
[A,B,C,D] = tf2ss(b,a);

x_0 = zeros(size(A,1), 1); % Initial state
N = 9;
u = randn(1, N); % Some input signal
y = zeros(N, 1); % Output signal

x = u_0;
for n=1:N-1
  y(n) = C*x + D*u(n); % System output
  x = A*x + B*u(n);    % Transition
end


I agree that one of the easiest ways to compute the Hilbert transform in this case is to use the analytic signal. This is most easily obtained via the Fourier transform. Note that the Fourier transform of

$$f(t)=\frac{\sin(at)}{at}\tag{1}$$

is

$$F(j\omega)=\frac{\pi}{a}\text{rect}\left(\frac{\omega}{2a}\right)\tag{2}$$

i.e., a rectangular function with value $\pi/a$ in the interval $[-a,a]$, and zero otherwise. The spectrum of the corresponding analytic signal equals the part of $F(j\omega)$ at positive frequencies, scaled by a factor of $2$:

$$F_a(j\omega)=\frac{2\pi}{a}\text{rect}\left(\frac{\omega}{a}-\frac12\right)\tag{3}$$

The inverse Fourier transform of $(3)$ gives the analytic signal:

$$f_a(t)=\frac{1}{2\pi}\frac{2\pi}{a}\int_0^{a}e^{j\omega t}d\omega=\frac{1}{jat}\left(e^{jat}-1\right)=e^{jat/2}\frac{\sin(at/2)}{at/2}\tag{4}$$

The Hilbert transform of $f(t)$ is the imaginary part of $(4)$:

$$\hat{f}(t)=\text{Im}\{f_a(t)\}=\frac{\sin^2(at/2)}{at/2}\tag{5}$$

as stated in your question.


general purpose microprocessor vs DSP(digital signal processor)?

On a basic level, in my opinion, a DSP chip must, as a minimum requirement, be able to optimally perform the convolution summation used to compute the output of an FIR filter.
$$ y[n] = \sum\limits_{k=0}^{K-1} h[k] \, x[n-k] $$
$y[n]$ = FIR filter output
$x[n]$ = FIR filter input
$h[k]$ = impulse response of FIR filter (also FIR coefficients)
To do this in $K$ instructions, one must be able to in one instruction:

Do a multiplication and addition, implemented as a single multiply-accumulate (MAC) instruction.

Fetch both the delayed sample, $x[n-k]$ and the FIR coefficient $h[k]$ for the next MAC.  To do both, this requires at least two data memory spaces: one for coefficient array $h[k]$, one for the signal array $x[n]$.  Both spaces can be accessed in parallel (simultaneously).  This is called a "Harvard Architecture".

Be able to saturate the results to defined maximum positive or negative values if the result exceeds that range.

Besides being able to saturate to the "rails", the Arithmetic Logic Unit (ALU) in the DSP should be able to round to the nearest LSB of the high-order bits (that will be output or saved to memory) based on the low-order bits of the accumulator.

To address the delayed samples in memory, the ability to do modulo arithmetic on the pointer to $x[n-k]$ is necessary.  This implements a "circular queue" or "circular buffer" or "first-in-first-out" (FIFO) buffer which, in the audio DSP world, is just called a "delay line".  The DSP programmer should not have to worry about when a pointer increments beyond the boundary of a circular queue requiring the pointer to be adjusted or wrapped around.  The DSP chip should make that modulo adjustment to the pointer automatically.

In addition, for the FFT, the ability to do bit reversal or bit-reversed addressing is necessary along with some other handy instructions that make for efficient coding of FFT butterflies.

Now, in order to be a "DSP chip", it really should be able to do all of these operations without explicit instructions in the inner loops.  They build it into the hardware of the ALU the means to perform those operations without explicit effort from the programmer (nor extra machine instruction cycles).  Nothing prevents you from using a normal CPU to do all of these operations, even in a "high-level" language like C/C++ .  And the general purpose CPU chip might be running at 5 GHz while the DSP chip is running at 500 MHz, so the 10x instruction speed might more than make up for the cost of all those extra instructions that the algorithm programmer must code.


I know of no "standard" way of incorporating the output of a set of Kalman filters into another. The reason is that each application will have a different "signal model" (way of combining the multiple outputs into different measurements).

For example, suppose you have two Kalman filters that give you state estimates:
$$
\mathbf{x}_k^1 = \mbox{state of first Kalman filter}\\
\mathbf{x}_k^2 = \mbox{state of second Kalman filter}\\
$$
then you'll need to decide the parameters in your new model:
$$
\mathbf{x}^{\tt tot}_k = [ {\mathbf{x}^1_k}^T {\mathbf{x}^2_k}^T ]^T\\
\mathbf{x}^{\tt tot}_{k+1} = \mathbf{A} \mathbf{x}^{\tt tot}_k + v_k\\
\mathbf{z}^{\tt tot}_k = \mathbf{C} \mathbf{x}^{\tt tot}_k + w_k
$$
where $\mathbf{A}$ is the state update matrix, and $\mathbf{C}$ is the output matrix, $v_k$ is the process noise driving the interactions between the different states and $w_k$ is the measurement noise.


Human eye perception is a complicated thing, the human perception tend to fail in different situations even in large changes of gray tones, also is capable of note minimal changes in the gray levels but all depends on the shape of the figure rather than its histogram.

In some cases you can enhance an image by applying histogram equalization techniques but this is only to emphasize some gray tones in order to observe more details in an image. This is related to gray images.


What Is the Meaning of the Term "Frequency Response'' When We Design a Pitch Detector

Taking a periodic function with many harmonics, and then filtering it, is a common way of generating (or modeling) periodic signals.

It is how you can model human speech, where S is the sound produced by vocal cords, and H is the frequency response of the cavities formed by larynx, mouth, lungs, etc. And this is how many vocoders analyze and synthesize signals.

Imagine the signal S to be a periodic signal with all harmonics; ideally, S may be an impulse train, so that you have harmonics present at all integer multiples of the fundamental frequency.

Then, H is a filter which will affect each harmonic independently, in order to produce (or model) the desired periodic signal.

Since the variable $k$ corresponds to frequency, we can assume that they are actually DFT coefficients. Applying a filter H to a signal S is a multiplication in frequency (or a convolution in time).


Correlation of independent random processes

No. Quoting Wikipedia's article Independence (probability theory):

  If $X$ and $Y$ are independent random variables,
  then the expectation operator $\operatorname{E}$ has the property

  $$\operatorname{E}[X Y] = \operatorname{E}[X]\operatorname{E}[Y].$$

Consider your $X(t_1)$ and $Y(t_2)$ as $X$ and $Y$ in this answer. If both $\operatorname{E}[X] \ne 0$ and $\operatorname{E}[X] \ne 0,$ then the product  $\operatorname{E}[X]\operatorname{E}[Y]$ will be non-zero.


recognize unwanted sounds, noises, and deleting them algorithm

Provided that the noise is uniform across your sound and that you have regions where only noise is present, you can make use of a technique called spectral subtraction. The whole process consist of 4 steps.

  Find a noise-only region and apply an FFT on it in order to obtain the noise profile (noise spectrum). The longer the noise region, the better the noise profile (resolution) you'll get (use longer FFT lengths)
  Apply FFT on your entire sound. This way, you'll get the spectrum for both the signal (people talking) as well as the noise.
  Subtract the noise profile from the overall sound spectrum (2-1). Ideally, you'll end up with the signal-only spectrum (of course, some
  noise will still remain).
  Perform an inverse FFT on the signal spectrum obtained in step #3 to get a time-domain signal (waveform) .

If the filtered signal still has some unwanted artifacts, you may try to interpolate the signal spectrum prior to performing the IFFT (step 4)

Note that you don't have to break up your signal into smaller chunks (this is not STFT)


Group Delay for Hilbert Transformer and Resulting Dispersion

A real-valued system that doesn't distort the shape of the input signal must have the following input-output relation:

$$y(t)=Ax(t-t_0)\tag{1}$$

with arbitrary real-valued constants $A>0$ and $t_0$.

In the frequency domain, Eq. $(1)$ corresponds to

$$Y(\omega)=Ae^{-j\omega t_0}X(\omega)\tag{2}$$

Consequently, the corresponding system is an LTI system with frequency response

$$H(\omega)=Ae^{-j\omega t_0}\tag{3}$$

i.e., the magnitude response must be constant and the phase response must be linear:

$$\phi(\omega)=-\omega t_0\tag{4}$$

Here, linear is meant in the narrow sense of the word, i.e., it must not have an additive constant term (hence, it must not be affine). Note that for real-valued systems, the phase must be  an odd function of frequency, so the sign of a possible additive constant in the phase would need to change depending on the sign of the frequency. Consequently, if there's such a constant phase term (as is the case for a Hilbert transformer) then the phase has a discontinuity at $\omega=0$, which implies that the group delay is undefined / infinite at $\omega=0$.

The phase response of an ideal Hilbert transformer with a possibly non-zero delay $t_0$ does not satisfy $(4)$ because it is given by

$$\phi(\omega)=-\frac{\pi}{2}\textrm{sgn}(\omega)-\omega t_0\tag{5}$$

As mentioned above, the phase of an ideal Hilbert transform has a discontinuity at $\omega=0$.