Confusion with regards to STFT and MFCCs

I am going to break down this answer into two main parts, one for a step-by-step computation of the STFT, and the other in regards to how to compute the MFCC's, once you have the STFT.

STFT: Let us say you wanted to compute the STFT of your signal $x[n]$, of total length $N=100$, and block size $M = 10$, with $50$% overlap, and using a hamming window. I will give you a step by step, (and redundant) flow, to get you started:

Take the first $M=10$ samples out, that is, samples $x[0], x[1], x[2] ...  x[9]$.
Window them, by multiplying them, element by element, with your window coefficients.
Take the DFT of this result.
Take the next $M=10$ samples out, that is, samples $x[5], x[6], x[7], ... x[14]$.
Window them, by multiplying them, element by element, with your window coefficients.
Take the DFT of this result.

etc, etc. You get the idea. Let us say that you took a (same size) DFT in the above example. Then, you will end up with a $10$ x $19$ complex matrix. If you then simply take the absolute magnitude of the above matrix, you will end up with a real matrix, corresponding to the absolute magnitude of your complex matrix. This is called the Short-Time-Fourier-Transform, or STFT. We will use the STFT as the input into computation of the MFCCs.

Before that, let us get some formulas and terminology straight: Let us call our STFT matrix $S$. This matrix is composed of $B=19$ blocks, or columns, and $K=10$ frequency bands, or rows. To refer to a time-frequency point in $S$, we will say $S[k,b]$, where $b \in 1, 2, ... B$, temporal blocks, (columns) and $k \in 1, 2, ... K$ sub-bands, (rows).

MFCC:

Once you have the STFT computed, you can go ahead and use that as a stepping stone for computing the MFCC's. Regarding your question, you seem to be confusing the time domain windowing done on each block for the STFT, and the frequency domain windowing done in the MFCC. They are completely separate and unrelated.

Let us say that you computed your STFT using a hamming window, and you have a $K=500$ x $B=1000$ real matrix. Thus using our prior terminology, you have a real matrix $S_{K \text x B}$.  Thus, you have $500$ sub-bands, (frequencies), and $1000$ blocks, (time).

Let $Y_{b}[k] = |S[k,b]|^2$. That is, we are letting $Y_{b}[k]$ represent the power-spectral density of a particular column/block of the STFT matrix $S$. Obviously, each $Y_{b}[k]$ is a $500$ x $1$ column vector.

Clearly, we have $1000$ $Y$'s. That is, we have $Y_{1}[k], Y_{2}[k], Y_{3}[k], ... Y_{1000}[k]$.

Now, you want to compute the MFCC of each column. Take the first column. This is $Y_{1}[k]$. Forget about all the other columns, and forget about the hamming window. The first column is now used in your MFCC algorithm. The first step in MFCC was getting a PSD column, which you have, because we have $Y_{1}[k]$.

The second step is the computation of, say, 10 power estimates. This is computed by windowing this column, by 10 overlapping windows, as per the mel-scale, like so:

Note how each (triangular) window is mostly zeros. So now you literally do a dot product of each window, with $Y_{1}[k]$. You will get $10$ numbers, for this column. Now take the log of those $10$ numbers. Now take the DCT of this result. Congratulations! You just computed the MFCC for this column. You just computed the MFCCs for $Y_{1}[k]$.

Now, go back, and repeat for all the other $Y_{2}[k]$. Then repeat for $Y_{3}[k]$, etc.


Similarity or Relation between Walsh Hadamard Transform and Slant Transform

On a 2-sample vector, or $2\times2$ pixel images, Haar, Walsh-Hadamard, or Slant are exactly the same. They are all separable, orthogonal transformations with fast $N\log N$ algorithms. According to Slant Transform Image Coding, 1974, it is superior in appearance to the above, only inferior to Karhunen-Loève,  which has no fast enough algorithms.


creating a seamless signal / loop using interpolation

An alternate method is to create an array of two copies of your input signal

y = [ x , x ]         % MATLAB
y = hstack(( x , x )) # NUMPY

Run $x[k]$ through a LPF.  Then slice $N$ samples from the "middle" of $y[k]$ where $N$ is the length of $x[k]$.

xsmooth = y( o : o + N - 1) % MATLAB
xsmooth = y[ o : o + N ]    # NUMPY

Where $o$ is an appropriately chosen offset.  Be mindful of filter "tails" if you carry out the LPF through conventional convolution.


Energy normalization across wavelet subbands

The report under reference is A Really Friendly Guide to Wavelets by Clemens Valens, 1999.
The title of the question might be misleading. The energy normalization concerns the wavelets themselves, required to have the same energy, whatever the scaling and the shift. This is useful for theoretical and practical reasons (algorithmic/computational stability,  and several statistical properties).
If function $\psi(t)$ is of finite energy $\|\psi\|_2^2$ (or square integrable), then by a change of variable $t\mapsto st$, you will find that  $\psi\left(\frac{t-\tau}{s}\right)$ has finite energy as well, and that $\|\psi\left(\frac{t-\tau}{s}\right)\|_2^2 = |s|\|\psi\|_2^2$. For the wavelets to have the same energy, you should divide by $1/\sqrt{|s|}$.
An  interpretation is the following. The energy of function is homogeneous to a rectangle with a time span horizontally, and a square in amplitude vertically.
You want to keep the area of the rectangles constant for each wavelet. So if you dilate the time by $s>0$, you should contract the square in amplitude by $1/s$, then the amplitude by $1/\sqrt{s}$.
For question 2, the spectrum will be stretched/shifted accordingly, if you compute the Fourier transform of $\psi\left(\frac{t-\tau}{s}\right)$. The scaling property of Fourier is mentioned in Clemens Valens tutorial at Equation 14:
$$\mathcal{F}\left(f(at)\right) = \frac{1}{|a|}F\left(\frac{\omega}{a}\right)\,.$$


Ideal Reconstruction of Upsampled Signal

There are some mistakes in the question and your description.

When we sample a signal of frequency $f_m$ with a sampling rate $f_s$, the sampled signal contains the frequencies $f=f_m \pm nf_s$, where $n \in \mathbb{Z} $.
Here the frequencies available in the sampled signal are calculated as,

$f_m=14100 Hz$
 and $f_s=400 Hz$.
we know that,

$ f=f_m \pm nf_s $, where $n = 0, \pm 1, \pm 2, \pm 3, ... $

ie., $ f= ..., -1500, -1100, -700, -300, 100, 500, 900, 1300, 1700, ... $

When we up sample the signal by 3, the frequencies will become $f/3$.

So here the new frequencies are, $ f_u = ..., -633.33, -500, -366.67, -233.33, -100, 33.33, 166.67, 300, 433.33, 566.67, ... $

After upsampling the signal is reconstructed by using an LPF of cutoff frequency 500.
When we apply LPF to the above signal only the frequencies whose absolute value is below 500 remains and other frequencies will loss.

there for after applying LPF the frequency components available in the signal are,  $ -500, -366.67, -233.33, -100, 33.33, 166.67, 300, 433.33 $.

I hope this solves your problem.


JPEG is far simpler.  It divides the image into 8x8 pixel blocks, and processes each using a Discrete Cosine Transform.  The results are quantised and then encoded.  The quality is fixed by the encoder.

JPEG2000 uses a 2D wavelet function, the output of which is four "images", each a quarter the size of the original.  One of those is actually an image, while the others are high-frequency components that can be added to it to re-construct the full-resolution image.  The wavelet process may be repeated multiple times.  The result is a tiny image, and a series of high frequency components that may be combined with it.  Each resulting component is quantised and encoded.

JPEG is fine for high quality and modest compression, which is why it is still very widely used.

JPEG2000 offers several advantages:

To achieve very high compression, it is possible to throw away or heavily quantise the high frequency components.  This gives a poor quality - but usable - image where JPEG would fail completely.
Images can be re-constructed progressively at ever improving quality.  This can be either in terms of increasing resolution or bit depth as required.
It supports the JPIP protocol for progressively transmitting images to a client.  The client may retrieve low-resolution thumbnails and then just the parts of the image they want at better resolutions.


The question is not so clear. I am not sure what other types of examples you consider. In this particular case, if you didn't get the hint notice that:

$$\begin{align}
c(1+\alpha D)(1+\alpha D^{-1}) &=  (c+\alpha^2c)+\alpha cD^{-1}+\alpha cD \\&=0.371+0.09D^{-1} +.09D
\end{align}$$

