About Fourier transform of periodic signal

We can figure out what's going on if we first understand a simple identity and then just compute the Fourier transform of the periodic function.

A useful identity

First let's prove that
$$D(\omega - \omega') \equiv \int_{-\infty}^\infty dt \, e^{i (\omega-\omega') t} = 2\pi \, \delta(\omega - \omega')\,. $$
We just use a test function $\tilde{g}(\omega)$:
\begin{align}
\int_{-\infty}^\infty \frac{d\omega}{2\pi}
\underbrace{ \left(\int_{-\infty}^\infty dt \, e^{i (\omega-\omega') t} \right)}_{D(\omega - \omega')} \tilde{g}(\omega)
&= \int_{-\infty}^\infty dt \, e^{-i \omega' t}\left( \int_{-\infty}^\infty \frac{d\omega}{2\pi} \tilde{g}(\omega) e^{i \omega t}\right) \\
&= \int_{-\infty}^\infty dt \, g(t) e^{-i \omega' t} \\
&\equiv \tilde{g}(\omega') \, . \\
\end{align}
This shows that by definition $D(\omega-\omega') = 2\pi \, \delta(\omega - \omega')$.
We didn't have to say anything weird about distributions or use any magic.

Fourier transform of periodic function

Now consider a periodic function $x(t)$ with period $T$.
Since $x$ is periodic we can write it as a Fourier series
$$x(t) = \sum_{n=-\infty}^\infty \hat{x}_n e^{i 2\pi n t/T} \, .$$

Now let's compute the Fourier transform,
\begin{align}
\tilde{x}(\omega)
&\equiv \int_{-\infty}^\infty x(t)e^{-i\omega t} \, dt\\
\text{plug in Fourier series} \qquad &= \sum_{n=-\infty}^\infty \hat{x}_n \underbrace{\int_{-\infty}^\infty dt \, e^{i(2\pi n/T -\omega) t}}_{D(2\pi n / T - \omega)} \\
&= \sum_{n=-\infty}^\infty 2\pi \, \hat{x}_n \delta(\omega - 2\pi n / T) \, .
\end{align}
So there we go: the Fourier transform of a periodic function $x(t)$ is a sum of delta functions appearing at $x(t)$'s Fourier series frequencies (i.e. the frequencies $2 \pi n / T$).
The weight of each delta function is $2\pi$ times the Fourier series coefficient $\hat{x}_n$.
The factor of $2\pi$ is there simply because we used the Fourier transform convention with angular frequency $\omega$ but the Fourier series convention with regular frequency.


If the image is given in the conventional format of sampled pixel values of an sensor matrix, then no it's not possible to apply any software means of EM wavelength/frequency processing on those pixel data. The reason is clear: the sensor pixels are not capturing time variation of electromagnetic wave, they only capture the static (or average) intensity at a point and collectively the intensity distribution along the sensor surface points. Wavelength (and frequency) sensing at a point requires a time domain capture, such as performed by the RF radio antenna or a microphone (sensing the acoustic pressure wave), which therefore enables subsequent frequency/wavelength processing possible.

Therefore practically there remains only a hardware lens filter possbility which is indeed a trusted parctice in the field of photography.

Look for image filters from optics part producers, they can supply any type of filter you want (they are a litle pricey however)


I am not aware of a direct comparison, other than the agast paper:
http://www6.in.tum.de/Main/Publications/Mair2010c.pdf

However, Agast is also based on AST, which is developed by the authors of FAST -Edward Rosten (et al.). It is just a better training. So it might be safe to assume the improvements, based on the original Agast paper.

Moreover, here is a tracker I implemented based on Agast features:

https://www.youtube.com/watch?v=rZWu4LxmsME

This works fairly well in many cases and was able to convince me on the power of Agast (Doesn't work as good with FAST features).


The interval resolution should be determined according to data you want to analyze.

It is equivalent to the question, what magnification of your Micro Scope would you use?
Well, it depends on the size of the things I want to look at.

If your data within the sampling frequency you set has different energy in each frequency bin and you'd like to see it, use small intervals.

Though,
The real question, how do you intend to use the Frequency Vector?


if the two audio signals are totally uncorrelated, the the squares of the two crossfade gains should add to 1.  so, to prevent that dip, your crossfade function will look more like $\sqrt{\frac{1}{2}(1+t)}$ and $\sqrt{\frac{1}{2}(1-t)}$ rather than $\frac{1}{2}(1+t)$ and $\frac{1}{2}(1-t)$  for $-1 \le t \le +1$.  note, for constant power crossfades, at the midpoint of the crossfade ($t=0$), both crossfade gains are at $\sqrt{\frac{1}{2}} \approx 0.707$ rather than at $\frac{1}{2}$.  that's what addresses the dip you want to avoid.

if your two audio signals are somewhere between totally uncorrelated and totally correlated, you need to determine the degree of correlation (that's what cross-correlation is all about), and then the crossfade functions can be determinied analytically.  i have posted to music-dsp a little theory about it.  it's not perfectly rigorous, but i think it works.


So I took a look at the slides and it seems to be giving a cursory overview of how to measure Doppler with multiple pulses, so it's very generic. Hopefully this is a pretty straight forward thing to discuss given that your questions are suspicious of their statements, that's good!

First Question: why is the Nyquist theorem applied to the doppler frequency $f_D$? It should be applied to the sampled sine wave, which
is at frequency $f_T$.

The received frequency $f_R$ is the transmitted frequency $f_T$ with the added Doppler shift $f_D$ as $f_R = f_T + f_D$. So in order to measure Doppler unambiguously, you need to be able to sample $f_R$ above Nyquist. But the question still stands as to why the text states that we must appropriately sample $f_D$ instead of $f_R$?
As is done in modern radar systems, and as implied in the text, the received pulses are mixed-down to baseband. That is, the carrier is centered around zero frequency. This means that  $f_T = 0$ and so
$$f_R=f_D$$
Where the Nyquist criterion for the PRF on sampling Doppler is now more apparent.
It is possible to undersample around carrier frequencies that are not zero and still extract the desired signal, but I don't think that's necessarily important for this discussion. It's a slightly more advanced technique that simply achieves the same thing. You can look up "bandpass sampling" for more information.

Second Question: why does aliasing cause folding? As far as I'm concerned, spectral aliasing means that the original (correct) spectrum is replicated and those spectra overlap (and hence it is not possible to separate them correctly).

As you already stated, sampling produces spectral copies that are separated apart by your sampling frequency, in this case which is the $PRF$. Aliasing causes "folding" because as the frequency being sampled goes past Nyquist, in this case $\frac{PRF}{2}$, there is a copy incoming from the left into your un-aliased window that spans $[-\frac{PRF}{2}, \frac{PRF}{2})$.
Let's do a simple example where we have the following properties:

$PRF = 50 \space kHz$

Unambiguous spectrum span: $[-\frac{PRF}{2}, \frac{PRF}{2}) = [-25 \space kHz, 25 \space kHz)$

Doppler frequency $f_D = 10 \space kHz$

Everything is fine here, we're well below the Nyquist rate and we can recover the Doppler frequency nicely:

Due to sampling, the adjacent copies to the left and right this sinusoid exist at $f_D \pm PRF$. Meaning we have two tones that exist at the frequencies
$$f_{left} = 10 \space kHz - 50 \space kHz = -40 \space kHz$$
$$f_{right} = 10 \space kHz + 50 \space kHz = 60 \space kHz$$
As expected, the copies are outside our Nyquist zone achieved by properly sampling.
Now, what if the Doppler frequency was beyond Nyquist? This is where aliasing causes "folding". Let's make the Doppler frequency be $f_D = 35 \space kHz$. The two adjacent frequencies are now
$$f_{left} =  35 \space kHz - 50 \space kHz = -15 \space kHz$$
$$f_{right} = 35 \space kHz + 50 \space kHz = 85 \space kHz$$
We can see here that the left copy has now moved, or folded, into our Nyquist zone! Then due to this aliasing, we measure the wrong Doppler:

A note on your last image:
This is a generic representation of signals that have some finite instantaneous bandwidth. In our case, we're dealing with sinusoids which have some peculiarities when it comes to aliasing since they are ideally delta functions in the frequency domain that don't have a "width" per se.


Energy Spectrum of a signal after FFT in Matlab

The cited passage gives a method for obtaining the energy in the harmonics of a specific frequency.  If your frequency of interest is $f_0$, then take a section of samples of length $N$ where $N = \frac{4 F_s}{f_0}$ and compute the $N$-point FFT.  Now you can find the energy of the second harmonic by squaring the magnitude of the $4k+1 = 9$th FFT bin. To account for errors in estimating $f0$ you can sum the squared magnitudes of bins 7 through 11.

Below is some example code that might give an idea of what the algorithm in question does:

%% Generate an example signal
Fs = 5e3;                     % Sampling rate
f0 = 151;                     % Frequency of interest
t = (0:Fs-1)/Fs;              % Time samples
x = cos(2*pi*f0*t);           % Pure sinusoid @ f0
noise = 0.1*randn(size(x));   % Gaussian noise
y = (x>0)-(x<0) + noise;      % Introduce harmonics and noise

%% Compute energy in the harmonics of f0
N = round(4*Fs/f0);           % Determine record length
Y = fft(y(1:N))/N;            % Compute N-FFT (normalized)
Y2 = Y.*conj(Y);              % Magnitude squared
K = floor(N/2/4);             % Number of harmonic "bins"
for k = 1:K
    idx =  4*k+1+(-2:2);      % Five indices centered at 4k+1
    Pyk(k) = sum(Y2(idx));    % Sum of energies around 4k+1
end

%% Plot the results
[Pyy,Fyy] = pwelch(y,[],[],[],Fs);
plot(Fyy,10*log10(Pyy));
hold on;
F = linspace(0,Fs-Fs/N,N);
Fk = F(1+4*(1:K));
plot(Fk,10*log10(Pyk),'ro');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/rad/sample)')
data2 = sprintf('Energy in Harmonics of %d Hz',f0);
legend('Power Spectral Density Estimate', data2);
hold off;

As far as finding a "signature" in the data, you can see that the algorithm could be very useful for that.  However, it is very important to note that the algorithm presumes some fundamental frequency of the signal in question.  If the signal in question has a fixed period, but the shape of the signal over that period changes with time, then this algorithm would be useful.  However, if the frequency varies with time, then this algorithm is not a good fit for your application.

If you are simply looking to estimate the spectrum of a signal, then maybe this page will be of help: FFT for Spectral Analysis (mathworks.com)


why is the DFS of a delta function equal to 1

$X[k]= \sum_{n=0}^{N-1} \delta[n]W_N^{kn} \quad$, where $\>W_N^{kn}=e^{−j\,2\pi k\,n/N}$

HINT:

What is the value of $\delta[n]$ when $n \neq 0 \>$ ?


Upsampling before performing pulse shaping

The reason is that you want to space the raised-cosine pulses by the symbol interval, $T_p$. Consider the signal you want to create: $$s(t) = \sum_k a_k p(t-kT_p),$$ where $p(t)$ is  your prototype raised-cosine pulse, and $a_k$ are the symbols. Notice that pulses are spaced $T_p$ seconds.

To re-create this with a filter, you can write the equation for the signal as follows: $$s(t)=\sum_k p(t) \ast a_k\delta(t-kT_p),$$ where $p(t)$ is interpreted as the filter's impulse response. Notice that the impulses are spaced $T_p$ seconds, and this results in the raised-cosine pulses being spaced by the same amount.

If you want to simulate this process in discrete-time in Matlab or similar tools, then you can create a vector p for $p(t)$, and a vector d for $a_k\delta(t-kT_p)$. Let's say your symbols are stored in vector a. Then, you can create d by upsampling a in such a way that there are $T_p$ seconds between symbols:

d = [ a[1], 0, ..., 0, a[2], 0, ..., 0, a[3], 0, .... ]

You need to insert as many 0s as necessary, according to your sampling frequency, to have $T_p$ seconds between symbols.

Then, you can generate the vector s corresponding to $s(t)$ with s=conv(p,d).


Octave is a good option. Some people call it as a free version of Matlab.
I found Wolfram Mathematica to be of great use. It seems to be much more powerful than Matlab, not to mention, simpler to use.


Why is the number of frequencies decomposed in scipy.signal.stft() equal to the hop size?

The number of frequencies 'measured' in the Fourier transform of each time frame is exactly equal to the hop size,

The number of frequencies is not determined by the hop size, it's determined by the FFT size.  If you don't specify the FFT size, it uses nperseg as the FFT size.

  and it's always evenly spaced between 0 Hz

Yes, FFT is always evenly spaced, and always starts at 0 Hz.  It doesn't implement Zoom FFT (since it's meant primarily for reconstruction with ISTFT, and Zoom FFT would throw away the required information).

  and half the audio sampling frequency.

Yes, your input is real, so the negative frequencies are redundant, so they are not returned.   nperseg is 100, so nfft is 100, so it returns 100//2+1=51 frequencies (including 0 and Nyquist).

If you give it a complex input, then len(f) will equal nperseg because it returns both the positive and negative frequencies.

  What if I want to have a hop size of 706 samples and 149 measured frequencies between 50 kHz and 65 kHz?

It doesn't support Zoom FFT, so, assuming your sample rate is 65*2 = 130 kHz, you'd need to calculate how many samples you want across the entire bandwidth (= 1292?) and enter that as nfft, and specify noverlap = nperseg - hop, and then ignore the frequencies below 50 kHz that are returned.

Note that specifying nfft > nperseg returns redundant interpolated information that isn't necessary for reconstruction.


In the stopband(s), Eqs $(1)$ and $(2)$ are equivalent, because in the stopband the desired response equals zero: $D(\omega)=0$. The reason why you might want to split the error in passband and stopband error is to apply different weights, as shown in your Eq. $(3)$.

The function freqz computes the frequency response of a discrete-time filter on a grid of frequencies, given the filter coefficients. I'm not sure what you mean by "achieve the transpose operation".


In the paper Haar-Based Stability Analysis of LPV Systems, the Haar wavelet transform theory have been used to design linear matrix inequalities (LMIs) to analyze the stability of LPV systems.

As the resolution level of the Haar wavelet increases, the number of variables and rows of the designed LMI increases and the feasibility of this LMI becomes a less conservative condition for the stability of the LPV system. Although the stability of LPV systems using LMIs that become less conservatives with the increase of variables and/or rows have already been proposed in the literature, the Haar-based approach can handle a larger class of parametric dependencies as well as non-convex parametric domains.


mathematical equation for Power Line communication channel

Agreed to others. You need to narrow down your interest to a frequency band as PLC can be achieved as narrow-band or wide-band. Your channel and noise model would depend on the frequency band you are interested in. For instance, if you are interested in broadband communication (e.g. OFDM), then you need to take transmission/reflection coefficients of the network of interest / topology, attenuation factor (which is a function of distance and frequency) and delay of all the paths in the receive path into account. Of course, the channel is also time-varying and the noise is always assumed to be non-Gaussian.

A pretty common channel model assumed for broadband communication would be as the following:

\begin{equation} \label{eqn_PLN_channel}
H(f)=\sum^L_{i=1}\left\{\prod_{k=1}^O \Gamma_{ik} \prod_{m=1}^P
T_{im} A(f, d_i) e^{-j2\pi f\tau_i}\right\}
\end{equation}
where $\Gamma_{ik}$ and $T_{im}$ denote the reflection coefficient
and transmission coefficient, respectively, $A(f, d_i)$ is the
distance and frequency dependent attenuation and $\tau_i$ represents
the delay of $i$th path out of all the receive paths ($L$) at the point of signal reception.

Derivation of this equation is simple and stems from what's physically happening in the network. Though, note that the channel depends both on frequency and time.

On the other hand, usually there are no closed-form estimates of the noise and the noise is numerically estimated based on a model such as Nakagami-m or Generalized Gaussian Distribution.


What is facebook's image compression algorithm?

If you have access to the bitstream, you should be able to figure out which codec they are using.
Assuming that it is good old jpeg, possibly with a tweaked encoder implementation, I would parse the stream using JPEGSnoop to learn as much as possible. It can even try to match «known»
encoders:
https://www.impulseadventure.com/photo/jpeg-snoop.html
If you can upload uncompressed bmp images to facebook and observe the first-generation compressed jpeg, you might learn even more as you have a reference.
Reverse engineering an encoder black box like this may be tedious unless they went with a plain vanilla implementation.
JPEG is actually a flexible standard, and «quality=85» is a simplified way of expressing how 8x8 DCT coefficients of luma (and chroma) are quantified. Sniffing the stream will allow you to read out those 2x64 quantization scalings, and also (if any) chroma subsampling, and (unlikely) non-bt601 color primary matrixing. An encoder is allowed to set any (small) dct coeff to zero if it likes (in the interest of smaller files), besides the regular quantization. Also, the encoder may use non-standard huffman tables but that should affect the bitstream not pixels.


Efficient correlation of a low duty cycle training sequence

One way to approach this would be to have a bank of $L$ subcorrelators, each of which use the corresponding $N$-sample long template. Then, in order to get the effect of correlating against the longer gated sequence, one need simply to sum the outputs of the subcorrelators with the appropriate delays applied. This is equivalent to breaking the long matched filter into a parallel formulation with $L$ subfilters.

For instance, assume the output of subcorrelator $i$ is $c_i[n]$, and that the notional output of the full $ML$-sample corelator is $c[n]$. The full correlator output is just equal to the sum of the shorter subcorrelators, after taking into account the $M$-sample delay between the expected position of each template in the input signal:

$$
c[n] = \sum_{i=0}^{L-1} c_i[n-(L-1-i)M]
$$

Using this structure, for each output sample, you effectively implement an $LN$-tap FIR filter instead of an $LM$-tap one. With that said, if any of these numbers are large (greater than 50-60 or so), it's likely more efficient to use FFT-based fast convolution techniques instead. This method is compatible with that; you would simply use fast convolution to calculate each of the short subcorrelator outputs, then perform the delayed sum as described above. I've implemented this very technique before with good success.

There are a couple other things to think about with this approach:

The OP suggests that some of the $L$ frames may not have distinct patterns from each other. If this is the case, you can achieve some complexity reduction by only implementing the number of unique subpatterns, then tapping the outputs from each subcorrelator appropriately.
This sort of technique can run into issues with scalloping loss if the overall length of the correlator ($LM$ samples) is large with respect to any carrier frequency uncertainty on the input signal. This is inherent with any kind of coherent correlator structure; the longer you correlate (thus giving you more potential gain), the more frequency selective the correlator is (the more loss you have for a given frequency offset in your signal).


Ratio of expected values of squares of errors in power quantities, in dB

Your question raises many very present concerns about measuring differences, and optimizing, in signal/image processing. Measuring can help:

to compare different outcomes from two different processes,
to optimize processing per se (for denoising, restoration).

Most use cost functions related to  $p$-norms (when $p\ge 1$) or quasi-norms (when $0<p\le 1$): $$\mathcal{C}_p(a,\hat{a})=\left(\sum_k|a_k-\hat{a}_k|^p\right)^{1/p}\,.$$ Those functions are symmetric, and are consistent with the physical units of $a$. The behavior of $\mathcal{C}_p$ depends on power $p$: classical norm inequalities tells you that $\mathcal{C}_p(a,\hat{a})\ge \mathcal{C}_q(a,\hat{a})$ when $p\le q$. Those (quasi-)norms are $1$-homogeneous. Using different $p$ or $q$ can be useful, especially in a sparsity context, see for instance references in: Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed \ell_1/\ell_2 Regularization. Here is a rendition of classical penalties:

Thus, quasi-norm or norm ratios are used to compare estimates, in a unitless fashion. The $\ell_p$ signal-to-noise ratios are typically of the shape:

$$\frac{\mathcal{C}_p(a,0)}{\mathcal{C}_p(a,\hat{a})}\,.$$

They are unit-less (ratios of same units). To better express orders of magnitude, one usually take logarithms. The logarithms can also been interpreted as unit-less measure. Interestingly, either $\log$ or $\ell_p/\ell_q$ ratios (like $\ell_1/\ell_2$) have been used a non-convex penalties in signal restoration.

Your functions with power $\alpha$ can be written as:

$$B\log_b \frac{\mathcal{C}_p(a^\alpha,\hat{a}_2^\alpha)}{\mathcal{C}_p(a^\alpha,\hat{a}_1^\alpha)}\,.$$

Now, those logarithms of ratios are unit-less, and are usually employed in a relative way, for comparison only. With $B>0$, $b>1$, $t \to B \log_b (t)$ is a monotonic increasing function. So, for value comparison only, we don't really care about the values for $b$ and $B$. However, for consistency in dB scale (as Olli Niemitalo already answered properly), signal-to-noise ratios tend to be computed in $10\log_{10}$ of energy measures, so homogeneously:

$ 20\log_{10}(|a|)$
$ 10\log_{10}(|a|^2)$
$ 5\log_{10}(|a|^4)$

make sense. However, it is not easy to compare them objectively, but they can be plotted on similar graph.


PSD using MATLAB's $\tt pwelch$ algorithm window length and measure units

Regarding your questions:

  Is there a good rule of thumb for determining the length $L$ of the
  sections (second variable for the function)?

$L$ (or window argument) specifies a window function to apply, prior to FFT computation of the signal. If $L$ is a vector, that vector is considered as window coefficients and is applied on (multiplied to) signal, for example:

Pxx=pwelch(x,[0 0 1 2 3 3 2 1 0 ]); % Means compute spectrum of x.*[0 0 1 2 3 3 2 1 0]

if it is an integer, then a Ham. or Han. window with length of $L$ is generated and then this Hamming window is applied to the signal:

Pxx=pwelch(x,10); % Means compute spectrum of x.*(10 point hamming window)

pwelch is based on FFT, and since the spectrum (DFT) of segmented signals shows frequency leakage, I think it would depend on how much resolution you you want and how much frequency leakage, you can tolerate, but generally The more longer the window for frequency static signals is the more leakage is reduced.
See this great tutorial for more detail:

http://www.fhnw.ch/technik/ime/publikationen/2012/how-to-use-the-fft-and-matlab2019s-pwelch-function-for-signal-and-noise-simulations-and-measurements

A side note for example codes above is that if the length of your signal x is greater than specified length for the window, then your signal is segmented to vectors with length of $L$

  assuming my input is a voltage signal , what are the measure units I get in the output?

Since your signal is pure voltage (unknown impedance and hence no current (A) is associated), it is not actual "power" or "energy", however you might use $\mathrm{V}^2/\mathrm{H}z$ for PSD.
This link might shed more light on it:
https://en.wikipedia.org/wiki/Spectral_density#Electrical_engineering
In case of logarithmic PSD, you've mentioned it is generally regarded as $\mathrm{dB}/\mathrm{Hz}$.


Why do we use "error quaternions" for control of attitude?

A possible explanation is that a set of 4 quaternions is not independent, so by converting them to error quaternions, we get 3 independent values that we can use to calculate 3 independent control actions (one for each axis).


Training and testing set for speech recognition

There are many open source datasets for speech available for free. You can find most at http://www.openslr.org

Beside that you can check http://voxforge.org


What would be the outcome of the convolution of two images?

You can use MATLAB to get the 2D convolution of these images. Before you do so, you'll want to convert image 2 to grayscale and convert both to double.

The result has a very large range, so I normalized the second image.

lena = im2double(imread('lena.jpg'));
color = im2double(rgb2gray(imread('color.jpg')));
out = conv2(lena, color./sum(sum(color)));

If you happen to be interested in blending images together, you might want to check out this paper.


How does resizing an image affect the intrinsic camera matrix?

Note: That depends on what coordinates you use in the resized image. I am assuming that you are using zero-based system (like C, unlike Matlab) and 0 is transformed to 0. Also, I am assuming that you have no skew between coordinates. If you do have a skew, it should be multiplied as well

Short answer: Assuming that you are using a coordinate system in which $u' = \frac{u}{2} , v' = \frac{v}{2}$, yes, you should multiply $a_x,a_y,u_0,v_0$ by 0.5.

Detailed answer The function that converts a point $P$ in world coordinates to camera coordinates  $(x,y,z,1)->(u,v,S)$ is:

$$ \left( \begin{array}{ccc}
a_x & 0 & u_0 \\
0 & a_y & v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

Where $(u,v,S)->(u/S,v/S,1)$, since the coordinates are homogenous.

In short this can be written as
$ u= \frac{m_1 P}{m_3 P} , v = \frac{m_2 P}{m_3 P}$
where $M$ is the product of the two matrixes mentioned above, and $m_i$ is the i'th row of the matrix $M$. (The product is scalar product).

Re-sizing the image can be thought of:

$$ u'=u/2, v'=v/2 $$

Thus

$$ u' = (1/2) \frac {M_1 P} {M_3 P}  \\
v' = (1/2) \frac {M_2 P} {M_3 P}
$$

Converting back to matrix form gives us:

$$
\left( \begin{array}{ccc}
0.5 & 0 & 0 \\
0 & 0.5 & 0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
a_x & 0 & u_0 \\
0 & a_y & v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

Which is equal to

$$ \left( \begin{array}{ccc}
0.5 a_x & 0 & 0.5 u_0 \\
0 & 0.5 a_y & 0.5 v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

For additional information, refer to Forsyth, chapter 3 - Geometric camera calibration.


When does the Amplitude Spectrum equal the Magnitude Spectrum?

You are right that also for the second signal you have $M(\omega)\neq A(\omega)$, but the point here probably is that in the fundamental interval $\omega\in (-\pi,\pi)$ the two are identical. It is only at $\omega=\pm\pi$ that the phase $\phi_M(\omega)$ jumps and that $A(\omega)$ changes its sign. So what they probably mean is

$$M(\omega)=A(\omega)\quad\text{and}\quad\phi_M(\omega)=\phi_A(\omega)\quad\bbox[5px,border:2px solid red]{\text{for }\omega\in(-\pi,\pi)}$$

which is correct (and which is not the case for the first signal).

Remember that such quizzes are made by fallible humans.


Generating a Sine Signal (for example a voltage sine wave) with changing frequency and changing sample rate (time)?

If you are doing this in realtime, then your generated signal WILL ALWAYS lag behind the original.
Your sampling rate will not change because of anything your calculation does.
You have a problem with discontinuities when the frequency of your measured signal changes - you are abruptly changing from the old to the new frequency.  This will cause skips in the generated signal.
You will not always be able to exactly match the real signal - the zero crossing may take place between two samples.  This will cause your calculated frequency to be wrong.
The calculated frequency won't necessarily fit well with your sampling rate, so that it won't end exactly on a sample - this will cause another skip in your generated signal.

You will have to gradually change f(t) in order to avoid the abrupt changes.  So, you measure f one time and it is (say) 1000Hz, then the next time it is 1100hz.  You cannot simply change from one to the other.  You will have to "walk" the value for f from 1000Hz to 1100Hz over the length of time it takes for one half cycle of 1100Hz.

That should improve things a good bit.


This is a very complicated topic. A good overview on how to do this in MATLAB is here:
https://ccrma.stanford.edu/~jos/vguitar/Fitting_Filters_Matlab.html

An FIR fit can typically done by interpolating up to a full FFT grid and adding either a minimum or a zero phase and then doing an inverse FFT to get the impulse response. You can then try different filter lengths to determine which lengths is required to create an acceptable fit.

IIR fitting is much more tricky. The above example uses MATLAB's invfreqz() which in some cases gives good results and in other only works very poorly. In general the methods are iterative, i.e. there are many design cycles and there are small modifications between design steps that are supposed to "improve" things. These often use "gradient" methods, i.e. you need to evaluate the derivatives of the error function with respect to zeros and poles (or whatever other representation you may choose). They also involve least square error methods. A lot of care has to be taken to correctly formulate the error criteria in terms of frequency spacing, weight functions, relative vs. absolute error etc. The formulation of the error function is a direct representation of what the specific application can tolerate (or not).

EDIT based on OPs question edit:
If you want to do frequency domain filtering, (i.e. FFT->mpy with Filter->inverse FFT) then you need to get all the framing, zero-padding and overlap handling correct. This due to the difference of circular vs. linear convolution and that the FFT "as is" assumes by definition periodic signals. Google "overlap-add" for more info.

You still need to decide on the phase for your filter. Leaving the phase zero creates a non-causal filter that will create lots of wrap-around problem. At the very minimum you need to shift the impulse response so it's causal.


Remove Shadows from Contour of an Object (Grapes)

This is what your code produces (using colormap(hot) with white boundaries that are obtained from IM2):

I think the shadows are nicely treated, why are you thinking that your code doesn't notice between shadow to grapes?

This is the code used for producing this image

bw=IM2<1;
[B,L] = bwboundaries(bw,'noholes'); % fill any holes, so that regionprops can be used to estimate
imagesc(GRAY); hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 1); hold on
end
colormap(hot(256))


How does the power spectrum of the amplitude envelope depend on its underlying carrier frequency?

Your intution is right, the carrier frequency has to be at least twice the highest frequency in your signal. Practically, much higher carrier frequencies are chosen, as you would need ideal or nearly ideal lowpass filters for carrier frequencies this low.


To elaborate on Dilip's answer (which is perfectly correct, though in practice the Ziggurat method is much more computationally efficient than Box-Mueller):

One key ingredient missing from your reasoning is whether you want your samples to be independent. It is not clear from your question, but this is the most common situation...

If you want your samples to be independent, then you will have to accept that their empirical mean and variance will not exactly be 0 and 1 - what is known is that the larger the sample, the closer the empirical values will be to 0, 1 - this is absolutely fine - this is a "feature" of the normal distribution rather than a "bug". It is a bad idea to try to select the samples that look more "random" or "well-behaved"; because this very selection process is making them less random! For example, if you come up with a process that always results in 2048-sample blocks with exactly $\mu = 0$ and $\sigma = 1$, it would mean that there are only 2046 possible degrees of freedom in your generation process, and thus that your samples are not independent (given the 2046 first samples, you can guess the values of the 2047th and 2048th using the $\mu = 0$ and $\sigma = 1$ equation... so these last two values are not random!).


Generating data with Gaussian distribution

Let $\Delta = [x_1, \dots, x_N]$ be your desired vector where $x_i$ are iid $\mathcal{N}(0,\sigma^2)$. Do you expect that $\mathrm{Pr} \left\lbrace ||\Delta||_2^2 < 0.5\right\rbrace = 0.9$ ?

By definition $y = ||\Delta||_2^2 / \sigma^2 = \sum_{i=1}^N \frac{x_i^2}{\sigma^2} \sim \chi^2_N(y)$ is Chi squared distribution of degree of freedom $N$.

Thus $z = ||\Delta||_2^2 = \sigma^2y \sim \frac{1}{\sigma^2} \chi^2_N(\frac{z}{\sigma^2})$

Denote the CDF $\mathrm{F}(z)$, you want something like $\mathrm{F}(z = 0.5) = 0.9$ with $\sigma^2$ is unknown variable.

A closed form expression of $\mathrm{F}(z)$ is available at wikipedia Chi squared distribution.

In MATLAB you can use the two functions chi2pdf() and chi2cdf() for numerical brute force solution for $\sigma^2$.

And after having the value of $\sigma^2$, just use the function randn() to generate your vector $\Delta$. For example with size N and $\sigma^2$ = sigma2, delta = randn(1,N) * sqrt(sigma2).


understanding time constant meaning in signal processing

There is no general formula to calculate time constants for all signals. Time constants are defined for some signals. For your exponential signal

\begin{equation}
e^{-a \cdot t} = e^{-t / c}
\end{equation}

it just means that always after the time $c$ has passed, the value of the function has been reduced by $1 / e$. A more common interpretation is half-life defined as

$$t_{1/2} = \frac{ln(2)}{a}$$

which is the time after the function has gone from 1 (at time zero) to 1/2 (at time $t_{1/2}$) and so on.

However, there is no general time constant that can be applied to any signal, because, as you already mentioned, it is unclear what it represents.


How to remove noise using mean and median filters?

You do not need loops. This implements a moving average filter, a median filter, and a median filter followed by a moving average:

%x = noise1;
x = randn(1,2000)+1;
lengthMean = 5;
lengthMedian = 7;
xFiltMedian = medfilt1(x,lengthMedian);
xFiltMean = filter(ones(lengthMean,1)/lengthMean,1,x);
xFiltMedianMean = filter(ones(lengthMean,1)/lengthMean,1,xFiltMedian);
figure(1);clf;hold on
plot(x(500:1500),'b');
plot(xFiltMedian(500:1500),'r');
plot(xFiltMean(500:1500),'g');
plot(xFiltMedianMean(500:1500),'c');

As said by @jojek the median filter is often more efficient than an averaging filter for spikes and outliers. And if you choose to apply both filters, it is generally more useful to apply the  median first (as done in the code above).

The choice of filter length should normally be related to properties of your signal and the noise (spectrum, distribution), and to the goal  you are assigned with. But the length of the signal is usually not a criterion.  For symmetry  properties, odd-length filters are practical.

Without further details, you can start with longer filters (e.g. lengthMedian= 37), see how they perform (how they reduce the  noise but also degrade the signal's properties), and reduce the length progressively to the equilibrium point, where shorter filters do not yield better results.

However, the noise in the  valley suggests it is not stationary, so you may need more adaptive filters (e.g. with a varying length), or a little more involved like a weigthed median filter that limits the location and flattening issues you get with long median filters. It is discussed on SE Cross validated.


Determination of periodicity in data and finding mean

To determine periodicity, do an autocorrelation. The autocorrelation will have peaks at time offsets consistent with the periodicity.

You mentioned using the Fourier Transform, if convenient you can do a circular autocorrelation using the FT as follows:

$$Corr = ifft(fft(x) fft(x)^*)$$

Where $fft(x)^*$ is the complex conjugate (if you signal is complex).


You can borrow from techniques used by the image stitching community, where they have to deal with registering translated and rotated images. Your mistake was that you calculated the correlation before finding the optimal translation. The Wikipedia article should give you enough information to help yourself.

Mathematica Example

First, I split the picture above into frame1.png and frame2.png. Then I shifted the second one to align with the first, and trimmed both to retain their intersection. Finally, I compute the correlation coefficient by flattening the matrix into a vector.

img1 = ColorConvert[Import["frame1.png"], "Grayscale"]
img2 = ColorConvert[Import["frame2.png"], "Grayscale"]
im2t = ImageAlign[img1, img2, "Transformation" -> "Translation"]
imgtr = FindGeometricTransform[img1, img2, "Transformation" -> "Translation"]
{trx, try} = Part[TransformationMatrix@Last@imgtr, {1, 2}, 3]
{trim1, trim2} = ImageTake[#, {1, Floor@try}, {1, Floor@trx}] & /@ {img1, im2t}
{cimg1, cimg2} = Flatten@First@ImageData[#, "Interleaving" -> False] & /@ {trim1, trim2};
Correlation[cimg1, cimg2]

The result I got was 0.967834, indicating a close match. If necessary I would have considered scaling and rotation to better align the images but it was not necessary here.

Here is a MATLAB example.


Relating phase shift to the ROL/ROR operations

If you mean phase shift of a sinusoid, then yes, a phase shift is equivalent to a translation.


Fast fourier transform for an audio file

Your question is so "basic", so "fundamental", that I may be misunderstanding its true meaning. In any case, hopefully you have some software routine available to you that performs FFTs and will accept your 'int input[N]' sequence as an input. I have no idea what software 'package', or software language, you're using but if variable 'x' is your input sequence then there should be some available command (or subroutine) like

Spec = fft(x)

you can use to compute your FFT spectral results. And there should be no need for you be concerned about defining your input sequence as complex-valued. Of course the length of your input sequence (number of samples) must be an integer power of two. Your FFT results will automatically be in the form of a sequence of complex numbers, so there's no need to worry about that.


Why does FFT output need to conjugate processing?

The real part represents cosine correlation, and the imaginary coefficient represent the sine correlation.  Cosines are even or symmetric functions in the DFT aperture, and sines are odd functions (or anti-symmetric) in the DFT aperture.  So you can only ignore the imaginary coefficients for exactly symmetric input.

Conjugation and multiplication is “just” a mathematical trick (identity) to compute sin^2 + cos^2 which is the hypotenuse squared (Pythagoras).


How can I find the actual frequency given a normalized frequency?

[Y,FS,NBITS]=wavread(FILE) returns the sample rate (FS) in Hz and the number of bits per sample (NBITS) used to encode the data in the music file.

After that, [S,F] = pmusic(Y,4,F,FS) to compute the pseudospectrum at frequencies F in Hz.


There's no difference in formula:
Try setting the $C=\text{diag}(\sigma_1^2, \sigma_2^2,\ldots,\sigma_n^2)$ in $(1)$ to see how $(1)$ simplifies to $(2)$. So, your formula $(2)$ is just a special case of $(1)$.
In the diagonal case, the elements of the multivariate distribution are jointly Gaussian, which is much stronger than them all just being Gaussian distributed.
Joint Gaussian distribution has a name of its own (unlike other multivariate distributions with identically distributed elements) because "joint gaussian multivariate random varaible" has the very interesting feature that simple MVG doesn't have: that uncorrelatedness implies independence. This is actually something that makes it useful in communications engineering (because rotation and normalization doesn't change joint gaussianness and doesn't introduce interdependencies between components - so, in an IQ receiver, the noise on I and Q are still independent, even after phase correction, for example).


How does harmonic distortion occur in ADCs?

There are many distortion mechanisms in ADC but the most important ones are:

"static" distortion due to INL
Track-and-hold distortion

The first one is caused by the fact that the voltage-in-code-out relationship is not perfectly linear (INL/DNL).

An approximation of harmonic distortion based on DNL is given by (Lee, JSSCC 12/2017):

$$
\mathrm{HD3} = -20\log\left( \frac{4}{3\sqrt{3}} \frac{2^B}{\mathrm{INL}_{\mathrm{max}}} \right)
$$

The second one is caused by track nonlinearity. The sampling switch is for example a weakly nonlinear CMOS switch. Its on-resistance is given by (using the simple square law model):

$$
R_{on} = \frac{1}{\mu C_{ox} \frac{W}{L} (V_{DD}-V_{in}-V_t)}
$$

From which is can be seen that the resistance depends on the input signal (see for example Razavi, Data Conversion System Design, page 16).

This distortion can be approximated as

$$
\mathrm{HD3} = 20\log\left(\frac{1}{4} \left(\frac{A}{V_{GS}-V_t}\right)^2 2\pi f_{in} R_{on} C \right)
$$

For a simple CMOS switch this is by far the dominant distortion source and is practically overcome by using a bootstrap sampling switch. However, also this technique has limitations.


If video uses color bars to test and align color values among TV sets, what could be the equivalent in Audio?

If you want the AI to relate to how we perceive sound, maybe something like

a linear chirp from 20Hz-20kHz
white noise in the same bandwidth
a single sine wave at 1kHz (roughly our greatest sensitivity)
a harmonic series from 55hz doubling up to 14.08kHz (all A notes to illustrate octaves)
other harmonic series and/or progression to illustrate important musical intervals (major and minor chords, pentatonic scale)
sound level variation of one or more of the above between 0dB SPL and some appropriate limit to establish dynamics

In terms of storage, a set of vectors is probably the best place to start. You can make one long vector with your test pattern and load it into and out of a WAV file I suppose.


Basic questions related to Kalman filter

A state has memory. The value of a state of a causal system reflects the values of all the previous states and inputs. A linear combination of all previous states and inputs, is a linear state equation.

A measurement only reflects the value of the current state and current inputs.  A measurement needs to be noisy. A Kalman Filter estimates the current state from the noisy measurement.

The presence of measurement noise is why it is called a filter. If there is no measurement noise, one uses a state observer, not a Kalman Filter. The KF equations don’t work without measurement noise.

If the measurement noise is a linear combination of the current state and input, which includes noise, the measurement is linear. When the state and measurements are linear, the KF is a linear KF.

Essentially a state has memory and a measurement only depends on the current state. We can write the state evolution as a recursive equation because it is a Markov process.

The dimensional issues are related to what older books called the controlability and observably matricies. More recent books introduce Grammians.

There are variations of the KF that include nonlinear states and  nonlinear measurements. The linear KF satisfies a number of optimality criteria if the matrices and initial conditions are accurate.  The other types aren’t optimal but typically good if the approximations are good.

One does not add error unless you are doing a simulation. you need to know the covariance to actually filter.

The KF can be derived a few different ways, (least squares, orthogonality, bayesian) and the solutions can be written out differently so books often don't exactly agree.  I prefer the innovations form of the KF because the innovations are an online diagnostic.  You can test if you have enough states or too many states.

The questions related to the actual state and measurement equations are really independent of the KF.  A book on linear systems using state variables is where you learn to write those, but many KF books have a chapter or 2 on the topic.  Usually, how to convert from continuous time differential equations to a set of difference equations.  Writing the state equations is really a prerequisite.  Most of the people who ask for help tend to fall in the category of not understanding state variables.

If you have problems understanding your state equations, you should include the specific equations.  No one is watching what you see on the internet. This isn't pornhub.

Kalman Filters tend to work even if you don't know those parameters well.  People tend to spend a lot of time tweaking their filters so expect to tweak a lot.  The theory can be intimidating and one can often build adequate filters without a deep understanding. If it works, it works.  If it doesn't work and you think it should, the theory becomes important.


Get problem to understand "Topological Structural Analysis of Digitized Binary Images by Border Following" Suzuki Abe

In the paper 4-(8-) refers to the 4 and 8 neighbor connectivity between the two pixels and in the algorithm NBD is the newest border found during the raster scan of the binary image and LNBD refers to the last newest border, NBD and LNBD is stored to obtain the surroundness relation between the borders.


First of all those sequences you provided must have a time index reference as well. I assume them to be causal sequences and $n=0$ to be the first sample.

Then you can use different approaches to find the impulse reponse for example using the Z-transforms provide: $$ H(z) = \frac{Y(z)}{X(z)} $$
Then you have to convert $H(z)$ back into time to find the impulse response $h[n]$.

But for this simple problem, I would instead use a time-domain approach directly and compute the samples of $h[n]$ by the following procedure:

Given $x[n] = -\delta[n] + \delta[n-1]$ of length $N_x=2$

and $y[n] = -\delta[n-1] -\delta[n-2] + 2\delta[n-2]$ of length $N_y=4$

(note that the first sample $y[0]$ of $y[n]$ is zero but here we consider $y[0]$ to be in the support-domain of $y[n]$ no matter what its value is, the result won't change for this problem)

1- Assuming an FIR $h[n]$ of length $L = N_y - N_x + 1 = 4-2+1 = 3$ samples.

2- Now propose the causal $h[n]$ to be $$h[n] = a\delta[n] + b\delta[n-1] + c\delta[n-2]$$

3- Now apply the convolution as the sums of shifted $h[n]$ rows as follows:

$$ y[n] = h[n]\star x[n] = x[0]h[n] + x[1]h[n-1] $$
$ $

$~~~~~~~~~~~~~~~ -a ~~ -b ~~ -c$

$~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~a ~~~~~~~b ~~~~~c$

$ + ~~~ ------------$

$y[n]=~ -a ~~(a-b) ~~ (b-c) ~~c$

$y[n]=~ -0 ~~~~(-1) ~~~~ (-1) ~~~~~2$

from which you can easily solve $a,b,c$ to be

$c = 2$ , $b-c=-1 => b = 1$ and $a=0$, therefore the impulse reponse is found to be $$ h[n] = 0\delta[n] + 1\delta[n-1] + 2\delta[n-2]$$


Signal processing - why compute imaginary part?

One use of Hilbert Transform is to recover the amplitude envelope of a signal.

Here is a practical example: Extracting Binary Magnetic-Strip Card Data from raw WAV

Often the shortest distance between two points is through complex analysis.  Even if the start and destination don't appear to involve complex numbers.

Regarding "What is the phase shift of 90°?": If it is single cisoid, i.e. a single point making a circular orbit around the origin on the complex plane, then applying a phase shift of 90° obviously sets it backwards or forwards by a quarter-cycle.  Now if you consider an arbitrary waveform to be a sum of cisoids with different amplitude, phase and frequency, (Fourier Theory shows that you can generally achieve this), you can just set each one back by a quarter cycle, and recombine. That's one way to accomplish the Hilbert Transform.

http://bl.ocks.org/jinroh/7524988 <-- this gives a visual example of combining cisoids into the target waveform.


What is the best method for measuring of room reverberation time?

Well, there is plenty of ways for approaching this problem. Which one is best depends on room itself, your resources and pursued accuracy. I think that best way to start is to look into ISO 3382-2:2008 standard, you will find lot of information. You should remember that post processing of your signals is very important here.

Nevertheless you should consider following methods:

Hand clap, balloon, starter gun. Cheap and most intuitive. You are recording the "impulse" itself from which you are then retrieving the reverberation time and other parameters. Basically almost every other method is relaying on getting the impulse response of your room. After that you filter it in frequency bands of interest and find the slope. Order is not random here. Hand clap doesn't have wide frequency response in low frequencies, balloon is slightly better, and starter gun will give best results.
Intermittent noise. You play the band-pass filtered pink noise (usually in octaves) and register it's level in dB scale. You turn off playback and record the decay. From the slope in dB scale you find the time at which energy curve will cross -60 dB from steady state. This is your reverberation time for this given octave band. It is important to use omnidirectional source and microphone. It is also good if you repeat results for different source-microphone arrangement, i.e. for 2 positions of your source and 6 positions of microphone.
Sweep sine. For that you presumably need a piece of software. You can either write your own (here is how to do it), use a free one (such us REW) or even buy some professional tool (i.e. Dirac, EASERA, etc.). Similarly to above you play signal from speaker, record for few positions and then analyse using favourite application - as simple as that. This is my method of choice.
Maximum Length Sequence. This is white-noise-like signal that correlates nicely with itself and can be used for extraction of impulse response. Personally I do not recommend it as sweep sine seems to work better for me in every case. Although some people like to use it for measurement of speaker cabinets.


What is a good way to locate butterfly in an image

The butterflies are really different and the backgrounds as well, then it makes it really difficult. I would try one of these ideas:

Train the Viola & Jones algorithm on your butterflies. It's a cascade of weak classifiers usually used for face detection (algorithm implemented in all the cameras), but it can be trained for other purposed.
The background seems to be pretty different of the butterflies. I would try to apply the thresholding by histereis according to the values of the pixels on the borders. Like that you select all the similar pixels that touch the border, and it will reduced a lot the region of interest.


Visualize visual words in bag of words model

Each SIFT descriptor corresponds to a region of the image. You take these from a bunch of images and group them into some number of clusters. I think what he's showing in the slide there is just a few samples from each cluster where he chose human-meaningful names for the clusters after the fact.


In this document you can find the coefficients of a seven-term Blackman Harris window. Ignoring the bizarre notation, it seems like the window is defined by

$$w[n]=\sum_{k=0}^6a_k\cos\left(\frac{2\pi kn}{N}\right),\quad n=0,1,\ldots,N-1\tag{1}$$

where $N$ is the window length and the coefficients $a_k$ are given by

$$\begin{align}
a_0&=0.27105140069342\\
a_1&= - 0.43329793923448\\a_2&=0.21812299954311\\a_3&= - 0.06592544638803\\a_4&=  0.01081174209837\\a_5&= - 0.00077658482522\\a_6&=  0.00001388721735\end{align}
$$

Unfortunately, they do not mention where those coefficients come from. Somebody must have run an optimization routine to minimize the sidelobe level for a given mainlobe width.

EDIT:

I think I found the relevant paper here. This is a corresponding Matlab function.


Please have a look at Page 5 of the DataSheet here.

It says the signal Bandwidth is 8500Hz only when output data rate is 25600Hz and Low Pass Filter setting = 1. And that the Bandwidth changes based on Output Data Rate.

So, I think you do not have to worry about sampling rate, since ODR(Output Data Rate) is 25600Hz and that means 16-bits(2Bytes) readings coming at this rate.

So, 1 second of recording will require : 1*25600*2Bytes = 51200/1024 KB = 50KB, of memory space. And, I am not sure but I think this is just for X-Axis measurements by the accelerometer.


The 8kHZ signal causes aliasing as we should sample at 2fmax = 2*8kHz = 16kHz (normally we should sample faster than 2fmax). It is also strange that you have 2 signal @4kHz and 8kHz but you want to use a LP at 10kHz.


Eventhough Laurent has given a broader sense of the answer, let me put here the communications theory sense ot it.

The concept of instantaneous freqency emerges when you consider Frequency Modulation or Phase Modulation systems, where the message is embedded into the change of the frequency or phase of a carrier signal.

This carrier is typically a single cosine wave or a complex exponential equivalent which posses a single well defined angle, whose time derivative gives you the instantaneous frequeny; hence the message. That's why they are also called angle modulation or exponential modulation systems.

The understanding is that there must exist a single carrier to represent the modulated message. Your example includes a general sum of three complex exponentials (unless you speciy those $\alpha$ pars) whose simplification into a single carrier seems not possible on the surface. Then there won't be a unique single instantaneous frequency definition for $u(t)$.

But you can consider three different instantaneous frequencies per complex exponential for your signal $u(t)$. Indeed this would be the case when multiple FM modulated signals (at different frequencies) were existing simultaneously as they are transmitted at different channels.


Optimal amplitude of an $m$-bit sinusoid

The phase of the sinusoid does not matter: A phase shift of a sinusoid is equivalent to shifting it in time, which results in a time shift of both the quantized sinusoid and the quantization error. The power spectrum is invariant to time shifts. We choose to work with sinusoid $A\cos(x)$.

Optimality by def. 1

Equivalently to maximizing signal-to-noise ratio (SNR), we can minimize the root mean square quantization error relative to the root mean square $\sqrt{1/2}A$ of the sinusoid, $\mathrm{RRMSE}$ $=$ $\sqrt{1/\mathrm{SNR}}$. It is enough to do the analysis in the first quarter-period of the cosine wave, because the rest of both the cosine wave and the quantization error are identical to their first quarter-period up to a sign flip and/or time reversal. Let $x$ be in the first quarter of the cosine wave, $0 < x < \pi/2$. Along the lines of Eq. 2 of my answer to a related question, the quantized cosine wave attains an integer value $k \in 0\ldots\mathrm{round}(A)$ when:

$$\begin{align}0 < x < \mathrm{acos}\left(\frac{\mathrm{round}(A) - 0.5}{A}\right),\quad &\text{if }k = \mathrm{round}(A),\\
\mathrm{acos}\left(\frac{k + 0.5}{A}\right) < x < \mathrm{acos}\left(\frac{k - 0.5}{A}\right),\quad &\text{if }1 \le k \le \mathrm{round}(A) - 1,\\
\mathrm{acos}\left(\frac{0.5}{A}\right) < x <\frac{\pi}{2},\quad &\text{if }k = 0.\end{align}\tag{1}$$

Each value of $k$ gives to the relative mean square error $\mathrm{RelMSE}$ an additive contribution of:

$$\frac{2\,\mathrm{MSE}_k}{A^2} = \frac{2}{A^2}\frac{1}{\pi/2}\int_{x_0}^{x_1}\big(A\cos(x)-k\big)^2dx = \frac{2\sin(x_1)\cos(x_1) - 2\sin(x_0)\cos(x_0)}{\pi} + \frac{8k\big(\sin(x_0) - \sin(x_1)\big)}{\pi A} + \frac{4k^2(x_1 - x_0)}{\pi A^2} + \frac{2x_1 - 2x_0}{\pi},\tag{2}$$

where $x_0$ and $x_1$ denote, as defined separately for each $k$, the bounds $x_0 < x < x_1$ given by Eq. 1. The total relative mean square quantization error is then:

$$\begin{gather}\mathrm{RelMSE} = \frac{2\,\mathrm{MSE}}{A^2} = \sum_{k=0}^{\mathrm{round}(A)}\frac{2\,\mathrm{MSE}_k}{A^2}\\
= \frac{2\,\mathrm{asin}\left(\frac{1}{2A}\right)}{\pi} - \frac{\sqrt{4A^2 - 1}}{2\pi A^2} - \frac{(2A^2 + 4\,\mathrm{round}(A)^2)\,\mathrm{asin}\left(\frac{2\,\mathrm{round}(A) - 1}{2A}\right)}{\pi A^2} - \frac{(6\,\mathrm{round}(A) + 1)\sqrt{4A^2 - (2\,\mathrm{round}(A) - 1)^2} - 2\pi(A^2 + 2\,\mathrm{round}(A)^2)}{2\pi A^2}\\{+\,\frac{1}{\pi A^2}\sum_{k=1}^{\mathrm{round}(A)-1}}\left(\left(2A^2 + 4k^2\right)\bigg(\mathrm{asin}\left(\tfrac{2k + 1}{2A}\right) - \mathrm{asin}\left(\tfrac{2k - 1}{2A}\right)\bigg) + \frac{(6k - 1)\sqrt{4A^2 - (2k + 1)^2} - (6k + 1)\sqrt{4A^2 - (2k - 1)^2}}{2}\right),\quad A > 0.5.\end{gather}\tag{3}$$

For optimality by definition 1, the task is to find the value of $A$ that minimizes the relative root mean square quantization error RRMSE $= \sqrt{\mathrm{RelMSE}}$. under the constraint $A \le 2^{m-1}-0.5$. Eq. 3 can be evaluated in Python using the mpmath arbitrary precision library:

import mpmath as mp
def RelMSE(A):  # valid for A >= 0.5
    A = mp.mpf(A)
    return 2*mp.asin(1/(2*A))/mp.pi - mp.sqrt(4*A**2-1)/(2*mp.pi*A**2) - (2*A**2 + 4*mp.floor(A + 0.5)**2)*mp.asin((2*mp.floor(A + 0.5) - 1)/(2*A))/(mp.pi*A**2) - ((6*mp.floor(A + 0.5) + 1)*mp.sqrt(4*A**2 - (2*mp.floor(A + 0.5) - 1)**2) - 2*mp.pi*(A**2 + 2*mp.floor(A + 0.5)**2))/(2*mp.pi*A**2) + mp.nsum(lambda k: (2*A**2 + 4*k**2)*(mp.asin((2*k+1)/(2*A)) - mp.asin((2*k-1)/(2*A))) + ((6*k-1)*mp.sqrt(4*A**2 - (2*k + 1)**2) - (6*k + 1)*mp.sqrt(4*A**2 - (2*k - 1)**2))/2, [1, mp.floor(A + 0.5)-1])/(mp.pi*A**2)

RRMSE appears to have a local minimum between each pair of successive integer $A$ (Fig. 1).

Figure 1. RRMSE (blue solid line and blue squares) and its approximation $\sqrt{1/6}/A$ (orange dashed line) based on the variance $1/12$ of a uniform distribution of width 1, for various ranges of $A$.

A selection of optimal $A$ are given in Table 1 together with the resulting RRMSE, also for some other common values of $A$. At larger $m$, RRMSE reduction by the optimal choice will be marginal. The table entries, which show only digits that are equal between two computations using different precision settings, can be generated by the following Python script (continued), the execution of which took days:

def approx_optimal_A(m):
    m = mp.mpf(m)
    return 2**(m-1) - 1 + mp.mpf("0.156936321399") + mp.exp(-1.2749819017 - 0.3464088917*m)  # Eq. 8
    # return 2**(m-1) - 1  # This less informed guess gives identical results but slower convergence

def to_max_digits(f, prec_1, prec_2, max_digits):  # return the at most max_digits digits of function f that are agreed about by both precision settings
    prec = mp.mp.prec
    mp.mp.prec = prec_1
    y_prec_1 = f()
    mp.mp.prec = prec_2
    y_prec_2 = f()
    digits = max_digits
    while mp.nstr(y_prec_1, digits, strip_zeros=False) != mp.nstr(y_prec_2, digits, strip_zeros=False):  # Beware: a possible infinite loop
        digits -= 1
    return mp.nstr(y_prec_2, digits, strip_zeros=False)
    prec = mp.mp.prec

double_digits = 15  # Print at most this many digits
dB_digits = 9

for m in range(2, 25):
    optimal_A = to_max_digits(lambda: mp.findroot(lambda A: mp.diff(RelMSE, A), approx_optimal_A(m)), 80, 100, double_digits)
    RelMSE_optimal = to_max_digits(lambda: 10*mp.log10(RelMSE(mp.mpf(optimal_A))), 80, 100, dB_digits)
    RelMSE_1 = to_max_digits(lambda: 10*mp.log10(RelMSE(2**(m-1)-1)), 80, 100, dB_digits)
    RelMSE_2 = to_max_digits(lambda: 10*mp.log10(RelMSE(2**(m-1)-0.5)), 80, 100, dB_digits)
    print(str(m)+"&"+optimal_A+"&"+RelMSE_optimal+"&"+RelMSE_1+"&"+RelMSE_2+"\\\\")

Table 1. Optimal $A$ by definition 1 for different $m\le24$ and the resulting RRMSE, with the RRMSE for some common choices of $A$ listed for comparison. $m=1$ has been left out because it cannot be handled by Eq. 3 and because a single bit does not represent any positive numbers as a two's complement representation. Note that RRMSE in dB is converted to SNR in dB by a sign flip, because $\mathrm{SNR}$ $=$ $1/\mathrm{RelMSE}$ $=$ $1/\mathrm{RRMSE}^2$.

$$\begin{array}{l|l|lll}m&\text{optimal }A&\rlap{\text{RRMSE (dB)}}&&\\
&&A=\text{optimal}&2^{m-1}-1&2^{m-1}-0.5\\
\hline
2&1.26827949461530&-11.1128053&-8.92729805&-10.1764645\\
3&3.23800942121037&-18.8206937&-17.9588863&-17.8882004\\
4&7.21658597929407&-25.5167673&-25.0903048&-24.7375654\\
5&15.2007181331955&-31.8127094&-31.5775537&-31.2013629\\
6&31.1888756714257&-37.9364880&-37.7977883&-37.4726954\\
7&63.1800835394190&-43.9858910&-43.9001114&-43.6414308\\
8&127.173613625523&-50.0049518&-49.9500053&-49.7527817\\
9&255.168894736361&-56.0134743&-55.9773382&-55.8307200\\
10&511.165479188890&-62.0200022&-61.9957638&-61.8884975\\
11&1023.16302205377&-68.0278661&-68.0113689&-67.9337180\\
12&2047.16126264484&-74.0380590&-74.0267095&-73.9708981\\
13&4095.16000722516&-80.0505958&-80.0427267&-80.0028089\\
14&8191.15911371601&-86.0651409&-86.0596541&-86.0312009\\
15&16383.1584789666&-92.0812822&-92.0774410&-92.0572079\\
16&32767.1580286428&-98.0986407&-98.0959437&-98.0815798\\
17&65535.1577094659&-104.116904&-104.115007&-104.104821\\
18&131071.157483397&-110.135830&-110.134493&-110.127276\\
19&262143.157323352&-116.155234&-116.154291&-116.149181\\
20&524287.157210089&-122.174984&-122.174318&-122.170701\\
21&1048575.15712995&-128.194979&-128.194509&-128.191950\\
22&2097151.15707326&-134.215150&-134.214818&-134.213008\\
23&4194303.15703317&-140.235446&-140.23521&-140.233931\\
24&8388607.15700481&-146.255831&-146.2557&-146.25476\end{array}$$

Optimality by def. 2

A periodic function such as a quantized sinusoid has a Fourier series; it is a sum of harmonic sinusoids, that is, sinusoids of harmonic frequencies of a fundamental frequency. Harmonic sinusoids are orthogonal. Therefore, the mean square of the periodic function equals the sum of the mean squares of the harmonic sinusoids. The mean square of the sum of the non-fundamental harmonics can then be calculated by subtracting the mean square of the fundamental from the mean square of the periodic function. For a quantized cosine wave $\mathrm{round}\big(A\cos(x)\big)$ this allows to calculate total harmonic distortion (THD) as:

$$\mathrm{THD} = \sqrt{\frac{\displaystyle\mathrm{MS} - a_1^2/2}{a_1^2/2}} = \sqrt{\frac{\mathrm{MS}}{a_1^2/2} - 1},\tag{4}$$

where MS is the mean square of the quantized cosine wave, $a_1^2/2$ is the mean square of the fundamental, and $a_1$ is the coefficient of the fundamental frequency cosine in the Fourier series of $\mathrm{round}\big(A\cos(x)\big)$, calculated using Eq. 3 of my answer to a related question and simplifying to:

$$a_1 = \frac{2}{\pi A}\left(\mathrm{round}(A)\sqrt{4A^2 - \big(2\,\mathrm{round}(A) - 1\big)^2} + \sum_{k=1}^{\mathrm{round}(A)-1}k\left(\sqrt{4A^2 - (2k - 1)^2} - \sqrt{4A^2 - (2k + 1)^2}\right)\right).\tag{5}$$

Mean square of the quantized cosine wave is calculated in its first quarter-period by:

$$\mathrm{MS} = \frac{1}{\pi/2}\int_0^{\pi/2}\mathrm{round}\big(A\cos(x)\big)^2dx = \frac{\mathrm{round}(A)^2}{\pi/2}\,\mathrm{acos}\left(\frac{\mathrm{round}(A) - 0.5}{A}\right) + {\frac{1}{\pi/2}\sum_{k=1}^{\mathrm{round}(A) - 1}}k^2\Bigg(\mathrm{acos}\left(\frac{k - 0.5}{A}\right) - \mathrm{acos}\left(\frac{k + 0.5}{A}\right)\Bigg)\tag{6}$$

THD is calculated and minimized by the following Python script (continued), using Eqs. 4, 5, and 6:

def a_1(A):
    A = mp.mpf(A)
    return 2*(mp.floor(A + 0.5)*mp.sqrt(4*A**2 - (2*mp.floor(A + 0.5) - 1)**2) + mp.nsum(lambda k: k*(mp.sqrt(4*A**2 - (2*k - 1)**2)-mp.sqrt(4*A**2 - (2*k + 1)**2)), [1, mp.floor(A + 0.5) - 1]))/(mp.pi*A)

def MS(A):
    A = mp.mpf(A)
    return mp.floor(A + 0.5)**2*mp.acos((mp.floor(A + 0.5)-0.5)/A)/(mp.pi/2) + mp.nsum(lambda k: k**2*(mp.acos((k - 0.5)/A) - mp.acos((k + 0.5)/A)), [1, mp.floor(A + 0.5) - 1])/(mp.pi/2)

def STHD(A):  # Square of THD
    MS_1 = a_1(A)**2/2
    return MS(A)/MS_1 - 1

for m in range(2, 25):
    optimal_A = to_max_digits(lambda: mp.findroot(lambda A: mp.diff(STHD, A), approx_optimal_A(m)), 80, 100, double_digits)
    B = to_max_digits(lambda: a_1(mp.mpf(optimal_A)), 80, 100, double_digits)
    THD_optimal = to_max_digits(lambda: 10*mp.log10(STHD(mp.mpf(optimal_A))), 80, 100, dB_digits)
    THD_1 = to_max_digits(lambda: 10*mp.log10(STHD(2**(m-1)-1)), 80, 100, dB_digits)
    THD_2 = to_max_digits(lambda: 10*mp.log10(STHD(2**(m-1)-0.5)), 80, 100, dB_digits)
    print(str(m)+"&"+optimal_A+"&"+B+"&"+THD_optimal+"&"+THD_1+"&"+THD_2+"\\\\")

Figure 2. THD (dB) as function of unquantized amplitude $A$. The local maxima are at amplitudes that cause a narrow protrusion to poke out at the extrema of the quantized sinusoid.

Figure 3. Amplitude $a_1$ (blue solid line) of the fundamental frequency in the quantization of a cosine wave of amplitude $A$, with the identity line (orange dashed line) plotted for reference.

Table 2. Optimal $A$ by definition 2 for different $m\le24$ and the resulting THD, with the THD for some common choices of $A$ listed for comparison. Surprisingly, THD is minimized by the same $A$ that maximize SNR (Table 1), also when tested with much higher precision than what is shown here. If the signal is considered to be a sinusoid of amplitude $a_1$ instead of amplitude $A$, then SNR values are obtained by flipping the sign of the THD values.

$$\begin{array}{l|l|l|lll}m&\text{optimal }A&a_1&\rlap{\text{THD (dB)}}&&\\
&&A=\text{optimal}&A=\text{optimal}&2^{m-1}-1&2^{m-1}-0.5\\
\hline
2&1.26827949461530&1.17011951948679&-10.7629578&-10.1492078&-10.5728562\\
3&3.23800942121037&3.19552705145171&-18.7633376&-18.2533980&-18.3370094\\
4&7.21658597929407&7.19632525078776&-25.5045573&-25.1895549&-25.0267366\\
5&15.2007181331955&15.1907044658090&-31.8098475&-31.6159563&-31.3681507\\
6&31.1888756714257&31.1838597476996&-37.9357894&-37.8138122&-37.5642711\\
7&63.1800835394190&63.1775601103885&-43.9857175&-43.9071306&-43.6902709\\
8&127.173613625523&127.172343338577&-50.0049084&-49.9531877&-49.7783344\\
9&255.168894736361&255.168255766598&-56.0134634&-55.9788184&-55.8439127\\
10&511.165479188890&511.165158147298&-62.0199994&-61.9964656&-61.8952456\\
11&1023.16302205377&1023.16286093052&-68.0278654&-68.0117064&-67.9371468\\
12&2047.16126264484&2047.16118185697&-74.0380588&-74.0268736&-73.9726322\\
13&4095.16000722516&4095.15996674788&-80.0505958&-80.0428071&-80.0036830\\
14&8191.15911371601&8191.15909344706&-86.0651409&-86.0596937&-86.0316405\\
15&16383.1584789666&16383.1584688212&-92.0812822&-92.0774606&-92.0574286\\
16&32767.1580286428&32767.1580235662&-98.0986407&-98.0959534&-98.0816905\\
17&65535.1577094659&65535.1577069262&-104.116904&-104.115011&-104.104877\\
18&131071.157483397&131071.157482127&-110.135830&-110.134495&-110.127304\\
19&262143.157323352&262143.157322717&-116.155234&-116.154293&-116.149195\\
20&524287.157210089&524287.157209771&-122.174984&-122.174319&-122.170708\\
21&1048575.15712995&1048575.15712979&-128.194979&-128.194509&-128.191953\\
22&2097151.15707326&2097151.15707318&-134.215150&-134.21482&-134.213010\\
23&4194303.15703317&4194303.15703313&-140.235446&-140.235212&-140.233932\\
24&8388607.15700481&8388607.15700479&-146.25583&-146.25567&-146.25476\end{array}$$

The equivalence of the two optimality definitions holds even when numerical precision is increased significantly, here with $m=4$ at least to 200 decimal places, in Python (continued):

m = 4
mp.mp.dps = 200
mp.findroot(lambda A: mp.diff(RelMSE, A), 2**(m-1)-1+0.157)
mp.findroot(lambda A: mp.diff(STHD, A), 2**(m-1)-1+0.157)

which outputs numerically identical optimal values for $A$ for the two optimality definitions:

mpf('7.21658597929406951556806247230383254685067097032105786583650636819627678717747461433940963299310318715204551609940031954265317274195597248077934451075855527')
mpf('7.21658597929406951556806247230383254685067097032105786583650636819627678717747461433940963299310318715204551609940031954265317274195597248077934451075855527')

Limit $m\to\infty$

At large $m$ it becomes difficult to directly optimize $m$ numerically, so another approach is desirable. A Taylor approximation of the sinusoid about its peak, where it most differs from a linear function, is a quadratic polynomial. This can be used to analyze the effects of quantization in the limit $m\to\infty$ $\Rightarrow$ $A\to\infty$. The difference between the MS quantization error of a sinusoid with an amplitude $A$ and the MS quantization error $1/12$ of a linear function is proportional to (Fig. 4):

$$\begin{gather}\mathrm{MS} - \frac{1}{12} \propto f(a) = \int_0^{\sqrt{4a + 2}/2}\left((x^2 - a)^2 - \frac{1}{12}\right)dx + {\sum_{k=1}^{\infty}\int_{\sqrt{4a + 4k - 2}/2}^{\sqrt{4a + 4k + 2}/2}\left((x^2 - a - k)^2 - \frac{1}{12}\right)dx}
= \frac{1}{60}\left(\frac{\sqrt{4a + 2}}{16a^2 - 4a - 1} + \sum_{k=1}^\infty\left(\sqrt{4a + 4k + 2}\left(16a^2 + 4a\,(8k - 1) + 16k^2 - 4k - 1\right)\\
- \sqrt{4a + 4k - 2}\left(16a^2 + 4a\,(8k + 1) + 16k^2 + 4k - 1\right)\right)\right),\end{gather}\tag{7}$$

when the amplitude $A$ is an integer $\to\infty$ plus a real number $-0.5 < a \le 0.5$. The sum in Eq. 7 appears to converge, whereas leaving out the term $-\frac{1}{12}$ would result in the sequence of partial sums growing without bound, indicating that no matter which value of $a$ is chosen, $\left(\mathrm{MS} - \frac{1}{12}\right)/\mathrm{MS} \to 0$ and $\mathrm{MS} \to \frac{1}{12}$ as $m\to\infty$.

Figure 4. $f(a)$ and $f'(a)$ of Eq. 7. The shape of $f(a)$ looks identical to the shape of RRMSE for large $m$ in Fig. 1.

By symbolically differentiating $f(a)$ defined in Eq. 7 with respect to $a$ (Fig. 4) and by finding the zero of $f'(a)$ near $a = 0.157$, the optimal $a$ at $m\to\infty$ can be computed to quite a bit of precision, in Python (continued):

def f(a):
    a = mp.mpf(a)
    return (mp.sqrt(4*a + 2)*(16*a**2 - 4*a - 1) + mp.nsum(lambda k: (mp.sqrt(4*a + 4*k + 2)*(16*a**2 + 4*a*(8*k - 1) + 16*k**2 - 4*k - 1) - mp.sqrt(4*a + 4*k - 2)*(16*a**2 + 4*a*(8*k + 1) + 16*k**2 + 4*k - 1)), [1, mp.inf]))/60

def Df(a):  # Derivative of f(a)
    a = mp.mpf(a)
    return mp.sqrt(2)*(16*a**2 + 4*a - 1)/(12*mp.sqrt(2*a + 1)) + mp.nsum(lambda k: (mp.sqrt(4*a + 4*k - 2)*(16*a**2 + 4*a*(8*k + 1) + 16*k**2 + 4*k - 1) - mp.sqrt(4*a + 4*k + 2)*(16*a**2 + 4*a*(8*k - 1) + 16*k**2 - 4*k - 1))/(12*mp.sqrt(2*a + 2*k + 1)*mp.sqrt(2*a + 2*k - 1)), [1, mp.inf])

to_max_digits(lambda: mp.findroot(lambda a: Df(a), 0.157), 63, 83, double_digits)

which gives $a \approx 0.156936321399$. This can be used to create an approximation of the optimal $A$ as function of $m$, intended for large $m \ge 20$:

$A \approx 2^{m-1} - 1 + 0.156936321399 + e^{-1.2749819017
 - 0.3464088917 m},\tag{8}$

where the coefficients in the exponent were calculated by a linear fit at $m \in \{21, 22\}$ to a linearization of optimal $A$ values of Table 1 or equivalently Table 2. The error from using the approximation is shown in Fig. 5.

Figure 5: Absolute error in the approximation of optimal $A$ by Eq. 8 as function of $m$. For $m \in \{21, 22\}$, which were used for fitting, and for $m \in \{23, 24\}$, the absolute approximation error was less than $10^{-8}$.

Out of curiosity, I also computed the worst-case $a$ that gives the largest MS quantization error as $m\to\infty$, by finding the zero of $f'(a)$ near $a = -0.43$, which turned out to be at $a \approx -0.433510875868$.

Conclusion

The two definitions of optimality appear to be equivalent, to convincing numerical precision. As $m\to\infty$, the optimal value of $A$ approaches approximately $2^{m-1}-1+0.156936321399$ (or more accurately for large finite $m$ the approximation of Eq. 8) and the quantization error reduction (by def. 1 SNR or def. 2 THD) in dB from choosing the optimal value approaches zero, compared to choosing another large value such as $A = 2^{m-1}-1$ or the nearly worst-case choice $A = 2^{m-1}-0.5$.

At the optimal amplitude $A$, the least squares (LS) sinusoid is of the same frequency and phase as the sinusoid being quantized, but has a somewhat lower amplitude $a_1$ given in Table 2. This is somewhat counterintuitive. To minimize THD (or to maximize SNR with the sinusoid being approximated as the signal) of approximating a sinusoid of amplitude $a_1$ using a waveform of $m$-bit numbers with values in range $-2^{m-1} \ldots 2^{m-1}-1$ that is constructed by quantizing (by rounding to the nearest integer) a sinusoid of amplitude $A$, one must choose optimal $a_1$ and $A$ that are not equal.

The results are applicable to continuous time without sampling, or to sampling in the limit $f\to f_s$, where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency. In general, optimality is not preserved by sampling. Random sampling with random times of the samples will preserve definition 1 SNR optimality, if the distribution of the sinusoid phase (modulo $2\pi$) at the samples is uniform. Also, for irrational $f/f_s$, the same optimality is preserved for the same reason, see equidistribution theorem. Definition 1 SNR optimality vanishes with rational $f/f_s$, because the distribution of sinusoid phases will not be uniform.


CP based Symbol Boundary detection algorithm which is fundamentally based on correlation will still work. You just need to do a timing recovery before FFT operation.

A small offset in sampling times will cause a small phase rotation in pilots and these small rotations are used to compute the timing offset.

Once you have a good timing offset estimate, you do timing recovery on the received samples of OFDM symbol. Timing recovery is nothing but doing some low cost polynomial interpolation to get the actual values of the samples at correct sampling times.

Or in OFDM once you have these timing offsets as phase rotation, derotate rest of the data symbols after FFT before slicing to get the complex QAM symbols.


Incremental Online Update of Tf-IDF Indices for Image Retrieval

Ok, apparently, it is not much different from an update of the mean, given the accumulated statistics (mode statistics). Basically, when the mean and the number of elements are known, the idf weights can be re-computed whenever a new element is inserted. This can be achieved in two ways:

Allocate your tree with the maximum number of features. Keep all the elements inserted. Whenever a new insert occurs, re-compute TF-IDF over the entire set.
Only keep the statistics which are required to update the weights. This way, an entire pass over the data is not necessary.

Normalization is again performed over the entire dataset once the idf is updated. Removal is not much different from the insertion, but one needs to use a suitable pointer-wise data structure if it is desired to have the computational complexity both for insertion and deletion.


How to detect the maximum resolvable spatial frequency of camera?

The image that you have obtained is already capturing the Modulation Transfer Function (MTF) of the optical system under test. It does not need to be further passed via an extra Discrete Fourier Transform (DFT) step. The spatial resolution limit is a matter of overall specification.

The test pattern's spatial resolution is measured in lines per millimeter and it varies in (a typical) $x$ direction (the width of the image) increasing towards the right.

As the resolution is increasing, the contrast (between the dark and light areas of the pattern) is decreasing.

Therefore, the simplest way to obtain the MTF is to create a scatter plot whose $x$ axis is spatial resolution in 1, 2,3,... lines/mm and the $y$ axis is the contrast (or modulation factor for that particular frequency) as $\frac{I_{max}-I_{min}}{I_{max}+I_{min}}$.

BUT, as you can see, the difference between $I_{max}, I_{min}$ in the test image is dropping to zero (blue trace in your diagram). Therefore, the test pattern (image) itself does not represent a $1$ input to the optical system across all spatial frequencies.

Therefore, it would be best to calculate the MTF of the test pattern and then the MTF of the system under test and present both of them on the same diagram if you want to depict how close to the "ideal" the optical system under test really is. This is the point of the MTF.

To find the smallest object that this optical system can resolve, you need to translate the size of that object into lines/mm and check that the MTF at that limit is above the noise floor of the sensing element.

Hope this helps.


Recovering signal from magnitude and phase of FFT

I tried running the same code, and it works fine, except for 2 minor points:

It's preferable to use exp(1i*phase) instead of exp(sqrt(-1)*phase)
The final result will deviate from the original signal due to rounding noise (typically around 1e-16). Also, some of this noise will appear in the imaginary part of the result (so, even if y was real, yinv will be complex, due to the small imaginary component). This will be particularly apparent if you try: plot(yinv), because Matlab will see that yinv is complex, and plot it on the real/imag axis


SNR contains $\log N$ term where $N$ stands for number of samples

In the additive model $y=s+n$, when the signal is deterministic, it adds coherently over the "realizations". Hence, its variance $V(\sum s_n) = V(N s) = N^2  V( s)$. And when the noise $w$ is independent identically distributed (IID), then $V(\sum w_n) = NV( n) $. This is a classical result on the  Variance of Uncorrelated Variables.

You can find this classical result detailed in Deriving the SNR for averaged signals. Warning, this does reach a limit when data is quantized, in other words the noise floor does not follow $1\sqrt{N}$, but plateaues.


From Fourier transform to Laplace Transform

So it sounds like X(s) is the analytic continuation of X(f), given the proper restrictions, is that right?

that's what i would say.  be careful about expressing the function and its argument.
use a different "$X(\cdot)$" if you need to.  like:
$$ X_\mathscr{L}(j 2 \pi f) = X_\mathscr{F}(f) $$
where
$$ X_\mathscr{L}(s) \triangleq \mathscr{L}\Big\{ x(t) \Big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-st} \, \mathrm{d}t $$
and
$$ X_\mathscr{F}(f) \triangleq \mathscr{F}\Big\{ x(t) \Big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2\pi f t} \, \mathrm{d}t $$


It just means "the transformation that turns $x$ into $y$."

You might also see $\mathbf{T}^{-1}$ which means the inverse: turning $y$ into $x$.


HPS Algorithm for detecting the fundamental frequency of a guitar string

OK Lets go!

I consider HPS a very simple algorithm, this image representation show exactly how HPS work !

Yeah is highly recommended apply a hann window for every framed data !

What do you need to do is apply one window function over your framed signal, then apply FFT, you need just the first half absolute values from FFT, its give you the Magnitude, now downsample it for N time, Multiply all downsampled and at  the end find the peak.

  How can I know how many times I need to downsample the signal?

Will depend, the more time you downsample you signal more accurate HPS can be, this has a cost that limits the frequency analysis, for example, imagine that you have a signal sampled at 44100hz, you use FFTsize = 4096 and you decide to make 5 time downsample, its tell me that you can not pick up frequencies above 4403Hz!

floor(4096 / 5) =  819
floor(819 /2) = 409
floor((44100/4096) * 409) = 4403


Why does knowing the impulse response allow you to determine the output for any LTI system?

(Adapted from this answer on dsp.SE)

The reason that the impulse response (also called the unit pulse response for discrete-time systems) determines the output for arbitrary input $x$ to an LTI system is that

  The output of a linear time-invariant system in response to input $x$ is
  the sum of scaled and time-delayed versions of the
  impulse response.

This is because of the linear property of the system: if the response to input signal $x_i$ is $y_i$, then the response to input $x_1+x_2+ x_3 + \cdots$ is $y_1+y_2 + y_3 + \cdots$.  So, let's decompose the input signal $x$ into a sum of scaled
unit pulse signals. The system response to the unit pulse signal
$\cdots, ~0, ~0, ~1, ~0, ~0, \cdots$ is
the impulse response or pulse response
$$h[0], ~h[1], \cdots, ~h[n], \cdots$$
and so by the scaling property the single input value $x[0]$,
or, if you prefer
$$x[0](\cdots, ~0, ~0, ~1, ~0,~ 0, \cdots)
=  \cdots ~0, ~0, ~x[0], ~0, ~0, \cdots$$
creates a response
$$x[0]h[0], ~~x[0]h[1], \cdots, ~~x[0]h[n], \cdots$$

Similarly, the single input value $x[1]$ or creates
$$x[1](\cdots, ~0, ~0, ~0, ~1,~ 0, \cdots)
=  \cdots ~0, ~0, ~0, ~x[1], ~0, \cdots$$
creates a response
$$0, x[1]h[0], ~~x[1]h[1], \cdots, ~~x[1]h[n-1], x[1]h[n] \cdots$$
Notice the delay in the response to $x[1]$.  We can continue further
in this vein, but it is best to switch to a more tabular form
and show the various outputs aligned properly in time. We have
$$\begin{array}{l|l|l|l|l|l|l|l}
\text{time} \to & 0 &1 &2 & \cdots & n & n+1 & \cdots \\
\hline
x[0] & x[0]h[0] &x[0]h[1] &x[0]h[2] & \cdots &x[0]h[n] & x[0]h[n+1] & \cdots\\
\hline
x[1] & 0 & x[1]h[0] &x[1]h[1] & \cdots &x[1]h[n-1] & x[1]h[n] & \cdots\\
\hline
x[2] & 0 & 0 &x[2]h[0] & \cdots &x[2]h[n-2] & x[2]h[n-1] & \cdots\\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots & \\
\hline
x[m] & 0 &0 & 0 & \cdots & x[m]h[n-m] & x[m]h[n-m+1] & \cdots \\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}$$
The rows in the above array are precisely the scaled and
delayed versions of the impulse response that add up to
the response $y$ to input signal $x$.

So, we see that if we know the impulse response of an LTI system, we can figure out the response of the LTI system to any arbitrary signal $x$. We don't need any additional information about the LTI system beyond the impulse response.


Harris corner detection shape of $E(u,v)$

$E(u,v)$ is zero for zero shift ($u=v=0$), where the two patches being compared are equal. As you increase the distance, the difference tends to increase. Depending on the local shape of the image, the $E(u,v)$ surface will have a different shape: at a corner or a small point, it will steeply increase from 0 in all directions. Along an edge, it will remain close to 0 along one direction, and increase steeply in the perpendicular direction. At flat image regions, the function will remain closer to zero in all directions.
The Harris corner detector (and other corner detectors too) look at the shape of this patch, if it increases steeply in all directions, the detector will produce a stronger response.


direct Frequency domain FIR filtering vs Overlap-add method

An FFT of twice the filter order is kind of short.  Take a look at the "Choice of FFT size" section from this article I wrote a while back.
Also your whole signal will fit into a reasonable size FFT, so I think the single FFT approach is fine in this case.


Calculating signal/noise ratio and the noise level

The decibels (dB) provides a logarithmic scale for the ratio of two power values. Mathematically it is defined as:
$$
  R = 10 \log_{10} \left(\frac{P}{P_0}\right)
$$
where $P$ and $P_0$ are power values in linear scale (such as Watts) and $R$ is the ratio expressed in logarithmic scale. Often a fixed reference power level (such as 1 Watt) is chosen for $P_0$. When that reference power level $P_0$ is 1 milliwatt (mW), the resulting unit for the $R$ quantity is called a dBm.

Given a known reference level, you can conversely obtain the power in linear scale from the power ratio value expressed in logarithmic scale by using the inverse formula:
$$
  P = 10^{R/10}\cdot  P_0
$$

So for a measurement of -46dBm (remember that the dBm units means that the reference power level is 1mW), the resulting power is:
$$
\begin{align}
  P_s &= 10^{-46/10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-5} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-8} \, \mbox{W} \\
      &\approx 2.5 \times 10^{-2} \, \mbox{µW}
\end{align}
$$
Similarly for a measurement of -96dBm, the resulting power is:
$$
\begin{align}
  P_n &= 10^{-96/10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-10} \cdot 1 \mbox{mW} \\
      &\approx 2.5 \times 10^{-13} \, \mbox{W} \\
      &\approx 2.5 \times 10^{-4}  \, \mbox{nW}
\end{align}
$$

Now for the signal-to-noise ratio (SNR), it is simply defined as the ratio of the signal power to the noise power. So in linear scale, that would be:
$$
\begin{align}
  \mbox{SNR} &= \frac{P_s}{P_n} \\
             &\approx \frac{2.5 \times 10^{-8}}{2.5 \times 10^{-13}} \\
             &\approx 10^{5}
\end{align}
$$
which when expressed in decibels would give us:
$$
\begin{align}
  \mbox{SNR}_{\mbox{dB}} &= 10\log_{10} \left(\frac{P_s}{P_n}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(10^{5}\right) \, \mbox{dB} \\
             &= 50\mbox{dB}
\end{align}
$$
We could of course also have taken a bit of a shortcut to compute that SNR directly from the power measurements in dBm by using logarithmic identities to get:
$$
\begin{align}
  \mbox{SNR}_{\mbox{dB}} &= 10\log_{10} \left(\frac{P_s}{P_n}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(\frac{\frac{P_s}{P_0}}{\frac{P_n}{P_0}}\right) \, \mbox{dB} \\
             &= 10\log_{10} \left(\frac{P_s}{P_0}\right) - 10\log_{10}\left(\frac{P_n}{P_0}\right) \, \mbox{dB} \\
             &= \left(-46\,\mbox{dBm}\right) - \left(-96\,\mbox{dBm}\right) \\
             &= 50\mbox{dB}
\end{align}
$$


Filter to eliminate a 60 Hz frequency range

The OP is trying to implement a notch filter with a 3 tap FIR filter (as restricted by the problem he is trying to solve). It must be noted that this would result in a very poor filter implementation: It will reject 60 Hz perfectly but will have unavoidable attenuation over a very broad section of the desired passband range. An actual practical implementation would be done with an IIR approach in order to provide a tight notch, as demonstrated by this link: Transfer function of second order notch filter
As a hint to determine the filter within this FIR restriction, as an exercise only and of no real practical value proceed as follows:
Take the Z transform of the difference equation.
Determine the transfer function form as a polynomial ratio $H(z) = Y(z)/X(z)$
From this know that the frequency response is determined by restricting $z$ to the unit circle ($z= e^{j\omega}$) with $\omega$ extending uniquely over the normalized radian frequency from $-\pi$ to $+\pi$ which corresponds to $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate.
Given this is an FIR structure, the transfer function will result in simply a 2nd order polynomial in the numerator (whos roots are the "zeros" of the filter). Place the zeros of this polynomial such that they corresponds to the 60 Hz locations on the unit circle (+/- 60 Hz for a real filter, mapped to the angular frequency as described above), which would provide complete attenuation at f = 60 Hz.


Impulse response of a causal system from transfer function in z-domain

Let's follow option 1. I'm using Python but same is directly available in matlab. The partial fraction expansion can be found via residue function

>>> import scipy.signal as sig
>>> b = [1, 1]
>>> a = [1, 1, 0.5]
>>> r, p, k = sig.residue(b,a)
>>> r
array([0.5-0.5j, 0.5+0.5j])

>>> p
array([-0.5+0.5j, -0.5-0.5j])

>>> k
array([0.])

This gives us
$$
H(z) = \frac{0.5-0.5i}{z+0.5-0.5i} + \frac{0.5+0.5i}{z+0.5+0.5i}
$$

Now either by computing manually, or using a Z-transform table, or using another computational tool you can show that these are two complex terms can be inversed transformed as

$$
\mathcal{Z}^{-1}\left[\frac{a}{z+a}\right] = \mathcal{Z}^{-1}\left[\frac{az^{-1}}{1+az^{-1}}\right] = -(-a)^k \mathbf{1}(k - 1)
$$
Here $\mathbf{1}(k)$ is the discrete time step function. But since we have complex numbers this looks like our response should be a complex function. However this is not possible since our transfer function has real valued coefficients. Hence even before I compute anything, I know already that imaginary parts will cancel out.

And since $|a| < 1$ this is going to be a decaying sinusoidal. Indeed if we check this with say sampling period $T=1$,

>>> import harold as har
>>> G = Transfer([1, 1], [1, 1, 0.5], dt=1)
>>> impulse_response_plot(G);

As expected, there is a one sample delay coming from $\mathbf{1}(k-1)$


Drop the first sample - apply the formulation as if x[2] was x[1].


Exercise Related to Frequency Resolution and SNR

Your calculation of the SNR is right and the DFT Bin Resolution is also OK.
One thing you're missing is the effective resolution due to "Windowing" effect.
The DFT of finite number of samples is interpolated by a Dirichlet Kernel (Like Sinc).  It means you resolution is also limited by the main lobe width of the Sinc which is proportional to the inverse of the window length.
Namely, using less samples means wider lobe which reduces the resolution and the ability to differ between two harmonic signals.


What is this sinc doing in my LP-turned-HP filter output?

You are doing a backwards-forwards filter for one, and a forwards-only filter for the other.  Thus the two filters will behave differently at one or both edges.

If your signal is strongly stationary and exactly periodic in the data window length, you could try wrapping the data around circularly and start your filters outset by their impulse response lengths.


By any reasonable metric, the dB value for zero amplitude is minus infinity.


Error plot between known and estimated data

Your formula/method for computing MSE between estimated and known inputs looks good to me. For symbol error rate you could use something like a Hamming distance which simply counts the number of times the estimated symbol is different from the actual symbol. In your 10 symbol example the error rate is 4/10 i.e 40%.

In Matlab you can do something like:

symbol_error_rate_percent = sum(u_hat ~= u)/length(u) * 100

You can compute this for $T$ trials and average them, like you did for MSE. Or you could just make a really long sequence of length $10 \times T$ and do it once.


What does this signal filter do? (coming from finance)

I thought you stock market tech analysts always used moving averages (rolling means) or "exponential moving averagers" for your smoothing operations. Anyway, your $macd[n]$ process is neither of those. As far as I can tell the frequency magnitude response of your $macd[n]$ process looks like the following:

Your $macd[n]$ process does not smooth its input data in the normal sense. It amplifies input data that has a periodicity (repetition rate) of 33.3 days per cycle. In the above plot you can clearly see the "mainlobe" that Fat32 was referring to. If your closing price data contains a periodicity of roughly 33 days then the output of your $macd[n]$ process will be a large value. Otherwise, the output of your $macd[n]$ process will be a small value. I don't know if such a frequency response makes stock market technical analysis sense, or not.

As for your $h[n]$ process, I believe its frequency magnitude response looks like the following:

Your $h[n]$ process appears to be a process that amplifies input data that has a periodicity (repetition rate) of 21.9 days per cycle. In DSP we use processes like your $macd[n]$ and $h[n]$ for filtering when we want to detect the presence of components in our input signal that have periodicities of 33.3 or 21.9 time units per repetition. We call such processes "resonators" or "bandpass filters."


BER performance of a space-time code over Rayleigh channel model

The classic approach is to evaluate the union bound of the pairwise error probability, taken over the code matrices.

If $X_0$ and $X$ are two different code matrices, then their pairwise error probability is the probability that the receiver mistakes $X$ for $X_0$, when $X_0$ is transmitted. It is written $P(X_0 \rightarrow X)$.

Assuming all code matrices are transmitted with equal probability, and that the code is $\mathcal{X}$, then the average error rate is upper bounded by
$$ P_{\text{err}} \leq \frac{1}{|\mathcal{X}|} \sum_{X_n \in \mathcal{X}} \sum_{X_k \in \mathcal{X},k \neq x} P(X_n \rightarrow X_k) $$
which is just an application of the union bound.

The calculation of BER from $P_{\text{err}}$ for specific channels can get complicated and in fact in many cases no closed form formula can be found. I would suggest the book "Space-Time Block Coding for Wireless Communications" by Larsson and Stoica for more details, in particular section 4, "Error Probability Analysis".


Matlab: xcorr: Having trouble when I try to calculate delay between two sinusoids

A sine wave is cyclo-stationary so is not the best waveform to use for estimating delay between transmit and received signals. It will correlate whenever the delay causes the signal to be an additional $2\pi$ in phase offset. Consider using waveforms with better autocorrelation properties for this purpose such as pseudo-random noise sequences or other generated white noise sources.
Regardless, the effect of frequency offset on cross-correlation is that the magnitude rolls-off as a Sinc function with frequency offset with the first null at $1/T$ where $T$ is the total time duration of the signal, and the phase decreases linearly from $0$ at $f=0$ to $-\pi$ at $f=1/T$. At a fixed frequency offset between the signals, a fixed phase offset would result, and for a narrow band signal this can be approximated as a fixed delay according to the carrier frequency of the signal.
Consider that correlation is the sum of the element by element product between two signals (the dot product), and observe how the product of two signals will be at the frequency difference between them (the beat note). Thus if your total time duration is exactly 1 complete cycle off in frequency, the result will average to zero (showing how the first null is at $1/T$).
The following code demonstrates this experimentally with a random white noise waveform of total time duration 0.1 seconds:
T = .1;
fs = 10e3;
N = T*fs;
t = [0:N-1]/fs;
x = randn(1,N);

index = 0;
foffset = 0:.01:20;
for f = foffset
    index++;
    y = x.*exp(j*2*pi*f.*t);
    # dot product is the correlation, and normalized
    corr(index) = x*y'/(N*std(x).*std(y));
end

figure
plot(foffset,abs(corr))
figure
plot(corr)

This creates the two plots shown below with the first plot of the magnitude of the normalized correlation versus the static frequency offset between the two signals, and the second plot as a polar plot showing both magnitude and phase components.


Generating a timeseries with an arbitrary power spectrum

You are reversing your phases okay.  You should be able to use an even value of N if you set the phase for N/2 to zero.

The problem you have is that your PSD spans the entire N.  Unless $ P_k = P_{N-k} $ you are not generating a DFT with complex conjugate symmetry which is what you need to produce a real valued signal.

Also, you should be using the "ifft" instead of the "fft" to do an inverse FFT.

Hope this helps,

Ced

Edit:  I didn't see your edit before I posted.  It seems you figured this out.

======================

Followup:

For even N:

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2-1)*2*np.pi
    phase[1:N/2] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase


It's generally not true that the average deviation from the mean is not important or much less important than the power of the deviation. For sampled data, the mean deviation would be defined as

$$\epsilon=\frac{1}{N}\sum_{n=0}^{N-1}|x_n-\mu_x|\tag{1}$$

where $\mu_x$ is the mean of the data. This quantity is definitely relevant. However, the absolute value in $(1)$ makes it much less convenient to use than the equally relevant standard deviation

$$\sigma=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}|x_n-\mu_x|^2}\tag{2}$$

which is differentiable, and, consequently, easier to handle analytically (such as in least squares fitting, etc.). (Note that in Eq. $(2)$ you might as well use a factor $1/(N-1)$ instead of $1/N$, but that is not relevant here).

Apart from the above argument based on mathematical convenience, it is true that in the case of uncorrelated noise sources, it is the noise variances (i.e., the squared standard deviations) that add up, and not the mean deviations.


The $k$-th FFT bin will contain the complex amplitude for frequency $\frac{k}{512} \times 44100$.