which gives you
$$\left\{\begin{matrix}
c(1+\alpha^2) = 0.371\\
\alpha c=0.09
\end{matrix}\right.
\Rightarrow \frac{1+\alpha^2}{\alpha}=\frac{0.37}{0.09}$$
and you will actually have two solutions

$$\left\{\begin{matrix}
\alpha=3.8634, c=    0.0233\\
\alpha=0.2588, c=  0.3477
\end{matrix}\right.
$$

See if the second answer is also valid...


Inverse $z$-transform of a transfer function in MATLAB

The function iztrans only supports symbolic expressions. This should do the trick:

[num,den] = tfdata(H);
syms z
H_sym = poly2sym(cell2mat(num),z)/poly2sym(cell2mat(den),z);
h_inv = iztrans(H_sym);


Have you seen these kind of DFT artifacts?

I think your sine have some sample(s) missing or is poorly sampled. Do you have a way to check that? You said that you use 65.1k to sample the sine, but what is its frequency?


How to do edge detection using Curvelet transform?

If you decide to stick to curvelets, you can find a precise description of an edge dectection algorithm is  Edge detection in microscopy images using curvelets, Tobias Gebäck and Petros Koumoutsakos, BMC Bioinformatics, 2009. The pdf is there. And you are very lucky: you can get, for free, an  implementation, CurveletUtils. It provides a Matlab source code for a GUI implementing the edge detection method. I hope you can get started with that.

If you care of other multiscale transforms providing different oriented  edges at different scales, you can take a look at a review paper  on 2D directional wavelets (including curvelets, shearlets, contourlets, etc.). It includes morphological and non-linear wavelets.

By the way, I do not understand why you call Canny and Sobel "morphological edge detectors".


Clarification regarding signal after pulse shaping

I assume in the following that your digital signal is a sequence of width-1 spikes or, equivalently, a sequence of delta pulses.

Your half-sine pulse shaping function gives you a lagged peak for every spike. Consider just your first digital signal spike at $t = 0$. If that were all you have, then the shaped output signal would be
$$g(t) := f(t) \star p(t) = \sin(\pi t / 2T_c),$$
which has a peak at $t = T_c$. That is, the original signal had a peak at $t = 0$ and the shaped signal has a peak at $T_c$. Is this the behavior you intended?

If you only had your first two peaks, then the shaped signal would be the sum
$$ g(t) = \sin(\pi t / 2T_c) + \sin(\pi t / 2T_c - T_c).$$ We can use standard trig identities to rewrite this. Set $\theta = \pi t / 2T_c$ and $\phi = \theta - T_c$, we have $\theta + \phi = 2\theta - T_c$ and $\theta - \phi = T_c$, so we can write
$$\sin(\theta) + \sin(\phi) = 2\sin(\frac{\theta + \phi}{2}) \cos(\frac{\theta - \phi}{2}) = 2\sin(\theta - T_c/2)\cos(T_c/2).$$
This function has only one peak, at $\theta = \pi/2 + T_c/2$. In other words, two peaks in the digital signal correspond to only one peak in the shaped signal. I think this is the behavior which surprised you.

Generally speaking, and without doing the math, I'd expect $N$ (edit) consecutive same-sign peaks in the digital signal to give rise to $N-1$ peaks in the shaped signal.

edit: Without doing the math, I'm not very sure what I expect out of a positive peak followed by a negative peak. I expect you'll get a positive shaped peak followed by a negative shaped peak.


Why are there comb-like hills in a periodogram?

A single spike (as you call it) appears theoretically only for a infinite-length sinusoid. Since your signal is 100 samples length, it is not infinite. You actually multiplied your infinite signal with a window which has a value of 1 over 100 samples, and 0 elsewhere. Since multiplication in time domain is equivalent to convolution in the frequency domain, your spectrum is a convolution of the single spike and the frequency response of the window (btw it is called rectangular window). This is the function you got.

I suggest you read about windowing: http://en.wikipedia.org/wiki/Window_function


Hide text in audio and got text by phone microphone when the audio is plaid anywhere

I'm not sure about how much experience you have in writting software, but if what you want is to hide some data information within an analog audio signal, it seems for me like the way to go would be to use spread spectrum modulation and encode the data digitally and add it to your audio signal as some background noise.
Spread spectrum is in itself already a form of encryption, so the bytes you add there will be from the get go, very difficult to recover for a hacker unless they: 1-know that are there, and 2-have the key for generating the pseudorandom sequence to de-spread the signal.
But if you want even more security you can choose to packet your data by SSH tcp packets, and have in your phone a program to extract those from the digital signal.

As for point number 2. My suggestion immplies that you add something to the analog signal, or if you will to the PCM raw waveform, which is the same for all practical purposes. So, taking into account that the job of an encoder/compressor like mp3 is to compres the signal, but to eventually recover the original waveform with as less RMS error as possible, we could say that the original noise waveform that you add will still be there, slightly modified, but still there, and hence when you recover the bits by demodulation, and apply suitable Error correction tecniques, I think your data should still be there intact as it was put into the sound.

Please take into account that i just suggested a possible technology, normally used with narrowband RF signals, but maybe you could do some research and check if it will work in your scenario.

Also Another technique that could be used, if you want to send Text (Like they do in the video), is to generate based on your text a signal that will look like the ASCII characters in the time-frequency plane. Then to recover you just perform time-frequency analisys and use any image recognition algorithm.
This just came to mind, so it probably sounds very wierd, but may be another way to add another signal to an audio file.

Hope this helps.


What is the no overload region of a quantizer?

It is stated in chapter 5 (Scalar Quantization I: Structure and Performance) of the book "Vector Quantization and Signal Compression" by Gersho-Gray that

  Associated with every $N$ point quantizer is a partition of the real line
  $\mathcal R$ into $N$ cells or atoms $R_i$, for $i = 1,2,\cdots, N$. The $i$th cell is given by $\mathcal R_i=\{x \in \mathcal R: Q(x) = y_i\} \equiv Q^{-1}(y_i)$, the inverse image of $y_i$ under $Q$. It follows from this definition that $\bigcup_i \mathcal R_i = \mathcal R$ and $R_i\bigcap R_j = \emptyset$ for $i \neq j$. A cell that is unbounded is called an overload cell. Each bounded cell is called a granular cell. Together all of the overload (granular) cells are called the overload region (granular region).

So in other words, the union of bounded regions (i.e. granular region) where saturation does not occur is referred to as the "no overload" region and is identified in the following picture


On the signal recovery in OFDM after estimation

There is a mistake int the conv function which you are using in your code. In ofdm, the channel must be convoluted with every symbol.
the link provided in the above comments are ok, but you need to modify them according to your parameters.  so replace the command of y = conv(x_ifft_p2s,h,'same'); by below command:

for jj = 1 : n_ofdm_sym
         y(jj,:) = conv(h(jj,:),x_cp(jj,:));
end

then add the noise into y.

NP: be careful if the n_ofdm_sym  is more than one, you need to build your channel $h$ as a matrix of size $nofdmsym$ x $nsubc$ .

Good luck


Is this code for doing a moving average correct?

Yes, this code is correct for implementing a moving average filter.

Nevertheless I do recommend to use the in-built smooth function in MATLAB.

y = smooth(x, 100, 'moving')


Licence plate detection with different backgrounds

Here are a few ideas:

First convert from color to grayscale. It looks like you have fairly good contrast already. There are various methods to perform this conversion; choose the simplest at first: gray = (red + green + blue)/3. Quite often you don't need anything better than that. For some applications, using just the green color plane is sufficient.
If you want the algorithm to operate on black text, but if some images have white or light text on a dark background, reverse the polarity of the grayscale image by swapping grayscale 255 for 0, 254 for 1, and so on until all images have black text on a light-color background. You'll need an algorithm to detect when this is necessary, OR you can run an algorithm on the original grayscale image and on the negative (reversed brightness) image and choose the better result. Keep in mind that the "close" morphological operation has a different meaning for light-on-dark text versus dark-on-light text.
Investigate some adaptive thresholding methods. Using a global Otsu threshold can be problematic--if possible, you should first have some sort of method to detect regions that are likely to be license plates, and then threshold and binarize those regions only. Thresholding an entire image will likely cause trouble.
Look into the Stroke Width Transform and improvements upon it. The SWT or equivalent can help either with segmentation, or with identifying likely license plate regions, or possibly even with providing edges useful for identification.
Since in your region there are only a few color combinations, you might be able to create some method that would could a template for every possible color combination and then try to find a match at each of several image scales. The license plate might be 400 x 250 pixels, or 100 x 65 pixels, or whatever.
For multiscale operations, consider holding an image pyramid in memory. In some cases you can perform an initial match in the smallest image of the pyramid and work your way down. With some modification you could create an image pyramid has a scaling factor of less than 2 from level to level.
Pull the sample images into Photoshop or GIMP and play around with the filters. You may be surprised to find a filter that is helpful.
Whenever possible, work on raw, uncompressed images rather than on JPEGs or other compressed images.
You may not have much of a choice in hardware, but given the fuzzy look of your "LEA 1956" sample, it'd be nice to have more resolution (more pixels per image) and/or sharper optics and/or a raw image instead of a JPEG.

Try to generate the best possible images that you can with pre-processing before your SVM goes to work. Thresholding is often the weak link here. Thresholding and binarization is a common approach, but when you see how noisy a binary image can be you should think very hard whether that is the kind of image you actually want to process (no matter how many people may claim that this is the standard method and hence the best).

Even if you're already starting with a framework and/or algorithms, I recommend doing a little background reading to gain an understanding of both human reading and the many approaches to OCR:

Character Recognition Systems by Cheriet and others is a good reference (though a few years old now) of the various techniques for Optical Character Recognition.

Reading in the Brain by Stanislas Dehaene gives a good high-level overview of human reading, and covers the range of scripts. It's a relatively fast read, and has references to academic works that are useful.


Understanding the frequency scale of a spectrogram

I think what you're saying is that you have a linear, time invariant (LTI) system model for both the propagation medium and the response of the measurement instrument itself. Given that you know the (at least approximate) magnitude response of the instrument, you would like to best estimate what the actual signal at the instrument's input looked like.

This is similar to the process of equalization used in communications systems. Since you have some knowledge of the system (in this case characterized by the impulse response $i(t)$ or the frequency response plot that you showed), you can attempt to undo the effects on the signal that are imposed by the measurement instrument. More generally, this type of operation is known as deconvolution: given an output of a convolution and optionally knowledge of one of the inputs to the convolution, one would like to calculate the other input. In this case, one of the convolution inputs is $i(t)$ and the output is the measurements reported by your instrument.

It appears that you don't know the exact impulse response of the measurement device, only its magnitude response. This isn't enough information in itself to execute the deconvolution, so you'll need to make some assumptions about the phase response of the instrument in order to move forward with that approach. This is a very powerful technique, however; as the Wikipedia article notes, the approach can be used to compensate for the effects of the seismic propagation medium also (since it's also modeled as a convolution with an impulse response), potentially allowing you to gain more insight as to the original seismic event.

Edit: And, as to your questions about how you're interpreting the plot, your understanding seems correct. The system has a highpass characteristic with a flat passband starting at approximately $0.02\ \text{Hz}$, with the magnitude response decreasing by a factor of 10 per decade of frequency below $0.02\ \text{Hz}$. This could be modeled well as a single-pole highpass filter.


What's the Difference Between LMS and Gradient Descent Adaptation?

The LMS algorithm is based on the idea of gradient descent to search for the optimal (minimum error) condition, with a cost function equal to the mean squared error at the filter output. However, it doesn't actually calculate the gradient directly, as that would require knowing:

$$
E(\mathbf{x}[n]e^*[n])
$$

where $\mathbf{x}[n]$ is the input vector and $e[n]$ is the filter output error at time instant $n$. The value of this expectation would typically be something that isn't known ahead of time, so it's instead approximated (in most cases) by the instantaneous value $\mathbf{x}[n]e^*[n]$ instead. This (described in more detail at Wikipedia) is the key feature of the LMS filter structure.

Gradient descent just refers to the method used to hunt for the minimum-cost solution; it doesn't force the use of any particular cost function. An LMS filter is a specialization of gradient descent that uses a mean-squared error cost function and the aforementioned approximation for the gradient at each time step.

Stochastic gradient descent would refer to techniques that don't directly compute the gradient, instead using an approximation to the gradient that has similar stochastic properties (e.g. the same expected value). LMS would be an example of an algorithm that uses a SGD approach, using the approximation I described above.


Frequency response of filter  signal mix

The Z transform is linear. So if
$$ Y(z) = H(z) \cdot X(z) $$ then
$$ Y^{'}(z) = a \cdot X(z) + (1-a) \cdot Y(z) = a \cdot X(z) + (1-a) \cdot H(z) \cdot X(z) = (a + (1-a) \cdot H(z)) \cdot X(z) $$

Hence
$$ H'(z) = a + (1-a) \cdot H(z) $$


Are all exponential functions eigensignals of LTI systems?

Complex exponentials are eigenfunctions of LTI systems because they are eigenfunctions of the convolution operator:
$$\begin{align}e^{j\omega_0t}\star h(t)&=\int_{-\infty}^{\infty}h(\tau)e^{j\omega_0(t-\tau)}d\tau\\&=e^{j\omega_0t}\int_{-\infty}^{\infty}h(\tau)e^{-j\omega_0\tau}d\tau\\&=e^{j\omega_0t}H(j\omega_0)\end{align}\tag{1}$$
where $h(t)$ is the system's impulse response and $H(j\omega_0)$ is its frequency response (evaluated at $\omega=\omega_0$).
Now try to prove in the same way that the convolution of an impulse response $h(t)$ with a real-valued exponential results in a scaled version of that exponential. Then draw your conclusions.


in my opinion, either of these two functions of time can be rightly labeled "sinusoidal":

$$ x(t) = A \sin(\omega t + \theta) $$

$$ y(t) = A e^{j(\omega t + \theta)} $$

where $A$, $\omega$ can be any non-zero real number (i wouldn't call DC a "sinusoid") and $\theta$ can be any real number.


The Lyons book is gold for someone in your situation.  Have patience with the topic of audio DSP -- the more time you spend trying to learn it, the less you will know that you know.  :)


Identifying the magnitude and impulse response from pole zero plot quickly

Some of the cues you want to look for are the following:

A zero on the unit circle causes a zero in the magnitude response at the corresponding frequency ($z=1$ is DC, and $z=-1$ is Nyquist)
For a zero at DC ($z=1$), the sum of all impulse response coefficients must be zero.
A zero near the unit circle causes a dip in the magnitude response at the corresponding frequency. The closer the zero to the circle the more pronounced the dip.
A pole close to the unit circle causes a peak in the magnitude response at the corresponding frequency. The closer to the circle the sharper the peak.
Complex conjugate poles cause oscillation. The closer to the unit circle they are, the slower does the impulse response decay.
If all poles are at the origin, the impulse response has a finite length (FIR).
If there are fewer (finite) zeros than poles then the first few samples of the impulse response are zero. The number of zero coefficients equals the difference between the number of poles and (finite) zeros.

For PZ-plot A, the magnitude response must be zero at DC and it must peak at $\pi/2$ (i.e., at $0.5$ in the magnitude plot). The impulse response (IR) must be oscillating and its values must add up to zero (which you can of course check only approximately). Furthermore, since there are two poles and only one zero, the first IR coefficient must be zero. So we are left with IRs 1 and 5, and with magnitudes b and c. Since we have a PZ-plot which is qualitatively the same as A (namely D), just with a smaller pole radius, we can easily guess that the peakier magnitude and the more slowly decaying IR correspond to A, whereas the others correspond to D. So we get

  A-5-b

  D-1-c

PZ-plot B must be FIR (all poles in the origin), and the first coefficients must be zero (3 poles, only 2 zeros). The magnitude must be symmetrical with respect to $\pi/2$ ($0.5$ in the magnitude plot), and there must be (equal) dips at DC and at Nyquist. Only magnitue plot "a" satisfies these criteria. The only matching finite IR is "4". So we have

  B-4-a

PZ-plot C is also FIR, and the first coefficient is not zero (because the number of poles equals the number of zeros). The magnitude must have a slight dip at $0.5$ and it must be symmetrical (because the PZ-plot is). This gives

  C-2-e

This would leave us with

  E-3-d

but the latter is non-sense, as I've mentioned in my comments. The PZ-plot E corresponds to a finite IR, but that IR must be complex valued. From the plot, let's estimate the zero at $0.8j$, which gives us

$$H(z)=c\cdot (1-0.8jz^{-1})\tag{1}$$

with an arbitrary scaling constant $c\neq 0$. This corresponds to an impulse response

$$h[n]=c(\delta[n]-0.8j\delta[n-1])\tag{2}$$

which is obviously complex-valued. The magnitude response must have a dip at 0.5, but no dip at -0.5 (note that for real-valued impulse responses the magnitude response is symmetrical). IR 3 is just an impulse at zero, so the corresponding PZ-plot must be empty, because there are neither poles nor zeros. IR 3 just shows the real part of the actual impulse response. The magnitude response "d" can be the response corresponding to E, but only for positive frequencies, which are not sufficient to characterize the system. So the PZ-plot E is misleading, to say the least. The same is true for the single impulse shown in 3, because it corresponds to none of the PZ-plots and to none of the magnitude plots shown in the figure.


Tone generation in Matlab (scientific paper)

The section on 12 added sidetones is referring to something different to the section on thirteen notes.

In regards to 13 notes: By diatonic I think the mean a diatonic musical scale, with the lowest note at 148.83Hz.

From this site:  http://www.phy.mtu.edu/~suits/notefreq446.html   148.83Hz is note $D_3$

So I would choose 13 consecutive notes in the D-scale starting with $D_3$. The notes are played according to Figure 1. Each note lasts 8.40s/72 = 7/60 s.

However the notes look equally spaced in frequency in Figure 1. So I think they might mean a chromatic scale with equal temperament. That would make more sense as the 13th note would be one octave above the first note, and each note would be $^{12}\sqrt{2}$ apart. Code then:

low_note = 148.83;
notes = low_note .* nthroot(2,12).^(0:1:12);

signal = [];

% First section
for k = 1:12
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end
% Second section
for k = 13:-1:2    % Reverse
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end
% Third section
for k = 1:12
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end


Autocovariance is defined on random signals.

$$
C_{XX}(t,s) = E\left[(X_t - \mu_t)(X_s - \mu_s)\right] = E[X_t X_s] - \mu_t \mu_s
$$

The expected value only makes sense for random processes. You could treat a deterministic process as a random process where everything except for the one signal you observed has probability zero.


How to choose pilot symbols count in OFDM?

It is a tradeoff. The more pilot subcarriers, the better then channel estimation, but the fewer subcarriers usable for data. Simulations and testing of performance can help in making the choice.

You probably have some numbers wrong. If you have 224 subcarriers, you cannot have 346 of them as pilot carriers.


Does Hamming encoding preserve orthogonality?

I don't think you cannot state generally that $y_1$ and $y_2$ will be orthogonal. I'll try to sketch out my thinking:

Since the Hamming code is a linear code, each parity bit can be represented as a linear combination of the information bits; that is, each bit can be represented as:

$$
p_i = \mathbf{g_i^T}\mathbf{x}
$$

where $\mathbf{x}$ is the information bit vector and $\mathbf{g_i^T}$ is the generator for the $i$-th parity bit (it's one column from the code's generator matrix $\mathbf{G}$); it is a binary vector that represents the parity checks that go into $p_i$. When formulated this way, the parity bit can be thought of as the projection of the information bits onto the respective generator vector.

Assuming that you have a systematic code (all Hamming codes that I've seen before are), the encoded vectors $\mathbf{y_i}$ look like (assuming the code has $k$ information bits and $m$ parity bits):

$$
\mathbf{y_i} = \left[x_{i_1}\ x_{i_2}\ \ldots \ x_{i_k}\ |\ p_{i_1}\ p_{i_2}\ \ldots \ p_{i_m} \right]^T
$$
$$
\mathbf{y_i} = \left[x_{i_1}\ x_{i_2}\ \ldots \ x_{i_k}\ |\ \mathbf{g_1^T}\mathbf{x_i}\ \mathbf{g_2^T}\mathbf{x_i}\ \ldots \ \mathbf{g_m^T}\mathbf{x_i} \right]^T
$$

In order for the vectors $\mathbf{y_1}$ and $\mathbf{y_2}$ to be orthogonal, their dot product must be zero. It has the form:

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^k \left(x_{1_j} + x_{2_j}\right) + \sum_{j=1}^m \left(\mathbf{g_j^T}\mathbf{x_1} + \mathbf{g_j^T}\mathbf{x_2}\right)
$$

Since the information vectors $\mathbf{x_1}$ and $\mathbf{x_2}$ are orthogonal, the first sum above is zero, so we are left with:

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^m \left(\mathbf{g_j^T}\mathbf{x_1} + \mathbf{g_j^T}\mathbf{x_2}\right)
$$

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^m \mathbf{g_j^T}\left(\mathbf{x_1} + \mathbf{x_2}\right)
$$

There is no requirement that this condition be zero for an arbitrary linear block code. I haven't taken the specific structure of the Hamming code into account here, but I suspect that it doesn't happen to meet this criterion either.


Convolution output signal is very large, do I need to scale inputs or outputs?

In practice you would make sure that your filters don't have very large gains, at least not over a large frequency range. The filter you're using in the given example has a maximum gain of almost $26$ dB (at DC), and the gain is no less than $12$ dB for any frequency. I.e., you're applying a flat gain of $12$ dB, and on top of that you boost a relatively wide band around DC by another $14$ dB. This is not what a typical filter for audio signals would do, unless you make sure that the input signal is small enough and the additional gain is really needed.


Problems with IFFT not being symmetrical

Here is an example using Matlab/Octave code. There is no need to conjugate.

x = rand(1000,1);
y = rand(1000,1);
z = conv(x,y,'full');
y2 = ifft( fft(z,1999) ./ fft(x,1999) );
norm(y - y2(1:1000)) / norm(y)

ans =

   1.5654e-13


Relation between Power Spectral Density and RMS of the signal

instantaneous power is:

$$ p_x(t) = |x(t)|^2 $$

mean power is:

$$\begin{align}
 P_x &= \lim_{T \to \infty} \quad \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} p_x(t) \, dt \\
\\
&= \lim_{T \to \infty} \quad \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} |x(t)|^2 \, dt \\
\end{align}$$

the relation between Power Spectral Density to mean power is:

$$ P_x = \int\limits_{-\infty}^{\infty} S_x(f) \, df  $$

the relation between Root Mean Square to mean power is:

$$ P_x = \lVert x(t) \rVert^2 $$


Banding in the derivative of a lock-in amplifier signal

It might not be so straight forward to simply take the derivative of the output from a lock-in amplifier. If you think about how a lock-in works, the final stage of a lock-in amplifier is the low-pass filter. So if you then simply take the sample-to-sample derivative, you will simply be undoing the work of the low-pass filter and getting noise. Instead, you need to decimate your time domain signal, and then take the derivative.


The low data rate comes from the need to transmit the data reliably over a noisy channel with a low power transmitter.

A GPS receiver has to be capable of picking out a signal under conditions that other radio systems would just call "total loss of signal."
In a typical two-way radio (analog police radio and the like) the signal is usually unusable if the received signal strength gets below -120dBm.  Digital signals usually get lost long before you get down that low.

GPS receivers work in ranges where a good signal is -155dBW.  In order to fish the data out of the noise, they use a faster signal (1023 million "chips" per second for the civilian codes) that is then modulated with the actual data.  You modulate the data on the chip signal, and modulate that onto the transmit frequency.  You use correlation on the chip signal to recover it, and demodulate the data from that.

The end result works like a very narrow band filter that only lets through data at a low rate, but very clearly.

The chip modulation scheme also has the advantage that you can transmit a bunch of them on a single carrier frequency, but be able to separate them cleanly on the receive side.  This is known as Code Domain Multiple Access (CDMA) and is the reason that the US GPS system only needs a single transmit/receive frequency despite having some 30 satellites.  Glonass uses 15 carrier frequencies (each used twice) to cover the globe with 30 satellites.  It makes the transmitters and the receivers just that bit more complicated.

I would suggest reading the Wikipedia GPS article.  It includes a lot of details I've glossed over in trying to keep this short.


Is there a difference what measurement units use in covariance matrix

Most certainly yes.  The units of the states are reflected in the state covariance, and state error covariance.

The units of the measurements propagate throughout the filtering and prediction steps.

Should units be appropriately scaled? Numerical accuracy typically benefits from appropriate scaling


Your first problem is memory - even the largest version of that processor doesn't have that much internal SRAM (in video terms - it's a stonking embedded processor compared to many though!) - a VGA frame is 300kB, vs the 100ishkB in the processor.  That means you'll have to process it as it comes in - say in 8 or 16 line chunks, which makes it much more of a "hard-realtime" problem as the deadlines have less slack in them.

It's not clear how you'll be capturing data, but I assume that you'll have some sort of DMA from either the ADC or some digital port, otherwise you'll be spending half your time shuffling data around!

Aiming for lowish compression ratios, you might just want to look at encoding the deltas between pixels along a line or a 3x3 square and then arithmetic code them.

See also Simple, streaming, lossless image compression - although I was asking for explicitly lossless compression, there may be some ideas for you there.

As a datapoint, one of my colleagues implemented a 320x240 8-bit mono JPG compression in a 60MHz 32-bit micro with floating point unit and a small cache.  He used some reference code (ie optimised for readability, not performance) and got 5 fps fairly easily.  Compression ratios were of the order of 10x IIRC. The image capture was done by external hardware bus-mastering the data into the micro's external RAM.


Discretize process noise in Kalman filter

Hi: I've never dealt with discretization in the kalman filter ( my models were already discrete ) so take the following  with a level of uncertainty ( no pun intended ). Also, you didn't show the original equations so I'll refer to them as the observation equation and the system equation. Based on your updating equations,  it seems that you're using the non-lagged formulation where the observation equation and the system equations are updated at the "same" time rather than the observation equation being one lag ahead.

So, assuming the above is correct ( without the original equations, it's a little hard to tell ),  the reason it's multiplied by $\triangle t$ twice is because there are two updates to the system noise due to discretization of the two equation process.  The first update is due to the $G_k$ matrix multiplying the variables in the observation equation at time $k$.  That's what causes the first $\triangle t$ multiplication. But then, you also have the update of the system equation at time $k$ due to the multiplication of the state variables by $F_{k-1}$. This is why
$\triangle t$ is used a second time.

But, rather than think of the discretization as using multiplcation twice, you can think of it as using multiplication sequentially. It looks as though once is multiplying twice because, in your formulation, the observation equation and system equation are not lagged with respect to each other. But it's really just sequential multiplication due to sequential updating. Note that if you used the other KF formulation where there is lagging, then the discretization process would probably be easier to see. I hope this helps some.

                                                 Mark


HINT:

It can be solved by convolution. Take two rectangles, one that is non-zero in $[-T/2,T]$ and the other being non-zero in $[0,T/2]$. Now you just have to choose the amplitudes right and the convolution of the two signals will look like the one in your question.


Convert Fourier transform of a signal to $\omega$ form?

In your example there is no sampling involved, so you simply have $\omega=2\pi f$, where $f$ is the frequency in Hertz, and $\omega$ is the frequency in radians.

Your first result for $X(f)$ looks correct. However, your result for $X(\omega)$ is wrong. What you need to know to be able to convert an expression with Dirac deltas from $f$ to $\omega$ is the following relation:

$$\delta(af)=\frac{1}{|a|}\delta(f)\tag{1}$$

With $a=2\pi$ and $\omega=2\pi f$ you get from $(1)$

$$2\pi\delta(\omega)=\delta(f)\tag{2}$$

You get the two following basic Fourier transform correspondences, depending on whether you use $f$ or $\omega=2\pi f$ as the frequency domain variable:

$$\begin{align}e^{j2\pi f_0t}&\longleftrightarrow \delta(f-f_0)\\
e^{j\omega_0t}&\longleftrightarrow 2\pi\delta(\omega-\omega_0)\end{align}$$

So your expression for $X(\omega)$ misses a factor of $\pi$, and the sums and differences of the two frequencies in the arguments of the Dirac deltas need to be multiplied by $2\pi$:

$$X(\omega)=\frac{5\pi j}{2}[\delta(\omega-2\pi9000)-\delta(\omega+2\pi9000)+\delta(\omega+2\pi11000)-\delta(\omega-2\pi11000)]$$


Inter aural time difference generation using sine waves

I understand the point.
If you are outputting signal[i] you cant jump to signal[i+n] instantaneously.
You should just fade it to the new point with some window function. Allow some period of fading time.

For example: If I currently output signal[i] and I want to shift it to signal[i+n] in 1 second. I will use half of Hann window to fade it in and out. So I time scaling Hann from pi to 1 second.
and output $signal[i+n]*Hann[i]+signal[i]*(1-Hann[i])$ for 1 second before I switch to signal[i+n] at full scale. This should not generate a click.

You may want to try other window function or different fading period.


Non Equispaced / Non Uniform DFT Bandwidth

Your formula isn't accurate.  Since you aren't trying for speed and internally  consider the interpolation between data points as a step/pulse then the formula should be.

$$
S(\omega)=\sum_{n=0}^{N-1}s_{k}sinc(\frac{\omega\triangle\left(t_{k}\right)}{2})e^{-wt_{k}}
$$
Having said that, what is "adequate"?  This is, more or less, subjective.  The original Nyquest criterion is based upon the result that any (freq < sam/2) (I think there is some small quibbling on that) can be perfectly reconstructed; if it and the sampling lasts long enough.  It says nothing about out of band signals; but the actual calculations will show aliasing.
It behooves you to generate the signals and interference that you envision and run it through your proposed algorithm. An alternative is to treat your system as a filter and try to "deconvolve" it. This is hazardous but can be done.
For a in depth look (with the hazards) read "Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements"  SIAM.
Now MRI/CAT scanners use various means of reconstructed signals.  All of the mathematicians seem to think that their technical success is an example of extreme luck; although I am sure that the designers might take exception to that.  You might try borrowing some of their mathematical techniques though.
Now in general instrumentation what we concentrate on is "noise temperature" and such.  This can be translated as: how far will your system signal/noise (or signal/distortion) be from the optimal number.  If you can get within 1db of the optimal S/N value then you should probably give up.  Of course this implies that you can calculate the optimal S/N.  In most instruments this can be done and this criterion gives you a stopping point; a point to give up and look for something else to do :)


Non-uniform acquisition to fixed sample rate

[WARNING, RUDE OPINION AHEAD]  Engine people often record data in a speed-invariant fashion using angular sampling. I have been working in that domain for a while, and I bear with you. I have been testing the performance of $0.1$ CA sensors, most of them doing some weird and undocumented interpolation. I also have tested $6$ CA sensors, at the other end of the camshaft.
My opinion so far is that they don't match well, for different reasons,
including the acyclic behavior of an engine:

some angular devices are crappy,

the analog filtering (if any) in the acquisition chain is not speed-dependent,

the speed cannot be assumed constant within a cycle, and even within a stroke.
And the measurements can be quite uneven at both ends of a camshaft, indicating that in between, at each cylinder, depending on the engine organisation, you might have various effects. This is all the more detrimental with noisy knock sensors :  their spec sheets are often imprecise, their calibration is often neglected, etc.

I have concluded that in a standard speed range (1000-4500 RPM), a precision below 0.2 CA is (often) meaningless. And at high speeds, 0.5 CA is doubtful.
So :

if you only have 0.1 CA data, do the simplest interpolation that copes with the processing you need. A simple LS interpolation (linear, splines) could do the work
if you have access to the angular timing for raw files, get rid of the spikes, and use it directly have a non-uniform sampling (but don't forget you only see it from one end of the shaft)
if you have data models, plug them into some super-resolution tool.
-if none of the above, well...

I'd suggest a track anyway, that I did not push to the end (so be careful): try to bound or to "probabilize" uncertainties: speed variation, quantization, filtering, etc. From the given measure, make a fuzzy measurement, ie a simulated multivariate object, that will go through further processing, which can give you a rough approximation of uncertainties pertaining to subsequent processing.


Why Cramér spectral representation and not DTFT for stochastic process

I will introduce some terminology and intuition that will be helpful when reading other references. It will be neither complete nor completely rigorous.

The measures that we first encounter in real analysis assign sizes (non-negative real numbers) to measurable subsets of $\mathbb{R}$; Lebesgue measure is the measure that agrees with the intuition we build in calculus (the measure of the interval $[a,b]$ is $b-a$, etc).

$Z$ is a measure, but it is a stochastic measure†. It does not assign numbers to measurable subsets of $[0,2\pi]$. Rather, it assigns a random variable to each such subset:
\begin{equation}
X_A = \int_{A}dZ(\omega).
\end{equation}
The convergence of the integral on the right-hand side is an issue that I would rather not try to explain (to you or to myself).

In particular, the $Z$ used for WSS processes is an orthogonal stochastic measure. One result is that random variables assigned to non-overlapping sets are independent of one another.

If $A$ is a Lebesgue-measurable set, then $Z(A)$ is a random variable, and the expectation of $\left|Z(A)\right|^2$ is
\begin{equation}
\mathsf{E}[\left|Z(A)\right|^2] = \textrm{Lebesgue measure of $A$}.
\end{equation}
Hence, Lebesgue measure is "under the hood" even if we stick with the notation $dZ(\omega)$.

Just as we can use Lebesgue measure to integrate functions over subsets of $\mathbb{R}$, we can use $Z$ to integrate functions over subsets of $[0,2\pi]$ (such as all of $[0,2\pi]$).

Let $\mu$ be Lebesgue measure on $\mathbb{R}$, and let $\nu$ be another measure on $\mathbb{R}$. $\nu$ is said to be absolutely continuous with respect to Lebesgue measure if there is a function $f$ such that $d\nu = fd\mu$, or the measure $\nu(A)$ of $A$ is equal to
\begin{equation}
\nu(A) = \int_{A}f(x)d\mu(x).
\end{equation}
The function $f$ is called the Radon–Nikodym derivative of $\nu$ with respect to $\mu$.

Not all measures are absolutely continuous with respect to Lebesgue measure. The example most familiar to electrical engineers is Dirac measure. Lebesgue measure assigns measure zero to any set consisting of a single point, and a measure that is absolutely continuous with respect to Lebesgue measure must do the same. But the Dirac measure $\delta_0$ assigns measure 1 to the set $\{0\}$ and to any set that contains $0$. Since $\delta_0$ is not absolutely continuous with respect to Lebesgue measure, $d\delta_0$ cannot be written as $fd\mu$.

There are also more exotic measures that are not absolutely continuous with respect to the Lebesgue measure.

I have found no evidence of the notion of absolute continuity of stochastic measures.

EDIT: While theoretical results about spectral representations of WSS processes are crucial for applications, the $dZ$ notation may be off-putting and perhaps even doubt-inducing. I suspect that writing $Y(\omega)d\omega$ for $dZ(\omega)$ is a useful abuse of notation that allows the user to manipulate symbols as though some analogue of the Radon-Nikodym derivative existed. Rigor can be added after the fact.

Note that rigor might arrive decades after the fact. Plenty of ideas seem to work just fine without complete mathematical rigor.


If I line in the sound from a video game console, can I detect when a sound is played?

Theoretically you could separate your desired sound by ICA method or other blind source separation methods from the background music, then design a matched filter by your separated sound. The matched filter have high response when the desired pattern (which is your desired sound) exist in the input even when it buried in noise (background music).

The matched filter has acceptable performance but i'm not sure about the ICA and other source separation methods and highly depend on your input signals (desired and background sound).


DTFT of sawtooth wave through DTFT of rect signal

The discrete Fourier transform (DFT) coefficients of a sequence $x[n]$ of length $N$ is defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

and the inverse transform is

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

The discrete-time Fourier transform (DTFT) of a sequence $x[n]$ of length $N$ is given by

$$X(e^{j\omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{3}$$

Plugging $(2)$ into $(3)$ gives

$$\begin{align}X(e^{j\omega})&=\sum_{n=0}^{N-1}\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{-jn(\omega-2\pi k/N)}\\&=\sum_{k=0}^{N-1}X[k]\frac{1}{N}\sum_{n=0}^{N-1}e^{-jn(\omega-2\pi k/N)}\tag{4}\end{align}$$

With the DTFT of a rectangular sequence $r[n]=1/N$ for $n\in [0,N-1]$ and zero otherwise

$$R(e^{j\omega})=\frac{1}{N}\sum_{n=0}^{N-1}e^{-jn\omega}=e^{-j(N-1)\omega/2}\frac{\sin(N\omega/2)}{N\sin(\omega/2)}\tag{5}$$

Eq. $(4)$ can be written as

$$\sum_{k=0}^{N-1}X[k]R\left(e^{j(\omega-2\pi k/N)}\right)\tag{6}$$

So $(6)$ is the correct version of the formula, and the DFT coefficients $X[k]$ are defined by $(1)$. If you use those formulas correctly it should work.


Shannon's entropy for non-equiprobable symbols

You need to distinguish between the entropy of your underlying distribution and the "entropy" of the realization.

Consider a set $\mathcal{X}\subset \mathbb{N}$. This set defines, which symbols are allowed to be contained in a random sequence $X=(x_1, \dots, x_N)$. In this sequence, each element $x_i\in\mathcal{X}$ is drawn from the distribution over $\mathcal{X}$ (i.e. there is a function $p(x)$ which gives the probability of occurence for each symbol $x\in\mathcal{X}$). As a constraint, you have $\sum_{x\in\mathcal{X}}p(x)=1$ (i.e. it is certain that one symbol out of $\mathcal{X}$ is chosen for each $x_i$.

Now, the entropy of the random variable for each symbol is defined by $$H=\sum_{x\in\mathcal{X}} -p(x)\log(p(x)).$$
Here, we use the convention $(p(x)=0) \Rightarrow (p(x)\log(x)=0)$, i.e. if the probability for one symbol to occur is 0, it does not influence the overall entropy. This is inline with the limit

$$\lim_{p\rightarrow +0}p\log p=0.$$

So, now to encode a given sequence $b$ (which is drawn from the distribution above), you take each element $b_i$ from the sequence $b$ and associate a bit-word of a length which is (roughly) proportional to its "realization entropy" (sorry, I forgot the common term for this) $p(b_i)\log(p(b_i))$.

So, on average, the encoding for each symbol has $H$ bits. Hence, to encode $N$ symbols, you'd need (on average) $NH$ bits.


Accurate description of the spectrum of sinusoid waveform added with WGN?

Notice that in your first image, there are peaks where they would be expected to be, so that's fine. I guess you think that "the spectrum obtained through MATLAB seems not like the theroetical one" because of the noisy spectrum you get, instead of a straight line at $\sigma^2$ on the vertical axis.

The random signal added to $s(t)$ is not random actually, it's just a realization you get when calling awgn(). So the FFT will return just a noisy spectrum for that signal. If you generate lots of realizations of $r(t)$ and plot the average of all their FFTs, you'll obtain what you were expecting for.

For your second question on how to calculate the threshold SNR, some criteria is needed. If you average lots of FFTs the variance of the noisy part of the spectrum will drastically decrease, so the SNR needed to see the peaks will be much smaller. If you use just one realization, then the SNR will be higher. There is not a fixed number for the needed SNR, it depends on what you refer to "sinusoid covered by noise". If the spectrum observed was just like the theoretical one (i.e. you averaged a lot of FFTs), a sinusoid would be "completely covered by noise" when its amplitude is such that the corresponding peak in the frequency domain is below $\sigma^2$.


Can't make sense from VHDL butterworth filter implementation

The filter structure is a digital leapfrog and the structure looks like this picture (note: the picture is a different order than the code):

These filters are discussed at some length on wikipedia and wikibooks.


Fixed point processing is not for the faint of heart. It typically requires managing scaling factors, clipping point, signal to noise ratios, underflow, limit cycles, rounding behavior at each single point of your processing chain.
Your best shot at dealing with different data ranges is to rigorously use the Q-number format for each data, state and coefficient in your processing (see https://en.wikipedia.org/wiki/Q_(number_format)) and do the correct multiply & shift  & round/truncate at each operation. This may require statistical or worst case analysis to determine the clipping points and best Qs for each variable.
In order of complexity, you can consider

Convert your fixed point data to floating point, run the algorithm in floating point, and then only limit and convert back to fixed point as the last step.
Same as item one but if you can't do float consider using 24-bit or 32-bit fixed point data instead. This creates a extra margin and headroom and makes life A LOT easier
Stick with 16-bit all around. That means you have to fight for every single dB of signal to noise ratio at every step of your processing. Depending on what you are doing you may consider selectively using double precision for critical state variables or operations.


Convert old computer cassette signal into binary

First of all, you'd go on an information scavenger hunt through the internet on that cassette's data format: http://beebwiki.mdfs.net/Acorn_cassette_format tells us you're quite spot on:

Signal
The signal may be in one of three states; zero, one or no carrier.
Breaks in the carrier are detected and used to reset the cassette
loading firmware.
Zero is represented by a sinusoidal wave at 1200 Hz. The exact
frequency is nominally 16 000 000 / 13 312, or 1201.9 Hz, but tape
decks vary in speed so small differences are tolerated. One cycle (at
1200 baud) or four cycles (at 300 baud) represent one zero bit.
Binary one is represented by a sinusoidal wave at 2400 Hz. The nominal
frequency is closer to 2403.8 Hz. Two cycles (at 1200 baud) or eight
cycles (at 300 baud) represent a one bit. An odd number of 2400 Hz
cycles can and does occur.
To allow the tape deck circuitry to settle, each data stream is
preceded by 5.1 seconds of 2400 Hz leader tone. This is reduced to 1.1
seconds if the recording computer has paused in the middle of a file,
or 0.9 seconds between data blocks recorded in one go.
At the end of the stream is a 5.3 second, 2400 Hz trailer tone. This
is reduced to 0.2 seconds when pausing in the middle of a file (giving
at least 1.3 seconds' delay between data blocks.) The timings are
derived from VSYNC interrupts so they vary between recordings.

So, yes, this is Frequency Shift Keying (FSK) and can be demodulated as such.
Considering modern PCs'/laptops'/phones'/fridges'/toaster ovens' computational power is basically huge compared to the data that fits on a cassette, we don't need to be very considerate on how efficiently we demodulate this, at all.
So, you'd go and

design two bandpass filters: one that's centered around 1201.9 Hz, the other one at 2403.8 Hz,
take their outputs, convert them to their absolute value, then
weight one of them so that they have the same amplitude when excited,
subtract them and finally
make a 0/1 decision based on the sign of the result to get the bit.

Sounds a bit taunting, but looks roughly like this in GNU Radio 3.9:

(You can ignore everything that says "QT" in its name, it's just for interaction / vizualization. Once you've figured out the right amplitude so that 0 and 1 bits come out the same (and that's really just necessary to minimize the probability of bit errors), you'd just multiply with that fixed value instead of having the Qt range.)
You'll have to add some discarding for the start and stop bits surrounding each byte, and you'll need to keep synchronized to these – I think the most elegant (and easy) way would be to write a quick few lines of python code in an "embedded Python Block" in GNU Radio Companion to detect the start of the start bit, then pass on through the 8 data bits, ignore the stop bit and go waiting for the next start bit. Pretty doable!


ORB Algorithm: What Is Intensity Centroid?

p and q - order of the moments, it is the analogue of the mechanical moments (i.e. of force or inertia). Image moments are very basic properties of image, invariant to rotation, could be used as simple descriptor.

Example,

Moment of inertia Ix= I2 = mass *x_coordinate^2.

Image moment Ix2 = intensity * *x_coordinate^2.

http://en.wikipedia.org/wiki/Moment_(physics)

http://en.wikipedia.org/wiki/Image_moment


No, energy is not the same as power.  Power is energy per unit of time.  An analogy- energy is like distance, while power is like speed.

The easiest way to determine energy in the frequency domain is to take the signal and do an FFT.  The Fourier Transform preserves energy, i.e. the transform has the same exact energy as the signal does in the time domain.  You can then find the energy of each frequency bin by multiplying its value with its complex conjugate.

Signal power represents how quickly it delivers energy.  In the context of audio it correlates with how loud the sound is, though there is not a 1-to-1 correspondence.


Is a notch filter for 50 Hz without delay possible?

You could use a 2nd order IIR notch filter as I describe in this post Transfer function of second order notch filter - That post demonstrates a 50 Hz IIR notch with 1 KHz sampling.

[Update: As @user47050 astutely points out in the comments, the IIR notch would also have minimal delay regardless of notch bandwidth, since the dominat delay in the IIR notch filter is specific to the signal components in the rejection bandwidth of the filter, not the signals that you want to pass through. The advantage of what follows is not necessarily then to eliminate a delay issue, but to provide tracking rejection of just the interference that is approximately 50Hz (but slowly varying in phase and amplitude) while minimizing impact to adjacent frequencies of interest.]

However to provide cancellation of the interference with zero delay, consider a feedback signal cancellation servo such as I depict in the diagram below. This would have no delay to your signal but will have a delayed response to changes in your interfering signal. Results will depend on the relative stability of the interference. If the interference is a relatively slowly changing sinusoidal tone, this should work very well.

This is an all digital nulling loop that can null down to the rms power level of your signal (at which point it can no longer estimate the cancellation coeeficient). The bold lines shown in green represent complex signals:The 50 Hz PLL is a NCO with sine/cosine outputs, the multiplier following it is a full complex multiplier (essentially a vector modulator providing complete control of amplitude or phase), the multiplier after the second BPF is essentially a quadrature down-converter providing a measure of the amplitude and phase of the residual.

The way this works can be described as follows: The input signal is assumed to have a strong 50 Hz signal, strong enough such that combined with the filtering of the first 50 Hz BPF the PLL can lock in phase to the tone (the tone needs to be 6 dB stronger than the noise at the PLL input).  This will provide a phase locked copy of the 50 Hz signal as a quadrature NCO output. The quadrature PLL is basically a very tight tracking filter, providing a clean replica for purposes of cancellation.

Next we have a cancellation loop provided with a post cancellation band pass filter to detect the residual and adjust a complex coefficient to adjust the gain and phase of the 50 Hz PLL output prior to subtraction. I show the approach of quadrature downconverting the 50 Hz BPF output to complex baseband (by multiplying with the complex PLL output), and accumulating the result which will find the null with the convergence gain coefficient k. This is a cancellation loop using the method of stochastic gradient descent. The accumulator accumulates residual error as measured on the I and Q output of the quadature downcoverter resulting in the correct setting to the "Vector Modulator" multiplier resulting in a real sine wave output that will be equal in magnitude and opposite in phase to the interfering signal at the input.

The diagram can be simplified in that the first 50 Hz PLL BPF can be eliminated and the PLL replaced with a constant 50 Hz sine/cosine generator, since the vector modulator provides complete control of rotating phase which serves to to adjust the frequency positive or negative as needed, as well as gain as needed to provide the cancellation required. The proposed use of the 50 Hz BPF + PLL is that it would be feasible of providing faster corrections to phase changes in the input signal, once the cancellation loop has converged (since the PLL will track the input and can be designed with wider loop bandwidth than the nulling loop)- however if the amplitude change in the interference dominates, then I don't see a strong advantage to having the PLL.

Without the PLL this would reduce to the following block diagram. Here everything is done at complex baseband, so a low pass filter instead of a bandpass filter is used (which I would likely implement as an exponential averaging filter).  k is the complex convergence gain and c is the determined complex correction coefficient that will scale the gain and phase of the generated 50 Hz cancellation signal (which as suggested earlier c would be a changing phase versus time to adjust for the small frequency offsets that would also exist).

The PLL implementation has a clear advantage over a 2nd order IIR notch filter in that it can track and null a tone interference over a relatively much wider bandwidth, within the tracking range of the PLL. Without further study of the implementation below, I am not confident that it would out-perform a 2nd order IIR notch filter, but it is interesting to consider given it is much simpler than the PLL approach.

As a fixed 50 Hz source with a 2KHz sampling rate, the tone generator is trivial to implement, but if you extended it for other applications and used an NCO for the fixed 50 Hz complex tone generator, you end up with a tuneable single tone adaptive canceller that will capture and cancel any single tone within the proximity to any frequency setting of the NCO, within range of the bandwidth of the LPF! The cancelling loop will null the 50 Hz interferor down to the noise floor coming out of the filter- so the tighter the filter, the deeper the null as long as the the signal is within the filter bandwidth.

The amount of cancellation will depend on the time rate of change of the amplitude and phase of the input signal. The approach as shown above has zero delay to the main signal of interest, but is estimating the characteristics of the interference from a lag time T in the past, where T is impacted by the bandwidth (delay) of the low pass filter and the loop BW of the nulling loop (set by k), so will have an upper frequency limit on how much of the noise bandwidth of the interference signal can be rejected (that noise bandwidth is visible in the power spectral density of the interfering signal, basically AM and PM noise side-bands). This nulling loop is a high pass filter to that interference; which is desired as the signal of interest is indistinguishable from noise.

If this approach proves useful, you can improve it further by instead using the method of recursive least squares (RLS adaptive algorithm) to determine the complex weight that is applied to the vector modulator based on the error terms out of the quadrature downconverter. This converges faster at the expense of computation complexity.


Basically, you can solve this kind of EDO using two different techniques (which also have c different meaning).

Classical Method - you will need the full EDO, containing input and outputs with initial conditions applied over the whole solution y(t), where y(t) = yhomogeneous(t) + yforced(t). The first term, you solve the homogeneous EDO considering NO CONDITIONS. Then, following the pattern of input you exactly find the forced contribution. Then, you add both contributions and apply the initial conditions...
In this technique you totally lose the feeling about the system itself.

zero state response and zero input response. In this case, the whole solution is written as y(t) = y0(t) (zero input response) + x(t)h(t)  (zero state response, as convolution of both signals. The text cutof the asterics).
To find y0(t) you do not consider the input x(t), and apply condition at t= 0^-. It happens because it is completely independent of the input. The system responds just because of the storage energy.
The second term, zero state response, it is the convolution of input and impulse response of the system. So, you need to find the impulse response of EDO. It is done considering the whole EDO , e.g., D4y -D3y +2D2y -2Dy +y = Dx -3x . First you have to solve the homogeneous equation and apply some special conditions generalized... y=Dy=D2y=0 and D3y=1. Lets call it yn(t).... Then you apply the input operator over this yn(t). It means,
D (yn(t)) -3yn(t)... After that, you finally have h(t)...
So, you do the convolution of h(t) and x(t) and find the zero state response
y


Using the linear operator to check for time invariance of a differential equation?

Let $y(t)=\mathcal{H}(x(t))$ describe the system by relating its input $x(t)$ to the output $y(t)$. Then, to check linearity simply check with $x(t)=1$ (i.e. a constant). You will see, that

$$\mathcal{H}(x(t)=1)+\mathcal{H}(x(t)=1)\neq \mathcal{H}(x(t)=1+1)$$

which is one requirement for the linearity. So, the system is non-linear. Actually, you can immediately see this from the fact that $x(t)$ counts quadratically, which is a non-linear relation.

Regarding time-invariance, you see that the time-argument is the same for all functions and the time-value is not used in the equation itself (such as $t\cdot x(t)$ would ne not time-invariant). So, the system is time-invariant.


Causal and non-causal systems in real world

In discrete-time systems, causality is a requirement only when processing (filtering) signals in real time; i.e. when the index $n$ corresponds to a physical time $n \times T_s$. In this case, a non-causal system is impossible to implement since calculating the current output would require inputs from a future time.

In practice, however, many systems are not real time.

Image processing, for example, has indexes for horizontal and vertical coordinates. So it's OK to apply a non causal filter to an image (usual examples being 0-centered FIR filters for blurring or sharpening, or IIR filters applied in both directions to obtain 0 phase distorsion).

In time sequences, many times the signal is stored and then processed in non real time, so you may apply non-causal filters there also.


Understanding a lowpass - comb filter implementation

Your diagram looks correct. Let's call the transfer function in the feedback loop $G(z)$. Consider the signal $w[n]$ at the input to the delay line. Its $\mathcal{Z}$-transform satisfies

$$W(z)=X(z)+G(z)Y(z)\tag{1}$$

where $X(z)$ and $Y(z)$ are the $\mathcal{Z}$-transforms of the input and output sequences, respectively. The output is just a delayed version of $w[n]$, i.e.,

$$Y(z)=W(z)z^{-M}\tag{2}$$

From $(1)$ and $(2)$, the transfer function of the complete system can be derived as

$$H(z)=\frac{z^{-M}}{1-G(z)z^{-M}}\tag{3}$$

Now it remains to find $G(z)$. It is quite straightforward to show that

$$G(z)=\frac{d_2f}{1-d_1z^{-1}}\tag{4}$$

Combining $(3)$ and $(4)$ gives the desired result with $d_1=d$ and assuming that $d_2=1-d_1$.


What is the difference between multiplying a delta and a step versus convolving a delta and a step?

First of all you need to see whether you are performing these operations for a continous time signal or discrete time signal.

Sampling theorem says that multiplication of a signal $x(t).\delta(t)$=$x(0).\delta(t)$ provided $x(t)$ is continous at $t=0$. But here in your question $x(t)$ is a unit step function which is not continous at $t=0$. Hence the multiplication of $\delta(t).u(t)$ is not defined.

In the discrete case the same property is $\delta[n].x[n]$=$x[0].\delta[n]$ and there is no question of continuity as signal is discrete. In this case your question can be written as $u[n].\delta[n]$=$u[0].\delta[n]$=$\delta[n]$ as $u[0]$ is equal to 1.

Now the convolution property says that $x(t)*\delta(t-t_0)$=$x(t-t_0)$ where $*$ is denoted here as the convolution operator. In your question this would be $u(t)*\delta(t)$=$u(t)$ which is nothing but unit step function for continous time.

The same is true for discrete signals i.e.$x[n]*\delta[n-n_0]$=$x[n-n_0]$. So for discrete case your question would be $u[n]*\delta[n]=u[n]$ which is a unit step function for the discrete case.

You can look at comments of Matt below my answer for better clarity.


I have solved the problem.

Since I was taking the sum of 4 pixels in downsampling, what I was doing in effect was convolving with a sum of delta functions (that is, one delta function for each pixel). In Fourier space this translates to multiplying with a cosine, which explains the increasingly strong (relative) discrepancy between input and observed power spectrum.

To mitigate this, Wikipedia helped me out. If we first implement a lowpass filter in Fourier space (filtering out high frequencies, to avoid aliasing later on), and then reduce the sampling rate (to 1 in 4 pixels), we get the desired result, and indeed: the observed power spectrum looks nice up until the cutoff:


Bandwidth expansion by computing missing harmonics in music?

Bandwidth expansion is actually a problem of missing audio data recovery, which is an active area of research.

Usually you assume that your signal follows some statistical model in the time-frequency domain, and the whole challenge is to infer the parameters of your model, based only on the available data.

In the case of bandwidth expansion, you would infer your model from the lower frequency part of your spectrogram, in order to recover the higher frequencies.

You can have a look at some of Paris Smaragdis' work on missing audio data and bandwidth expansion, using latent variable models. Some demos are available here.


Computing the image derivatives are very important for a lot of vision tasks. So let me first begin with answering your questions quickly and then get to more details.

Whether central difference can be applied for finding partial derivative?
YES
Whether the filter mask (which I have proposed above) can be used or not for finding partial derivative?
It can be used, but only the vertical one (in $y$ direction). The horizontal one would be the transpose of that.
If not, please suggest me a filter mask which can perform it better and also mention how is it performing the required operation.

For noisy images computing derivatives this way would not be very helpful. That's why Sobel operator exists (http://en.wikipedia.org/wiki/Sobel_operator). Instead of the mask you proposed, it suggests to use $0.125 \times [1, 0, -1; 2, 0, -2; 1, 0, -1]$. However Sobel lacked the rotational symmetry. So Scharr kernels were proposed (the same wiki page, on the very bottom). These were proven to be better alternatives. Yet, these developments did not stop the further research. Being not so old, Farid and Simoncelli proposed more robust differentiation operators. See the paper:

Hany Farid and Eero Simoncelli "Differentiation of Discrete
Multi-Dimensional Signals" IEEE Trans. Image Processing. 13(4): 496-508 (2004)

One implementation of this is presented in Peter's Computer Vision MATLAB functions:

http://www.csse.uwa.edu.au/~pk/research/matlabfns/

Check derivative5 and derivative7 functions.

Reading about Scale Space would also give you more insight on the nature of such linear operators: http://en.wikipedia.org/wiki/Scale_space

Hope these help. Cheers,


Compare two words from different speakers

There is no (known) magic solution, as the sound of the same word from very different speakers are often not actually similar to one another in term of any simple characteristics such as you list (frequency range, intensity range, pitch, etc.)

To often, humans can only guess at what word is said from an unknown speaker based on context.  And you don't have that with single word recordings.


what is probability density function (PDF) of periodic signal

There is no best pdf for a periodic signal. There is also no way to find the 'exact' pdf of a measured signal. What you have to do is to measure the pdf from the data. Use a histogram to approximate the pdf of your data. Define a number of intervals within the amplitude range of your signal and simply count the numbers of data samples per interval. This gives you the histogram which is an approximation of the pdf of your signal.


Detecting peak value on realtime IMU data stream

In order to detect peaks, the simplest method will be to take the difference of adjacent samples and check for zero crossing. You may also want to smooth out the data before taking difference to prevent spurious peaks due to noise. An outline of the procedure here would be like this

Smooth out the data using Moving average filtering. It is a simple filtering where you take average of, say, 5 neighboring samples to arrive at sample value.
For the above filtered signal take the adjacent sample difference ($y[n] = x[n] - x[n-1]) and save the sample indices where difference hits zero (or near zero).
To differentiate positive peaks from negative peaks, also check for sign of sample values around these indices. Eliminate negative peaks, and remaining ones indicate positions where positive peaks occur.
You may not get exact zero crossing, so you may want to keep a threshold. Also, the position may have offset due to filter delay. A 5-sample MA filter will delay your sample train by ceil(5)=3 samples.
There are ways to combine MA filtering and taking difference. Look for filters which do differentiation (https://www.dsprelated.com/showarticle/35.php)


If the discrete-time filter is described by the linear difference equation

$$y[n]=\sum_{k=1}^{K}a_ky[n-k]+\sum_{m=0}^{M}b_mx[n-m]$$

(with $a_K\neq 0$ and $b_M\neq 0$), the filter order is given by
$\max\{K,M\}$, i.e., by the maximum delay necessary to implement the filter.


How to apply band-pass filter and get the filtered back from frequency to time domain in python?

You shouldn't perform the filtering operation on the FFT data, but on the original time domain data. The functions filtfilt and lfilter both take time domain data as their inputs, not frequency domain data. Filtering can also be implemented in the frequency domain (by multiplication), but the functions you're using don't do that.


Low-Pass Filtering on Small Chunks of Audio

Filtering a data block usually results in more data than fits in the size of the block.  If you throw away this added data, that will produce artifacts across blocks.

For FIR filters, you need to pad each chunk or block with at least the length of the impulse response of your filter before filtering (>= N+M-1).  Then use overlap-add or overlap-save (FFT fast convolution algorithms) to process all the extra results from convolution of the added padding (usually carry them over to the next chunk or chunks).  Otherwise you are either getting circular convolution artifacts in each block, or are throwing away part of the filter response from applying your filter to each block.

For IIR filters, you need to carry over the filters internal state between chunks, and not start over from zero or other default initial state of the IIR.  That state still carries energy from the filter result, and if you throw it away, you will get glitches between chunks.


How do I examine if the signal can be z-transformed?

You can use the basic properties of Z-transform to find the result:

Given, $$ x[n] = \sin(\frac{2 \pi}{N} m n) u[n] = \sin(\omega_m n) u[n]  $$

Decompose the sine wave into complex exponentials:

$$ x[n] = \frac{1}{2 j} \left( e^{j \omega_m n} - e^{-j \omega_m n} \right)   u[n] \tag{1}$$

And observe the following :

1 :
$$ U(z) = \mathcal{Z} \{ u[n] \} = \frac{1}{1-z^{-1}} ~~~~,~~~ |z| > 1 $$
2 :

if  $x[n] \longleftrightarrow X(z)$ then
$$ e^{j \omega_m n} ~ x[n]  \longleftrightarrow X( z ~e^{ -j \omega_m} ) ~~~,~~~\text{ ROC the same} $$

Then apply these two observations to Eq (1) :

$$ \begin{align}
X(z) &= \mathcal{Z} \{ x[n] \} = \mathcal{Z} \{ \frac{1}{2 j} \left( e^{j \omega_m} - e^{-j \omega_m} \right)   u[n] \} \\
&= \frac{1}{2 j} \left( \mathcal{Z} \{  e^{j \omega_m} u[n] \} - \mathcal{Z} \{ e^{-j \omega_m} u[n] \} \right)  \\
&= \frac{1}{2 j} \left(  U( z e^{-j \omega_m} ) - U( z e^{j \omega_m} ) \right)  \\
&= \frac{1}{2 j} \left(  \frac{1}{1-e^{j \omega_m} z^{-1}} - \frac{1}{1-e^{-j \omega_m} z^{-1}}  \right)  \\
\end{align}
$$


The Stability of the RLS Algorithm for Equalization

Short Answer:
Basic RLS (no forgetting, no weird weighting, etc.) is ALWAYS Lyapunov stable. If the regressor sequence for the LS problem is persistently exciting--which is data and problem dependent, not algorithm dependent--then RLS is exponentially stable. So I don't know what you mean by "LMS is more stable than RLS"--more stable in what sense?
More Details:
Basic RLS has the discrete-time update equations
\begin{align}
%
\theta_{k+1} &= \theta_k + P_{k+1}\phi_k^T(y_k-\phi_k\theta_k), \\
%
P_{k+1} &= P_k - P_k\phi_k^T(I+\phi_kP_k\phi_k^T)^{-1}\phi_kP_k.
\end{align}
Note that the covariance update is decoupled from the estimate update and is also monotonically decreasing. If the regressor sequence is persistently exciting, then the origin of the covariance update is exponentially stable. If nothing can be said about the regressor, then the origin of the covariance update is only guaranteed to be Lyapunov stable.
Next, define the error $\tilde \theta_k = \theta_k-\theta$, where $\theta$ is the true value of the parameters. Then you can show that
$$\tilde \theta_{k+1} = P_{k+1}P_0^{-1}\theta_0.$$
In other words, the error is entirely dictated by the covariance, and hence is also exponentially stable or Lyapunov stable as the covariance is.


What is the relationship between group delay and propagation delay for a filter?

The quantity you might be looking for is signal-front delay, which is the delay of the beginning of a signal passing through a linear system. It is simply the largest value $\tau_{sf}$ (in samples) for which
$$h[n]=0,\qquad n<\tau_{sf}\tag{1}$$
is satisfied, where $h[n]$ is the system's impulse response.
If an input, or change in the input, starts at $n_0$, then the system's response to it begins at $n_0+\tau_{sf}$.
If the impulse response has a long period of pre-ringing, the signal-front delay is generally much smaller than the actual propagation delay. In that case, definition $(1)$ can be modified as
$$\big|h[n]\big|<\delta,\qquad n<\tau_{sf}\tag{2}$$
where $\delta$ is a positive constant which is small compared to the maximum value of $h[n]$.
Depending on the nature of the impulse response, another meaningful definition of signal delay introduced by a linear system could be the center of gravity of the impulse response:
$$\tau_{gr}=\frac{\displaystyle\sum_n nh[n]}{\displaystyle\sum_n h[n]}=\frac{jH'(0)}{H(0)}\tag{3}$$
where $H(\omega)$ is the system's frequency response, and $H'(\omega)$ is its derivative w.r.t. to $\omega$. Note that $\tau_{gr}$ is simply the group delay evaluated at $\omega=0$.
Using the center of gravity as an estimate for signal delay is only useful for systems with a lowpass characteristic. For other systems (such as highpass or bandpass filters), the denominator of $(3)$ could be zero or very close to zero. For that reason, for general systems one would exchange $h[n]$ in $(3)$ for either its magnitude or its square.
Note that unlike phase delay and group delay, neither signal-front delay (including its modified form $(2)$) nor the center of gravity are functions of frequency.


Variation of channel attenuation with frequency

No, there would be no local attenuation over the 10MHz subcarriers at 1800 MHz of any significance if you exclude the multipath effects. The standard path loss equation holds and this would be a non-dispersive (linear phase) channel with minimal difference in attenuation across the 10 MHz band. However, due to multipath fading you can see significant changes depending on the delay spread of the signal due to the fading environment (frequency selective fading).

See the figures below to get a sense of the delay spreads at 1.8 GHz for different environments and its effects on fading.

In the simple model of a direct path and a delayed path both of equal amplitude, you would have frequency nulls spaced at $f= 1/T$ where T is the delay in seconds and f is the frequency spacing in Hz. This gives you a sense if the delay spread involved will result in flat or frequency selective fading. For the three typical delay spreads listed relevant to your 10 MHz channel, the spacing is 5 MHz (for 0.2 us open areas), 2 MHz (for 0.5 us suburban) and 333.3 KHz (for 3 us urban). Of course the nulls won't be so evenly spaced as in the figure as multiple delay paths are involved, but this helps give a first order metric as to the frequency selectivity of the fading.


Voice masking for removing voice and not birdsong

Sound masking in time-frequency is a very complex problem, especially to do it in real-time while collecting data. For ensuring privacy I would suggest using a simpler approach.

Sparse random sampling

Instead of recording continuously for hours (potentially picking up whole conversations), record only short segments spaced from each-other in time. They need to be long enough to capture events in question (birdcall/birdsong), but short enough to not be able to capture much coherent conversation.

For instance record 0.1-1 second segments every 10-100 seconds. This can be done on a fixed schedule, or using reservoir sampling.

This approach has the advantage of not introducing any bias in the selection. It reduces the amount of data, and limits analysis to short time windows. It may also meet the privacy protection requirements.

Drop samples with Voice Activity Detection

While recording, consider segments of 1-60 seconds. Run a standard Voice Activity Detection (VAD) algorithm on each of them. If voice is detected, discard the entire segment. Need to use a low threshold for the VAD to get a low amount of false negatives (sample actually has speech).

Sampling only when the VAD does not trigger may introduce bias in the data. For instance if there is humans talk at certain times of day. Or if some birdsong is confused as human voice by the VAD.


Question about ASK, correlation, and frequency domain

You are mixing different notions.

What happens in basic digital transmission is :

Digital signal generation (0's and 1's)
Bit to symbol association (digital modulation), in our example ASK.
Shape filtering, which allows the signal's spectrum to be contained within a standardized shape, since we can't just throw a digital signal with very sharp transitions in the air, in would pollute the electromagnetic spectrum.
Finally, we mix the signal with a carrier frequency in order to shifts it's spectrum around the appropriate transmitting frequency

What happens at the receptor is the same operations in anti-chronological order :

Analog demodulation (the signal is mixed with a carrier of the same frequency)
Matched filtering, which is the same operation as in the emitter, but with a specific filter, which is ideally the complex conjugate and temporally reverted of the filter used in emission ($g^*(-t)$). Why is it this one and not another ? See Nyquist ISI criterion.
Symbol estimation
Digital demodulation of estimed symbol

About the filtering operation specifically : time domain and frequency domain operations are equivalent. If you do a time domain filtering, or any mathematical operation (under trivial notations) :
$$
y(t) = g(t)*x(t)
$$
It is equivalent to frequency domain multiplication :

$$
Y(f) = G(f)X(f)
$$

Analog amplitude modulation of a signal $r(t)$ is, like you say, the operation $r(t)\cos(2\pi f_0t)$ (multiplication) which in frequency domain is equivalent to $R(f)(1/2)*(\delta(f-f_0)+\delta(f+f_0))$ (convolution).

Demodulation takes $r(t)\cos(2\pi f_0t)$ and multiply it again :
$$
s(t) = r(t)\cos(2\pi f_0t)\cos(2\pi f_0t)
$$
Which gives after simplifications :
$$
s(t) = r(t)\frac{1}{2}\left(1+\cos(4\pi f_0t)\right)
$$
That's when your LPF comes into action : it filters the unwanted $2f_0$ component and keeps only $r(t)$ which is the demodulated signal, and to do this mathematically, we write $s(t)*lpf(t)$, $lpf(t)$ being the LPF's impulse response.

In short : analog demodulation is either a multiplication in time domain ($r(t)\cos(2\pi f_0t)\cos(2\pi f_0t)$) followed by a convolution in time domain ($s(t)*lpf(t)$)

OR

It is a convolution in frequency domain ($R(f)*M(f)*M(f)$) (with $M(f)=(1/2)*(\delta(f-f_0)+\delta(f+f_0))$ followed by a multiplication in frequency domain ($S(f)LPF(f)$)

BUT both methods outputs are the same, physically speaking. We work in frequency or time domain according to which domain is best suited to give easy calculations, as you may have noticed.


First off- I know you are doing you're best (which is why I am answering) but your question is very hard to read.  Regardless of your English skills, it would be easier to read if you would break it up into paragraphs and throw in some pictures to help show what you're talking about.

  These 1023 chips (vector’s elements) should be repeated every 1 ms in order for the chipping frequency to be achieved (10Mchips/s).

1023 chips repeated every 1 ms means you have $\frac{1023}{.001 s} = $1 Mchips/s, not 10 Mchips/s.

  My first problem is that in the first place in order to demonstrate the spectrum of the codes,spectrul nulls should be represented at frequencies multiples of 1.023 Mhz.The only way that I could somehow think how to do this is to use massive vectors where are not really supported by my laptops memory (RAM).

You should not need a "massive" vector to demonstrate this.  Increasing the number of samples that you FFT just gets you better frequency resolution- it doesn't change the frequency range.  I think that you should be able to see the nulls just fine with several repetitions of the spreading code.

  I say somehow because the sampling rate should change (for example wherever I have one ‘1’ I should add some more and the same applies to zeros).

I have no idea why you say that the sampling rate should change, unless you're talking about going from the unspread rate to the spread rate.


About how to increase the FIR filter sampling frequency in FPGA? And what is the tradeoff of increasing the sampling frequency?

Does the “one tap” delay equal to “1/Fs (sampling frequency)”?

Yes

  If I have N-tap, the total delay will be N/Fs?

Depends on how you define "total" delay, but in general the answer is no. For a minimum phase filter the delay will be 1. For a linear phase filter it would roughly be N/2

  If the Fs sampling frequency is increased, the “one tap” delay is decreased?

Yes. But it also changes the frequency response of your filter.

  Is there any trade off/drawback of increasing the sampling frequency?

Many. It really depends on what the requirements of your application are.


How to remove color cast from HSV information

It is not necessary that your device's color mismatch can be corrected with a simple shifting of the Hue component (except of course if you are absolutely sure that this is the case with your equipment).

In general, to correct the response of the device you somehow need to recover the mapping that it already performs between reference and observed colors and then see how it can be corrected by some curve (monotonic or not).

For this purpose you will need a colour chart (and ideally, an instrument such as this one).

(Please note that it is best that this correction is applied at the system level, i.e. your application will still be sending some "color demand" normally, without correcting it but that demand would undergo some correction prior to display automatically by the system.)

At the very basic level, you can send a simple reference image to your display and then try to adjust its R,G,B settings (if there are any) so that the colours seem to display properly (according to your eyes that is)

If you have a system that is not capable of applying such profiles but you still must apply some color-correction you could try matching the relative distances of the colours with a known target (or image in your case) and then (if required) try to correct towards absolute colours by shifting (with trial and error).

Take an evenly lit image of the colour chart with a good digital camera, color correct it (to reduce any effects of the lighting or sensor used to acquire your "reference" image) and feed it to your display. Take an image of the display (ideally in the same lighting conditions), apply the same color correction with the profile from the first step (and for the same reasons) and then examine the relative color correspondence between any two colors from the color chart image and two colors from the color chart display image to recover the mapping. (A piece of software that might be of help in this process is LPROF)


Why is there little emphasis on phase distortion?

Why do we not seem to care about the nonlinear phase response in certain applications?

Very likely because these certain applications don't care about phase!

  Is it possible to correct the phase distortion if we are filtering offline (i.e. not in realtime)?

With another filter that, combined with the distorting filter, has the desired phase properties, yes. The question whether such a filter exists is actually a complicated one and depends on the distorting filter!

  Given that we have a model of our filter, can we take the outputted phase distorted signal, and then individually shift the frequencies within the signal?

This is where the complications come from: If you pick "single" frequencies only, you could only represent perfectly periodic, infinite signals (limit of bandwidth-time product in the Fourier transform). So you'll need some kind of system that is "continuous enough".

  It also seems like for application where we can't allow phase distortion, we often implement FIR filters. What are some applications where FIR filters would be absolutely necessary? Essentially anything being filtered in realtime?

No; real time isn't the point here.

So, there's many applications where you want a linear system.

Which practically reduces your choices to IIR and FIR filters.

However, the question whether to choose either of these is often one of design possibilities (I personally find it relatively easy to design a FIR with the window methods, and relatively hard to design good IIRs, but that's me, who's grown up in a FIR-friendly environment); it's often also a question of how stable an IIR design can stay under limited bit width assumptions, or how long a FIR becomes.

Overall, the choice of digital filters is a along topic, that actually fills shelves in libraries. It's impossible to just say "we use that", without knowing a lot of the requirements set for a filter.

Having said that, in many cases you just really need something that is easily adaptable. FIRs lend themselves to that excellently, because you'd usually "just" have to exchange some tap values in software, and depending on the actual implementation of the FIR, you might even have something like a continuous response during that exchange.

Note that all questions regarding which type of filter you use typically involve tradeoffs – your filter response isn't perfect, but you can't implement infinite length FIRs anyway. Your latency is low, but you paid for that with a suboptimal suppression in the stop band. Your decimating filter has adjustable decimation, but in exchange its CIC response is rather round.

A lot of these tradeoffs involve the actual thing you need to implement your filter on. You really shouldn't assume that an experienced DSP hardware designer would just take your Matlab FIR and put it in his FPGA/ASIC design. First of all, there's typically no floating point in those designs, so he'd first need to check dynamic range, convert to fixed point, estimate the error incurred by that, adjust for that, and secondly, things that might be hard in matlab might be easier in hardware; changing the path that samples travel based on the sign bit at Megasamples per second? Will totally thrash CPU performance if you constantly exchange taps in CPU caches, but for FPGA designs, that doesn't matter, since the hardware needs to be there, anyway. Using an IIR in software? What's the point, linear access to memory is preferable and caches are fast, so if the filter's significantly long, it'll probably be fastest to use a FFT-based FIR implementation. In hardware, that'd be a horrible waste of space if the same can be done with an IIR of half the length.

Coming back to the question of phase distortion and filter choices: We know it's very simple to test whether a FIR filter has linear phase: iff its taps are symmetric (hermitian in the complex case), then it has linear phase. That is a strong tool, also for designing linear phase systems (don't confuse with linear systems). Approximating linear phase with an IIR can get pretty complicated / impossible, and getting linear phase with a nonlinear filter – I'm a communications engineer; non-linear filters are pretty much frequency mixers, so we, most of the time, avoid those, if we don't intend to drop phase information anyway.


Extreme value at start of FFT proccesed signal

This is just a quantity proportional to the 0-th frequency of your signal, in other words the average value. It is not abnormal, since it is necessary to have an invertible Fourier transformation.

If you want to remove it, you can just remove a central tendency (like the mean) from your signal (windowed or not).


Help with deconvolving photon reads distributed across neighboring pixels

There is some work on resolving point sources from their smoothed signal called super-resolution, along with some elegant methods to solve it using convex-programming. The idea is that if you know that your signal consists of spikes (single points of energy), you can resolve them from the smoothed signal, even when the smoothing causes interference between spikes.

I used this method in my PhD to find the amplitude of 1D spikes from low-pass filtered signals, using the research published here (not mine). I didn't look at 2D signals like yours, but I recall that some work has been done on it, so you might be able to find some code.

For your method, I'm guessing the negative pixels occur because the smearing isn't exactly modelled by your PSF

What about something simple:

Assume that the smear occurs within a 3 x 3 radius.
Convolve the image with a 3 x 3 filter of all ones. This will add all the values in that area.
Find the points of local maxima in a 3 x 3 radius in the original image.
Replace the value of those points with those of the convolved image and zero the rest.


Difference between "Discrete Wavelet Transform" and "Discrete Wavelet Decomposition"

The discrete wavelet transform should denote "the operations" that, applied to some data, yield a discrete wavelet decomposition. The first one can be seen as a matrix operator, while the second relates to the actual wavelet coefficients, or the structure of thereof, that you would obtain after the application of the first one.

In everyday language, they are often used interchangeably, by a form of metonymy:

  a figure of speech in which a thing or concept is referred to by the
  name of something closely associated with that thing or concept


If somebody squares a real signal in the time domain, he probably deals with the signal's energy. Since this is a pointwise product, this results in a convolution in the Fourier domain:

$$ \mathcal{F}(x\times x) = \mathcal{F}(x)\ast  \mathcal{F}(x) $$

In some cases, as stated by @MBaz, sometimes trigonometric identities can explain the frequency shift, as:

$\cos^2 x=\frac{1}{2}(1+\cos2x)$
$\sin^2 x=\frac{1}{2}(1-\cos2x)$

From the above formulae, one can see that:

they are positive, and the non-oscillatory term is $1/2 \cos(0x)$, explaining a peak at frequency zero,
the scalar before $x$ has doubled, explaining the peak moves.