Calculating colour temperature from two sources of varying luminous intensity

LEDs do not emit black-body (Planckian) radiation which has a specific spectrum, so instead of color temperature which concerns black-body radiators we have to resort to correlated color temperature (CCT). It is the temperature of a black body that emits light of a color that is perceptually closest to the color of light of the same brightness from, in this case, the LEDs. The CIE 1960 $(u,v)$ chromaticity space ignores brightness and is designed to be perceptually uniform. It is used to calculate correlated color temperature and is based on the CIE 1931 $(X, Y, Z)$ color space.

Your LED datasheet gives for the two LED types their nominal fluxes in lumens and also gives their CIE 1931 $(x, y)$ chromaticities. Flux in lumens is calculated from a light spectrum using the CIE 1924 photopic (normal color vision, not low-light) luminosity function which is identical to the CIE 1931 spectral efficiency for the $Y$ color component, by definition. So we can work with normalized $(X, Y, Z)$ color where $Y$ equals the flux in lumens. This normalization is compatible between LED types of different chromaticities. Knowing $Y$ and the chromaticity $(x, y)$, we can easily calculate the rest of the $(X, Y, Z)$ triplet describing a LED type's light color, using the CIE 1931 definitions of $x$ and $y$:

  $$x = \frac{X}{X+Y+Z},\quad y = \frac{Y}{X+Y+Z}\tag{1}$$

$$\Rightarrow\quad X = \frac{Y\,x}{y},\quad Z = \frac{(1 - x - y)\,Y}{y}.\tag{2}$$

The $(X, Y, Z)$ color space is a linear-intensity color space so the $(X, Y, Z)$ colors of the LEDs sum additively in additive mixing of their lights. The correlated color temperature is then calculated from the summed $(X, Y, Z)$ color by going back to $(x, y)$ using the above definitions and by converting to CIE 1960 $(u, v)$ by the definition of that chromaticity space in relation to CIE 1931 chromaticity:

  $$u = \frac{4x}{12y - 2x + 3},\quad v = \frac{6y}{12y - 2x + 3}.\tag{3}$$

These steps could also be combined to a single equation. What remains to do is to find the temperature of a black-body that gives $(u, v)$ coordinates (Planckian locus) nearest to this point by straight-line distance. Krystek, M. (1985), An algorithm to calculate correlated colour temperature. Color Res. Appl., 10: 38-40 recommends using a bisection search together with the following rational Chebyshev approximation valid for $1000 < T < 15000$ where $T$ is temperature in kelvins:

  $$\begin{eqnarray}u = \frac{0.860117757+1.54118254 \times 10^{-4}T + 1.28641212 \times 10^{-7} T^2}{1+8.42420235 \times 10^{-4}T + 7.08145163 \times 10^{-7}T^2},\\
v=\frac{0.317398726+4.22806245 \times 10^{-5}T + 4.20481691 \times 10^{-8} T^2}{1-2.89741816 \times 10^{-5}T+1.61456053 \times 10^{-7}T^2}.\end{eqnarray}\tag{4}$$

The Colour package (install using pip install colour-science) for Python includes the Krystek 1985 method but I just found a bug in their implementation in version 0.3.14 which will be fixed in version 0.3.15. They also have other similar methods in function uv_to_CCT and in other related functions. Example:

import numpy as np
import colour

# Subscript 0 denotes 1800K LED, 1 denotes 4000K LED.
# LED chromaticity center point:
xy_0 = np.array([0.5656, 0.4136])
xy_1 = np.array([0.3860, 0.3870])
# LED nominal flux (lumens) at nominal forward current:
flux_0 = 825
flux_1 = 1370
# Mixing weight (lumens)
weight_0 = flux_0
weight_1 = flux_1

# Proportional XYZ color at flux of 1 lumen (setting Y = 1):
XYZ_0 = np.array([xy_0[0]/xy_0[1], 1, (1 - xy_0[0] - xy_0[1])/xy_0[1]])
XYZ_1 = np.array([xy_1[0]/xy_1[1], 1, (1 - xy_1[0] - xy_1[1])/xy_1[1]])

# Correlated color temperatures of the two LED types
cct_0 = colour.xy_to_CCT(colour.XYZ_to_xy(XYZ_0), method='McCamy 1992')
cct_1 = colour.xy_to_CCT(colour.XYZ_to_xy(XYZ_1), method='McCamy 1992')

# Correlated color temperature of a mixture of the two LED types each at nominal current
colour.xy_to_CCT(colour.XYZ_to_xy(weight_0*XYZ_0 + weight_1*XYZ_1), method='McCamy 1992')

The obtained correlated color temperatures of the two LED types are 1746 and 3934 kelvins, close to their nominal CCTs of 1800 K and 4000 K. At nominal forward current for each LED type, the CCT of the mixture is 2715 kelvins. So the script seems to be working as intended.

We can also get a color temperature curve as function of mixing ratio, and see that it is not exactly linear, although this depends on the definition of a mixing ratio. Continuing in Python:

import matplotlib.pyplot as plt

ccts = np.zeros(101)
for percent_1 in range(101):
    weight_0 = (100-percent_1)*flux_0/100
    weight_1 = percent_1*flux_1/100
    ccts[percent_1] = colour.xy_to_CCT(colour.XYZ_to_xy(weight_0*XYZ_0 + weight_1*XYZ_1), method='McCamy 1992')

plt.plot(range(101), ccts, "-", range(101), cct_0*(100 - np.arange(101))/100 + cct_1*np.arange(101)/100, "--")
plt.ylabel("correlated color temperature (K)")
plt.xlabel("4000 K LED flux / nominal flux (%)")
plt.grid()
plt.show()

Figure 1. Correlated color temperature (solid blue) of a mixture of 1800 K and 4000 K LED types at different mixing ratios. 0 % and 100 % mean a single LED type operating at its nominal flux at its nominal forward current and the other LED type off. 50 % means each LED type is operating at half of its nominal flux. Calculations use manufacturer's LED data. A linear relation (dashed orange) is shown for comparison.


Bandwidth range for Fast Fourier vs principal component analysis?

From Karhunen–Loève theorem, when talking about stochastic processes:

  In the theory of stochastic processes, the Karhunen–Loève theorem
  (named after Kari Karhunen and Michel Loève), also known as the
  Kosambi–Karhunen–Loève theorem is a representation of a
  stochastic process as an infinite linear combination of orthogonal
  functions, analogous to a Fourier series representation of a function
  on a bounded interval.

Basically, K.-L. yields an adaptive representation, while Fourier provides a fixed representation with  sinusoidal functions.

There are many flavors of Fourier tools, so since you are talking about "fast" Fourier and PCA, I assume for now you are dealing with discrete data, thus the empirical version of the  Karhunen–Loève transform (PCA) and the discrete Fourier transform (DFT). The idea of a continuous set of frequencies does not fit gracefully here.

Both PCA et discrete Fourier can be cast into linear and orthogonal transforms, so they are applicable to any data, and may preserve all information.

From Principal Component Analysis (p. 44 sq.), you can find that you can derive DFT bases from PCA, when you are studying a process that follows a "correlated" Markov model.

Finally, when you talk about infinite sines with infinitesimal bandwidth, you are more in the context of continuous functions (or distributions) and standard Fourier analysis, which might be a cause for confusion. A solid book on the topic could help, like Mathematical principles of signal processing: Fourier and wavelet analysis, by P. Brémaud.


Estimate changing time lag between signals

If the signals are as you've drawn them (flat, with abrupt changes), then perhaps this approach might work:

Take the absolute value of the differences between successive time samples.
The non-zero entries are the transition times.
Find these times for both signals.
Subtract the two.

I've made an attempt to do this in the R code below.

Example output:

 > print(differences_in_delay)
  [1]  2 -5  2 -4 -1  5 -4  3 -3  0
 > print(t[which(abs(diff(x2))>0)] - t[which(abs(diff(x1))>0)])
  [1]  2 -5  2 -4 -1  5 -4  3 -3  0

It breaks down when the transitions are within a sample or two of each other, but otherwise seems to work OK.

Will also probably break down if there's noise on your signals.

R Code Below

T <- 1024
t <- seq(1,T,1)
Nchanges <- 10
times_for_changes_for_x1 <- sort(runif(Nchanges,T/10,9*T/10))
differences_in_delay <- round(runif(Nchanges,-5,5))
times_for_changes_for_x2 <- times_for_changes_for_x1 + differences_in_delay
changes_for_x1 <- runif(Nchanges,-5,5)

x1 <- 0*t
x2 <- 0*t
change_no_x1 <- 1
change_no_x2 <- 1

for (tm in t)
{
  if ( (change_no_x1 <= Nchanges) && (tm > times_for_changes_for_x1[change_no_x1]))
  {
    x1[tm] <- x1[tm] + changes_for_x1[change_no_x1]
    change_no_x1 <- change_no_x1 + 1
  }
  else
  {
    if (tm > 1)
    {
      x1[tm] = x1[tm-1]
    }
  }
  if ( (change_no_x2 <= Nchanges) && (tm > times_for_changes_for_x2[change_no_x2]))
  {
    x2[tm] <- x2[tm] + changes_for_x1[change_no_x2]
    change_no_x2 <- change_no_x2 + 1
  }
  else
  {
    if (tm > 1)
    {
      x2[tm] = x2[tm-1]
    }
  }
}

plot(x1,type='l')
lines(x2,col="red")

plot(abs(diff(x1)>0)*1, type='l')
lines(abs(diff(x2)>0)*1, col="red")

print(differences_in_delay)
print(t[which(abs(diff(x2))>0)] - t[which(abs(diff(x1))>0)])


As MattL already mentioned, they are the same block diagrams as Direct Form II (canonical) implementation of a third order (on the FIR side) LCCDE.

However, the adders in the first block diagram are two input adders, while the the ones in the second block diagram are multiple input adders. Theoretically there is no difference.

In a popular DSP text, only the former block diagram is displayed, with sum (sigma) symbol replaced with a + sign.


Understanding the meaning of amplitude in FFT

The continuous-time Fourier transform of a function $f(t)$ is in essence an integration of $f(t)$ multiplied with a complex exponential kernel:
$$ F(\omega)=\int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt \tag{1}$$
Since the exponential function is unitless, the unit of the Fourier integral will be the multiplication of the units of the function $f(t)$ and the differential $dt$.
Assuming that the function $f(t)$ had a unit of micro Tesla, and its argument $t$ is time (in seconds), then the unit of $dt$ will be seconds. As a consequence, the unit of the Fourier transform, $F(\omega)$, will be  micro Tesla  second  $$\mu T \cdot s \tag{2}$$
However, what you actually compute is the discrete-time Fourier transform, $F(e^{j\omega})$, of the samples $f[n]= f(nT_s)$ of the function $f(t)$, via the summation: $$F(e^{j\omega}) = \sum_{n=-\infty}^{\infty} f[n] e^{-j\omega n} \tag{3}$$ where $T_s$ is the sampling period in seconds.
Furthermore, instead of the continuous-argument function $F(e^{j\omega})$, you will compute its samples $F[k]$
$$ F[k] = \sum_{n=0}^{N-1} f[n] e^{-j \frac{2\pi}{N} n k} \tag{4}$$
through a DFT (discrete Fourier transform) of the samples $f[n]$ of length $N$, possibly implemented with an FFT algorithm.
The unit of the samples $f[n]$ is the same as the unit of $f(t)$, making the unit of $F(e^{j\omega})$  and $F[k]$ as micro Tesla. Therefore the unit of the FFT samples $F[k]$ of $F(e^{j\omega})$ will be micro Tesla.
Note that there is an (implicit) amplitude scaling by $1/T_s$ in the computed DFT samples $F[k]$, and when you want to display the continuous-time Fourier transform $F(\omega)$ from the samples $F[k]$, you multiply them with $T_s$, which corrects not only the amplitude scaling, but also the unit of it, by making it micro Tesla second as in (2).


Apply Transfer Function in Continuous Domain in Matlab

One approach could be to get a finite number of values of the impulse response of your system using the 'impulse' function, and then convolve this with your signal (taking care of the actual sampling frequency you are using).

% sampling time of 1ms, 20s of time signal
T = 0.001;
t = 0:T:20;

% low pass filter with cutoff frequency at sine frequency
wo = 1;
x = sin(wo*t);
s = tf('s');
H = wo / (s+wo);

% get the h(t)
h = impulse(H,t);

% plot the output
plot(t,x);
y = T*conv(x,h);
y = y(1:length(x));
hold all;
plot(t,y)


After seeing the code, I think the mistake may be

sinc((t(i)-n*Ts1)/Ts1)

Change it to

sinc((t(i)-(n-1)*Ts1)/Ts1)

An extra Ts1 delay was being applied to sinc. The index $n$ starts from 1 for MATLAB but it should start from 0.


What is the need for prediction filter in PCM and DPCM?

The DPCM works by the difference between actual samples and predicted samples.

  If we already have the actual samples, we can quantize it and encode it later.

  But why we use prediction filter?

the idea is to transmit the D (or $e[n]$) and let the receiver add that to the previous PCM sample it has (or can predict from the previous PCM samples using the same predictor).  and the reason was that, if the prediction was good, then $e[n]$ would be small and a smaller word width was needed.

  And by the definition "predicts the present samples with previous samples" how does the filter do that?

the idea is that a good prediction of $x[n]$ is a linear combination of $x[n-1]$, $x[n-2]$, $x[n-3]$, $x[n-4]$... $x[n-N]$.  so some guys in the sixties just suggested that the prediction is:

$$ \hat{x}[n] \triangleq a_1\, x[n-1] \ + \ a_2\, x[n-2] \ + \ a_3\, x[n-3] \ + ... + \ a_N\, x[n-N]   $$

the prediction error is

$$ e[n] \triangleq \hat{x}[n] - x[n] $$

and this $e[n]$ is what is transmitted and the next PCM sample can be constructed from the previous PCM samples and the error as

$$ \begin{align}
 x[n] &= \hat{x}[n] \ - \ e[n] \\
      &=  a_1\, x[n-1] \ + \ a_2\, x[n-2] \ + \ a_3\, x[n-3] \ + ... + \ a_N\, x[n-N] \ - \ e[n]
\end{align} $$

you can see that the prediction is a linear FIR filter and that the reconstruction has feedback and is a linear IIR filter (that happens to have all of the zeros at the origin, sometimes called an "all-pole filter" which is slightly imprecise).

then the idea is to, with knowledge of the statistical spectral characteristics of $x[n]$, namely its autocorrelation, find the values of $a_1$, $a_2$, $a_3$, ... such that the prediction is the best, that is such that $|e[n]|^2$ is minimized in a statistical manner.  this is what linear predictive coding (LPC) is at the most fundamental.  smaller $e[n]$ means that the word width needed for that signal is smaller and that the number of bits you need to send $e[n]$ over a transmission channel is less.  and, it means you can add Shannon and Huffman to the mix and send even fewer bits over the channel.  that is one way to do lossless coding.


For this you can use local features like SIFT, SURF, FREAK, etc.


You're almost there; you just need to connect a few dots.  Let your $\frac{10}{5 + i 2 \pi f} = G(f)$.  Then $g(t) = 10 u(t) e^{-5 t}$.

Now we get into that exponent part.  Your $h(t) = g(t - t_0)$ is correct, but you're applying it incorrectly.  You need to apply it to the part that I've labeled as $g(t)$: $h(t) = g(t - t_0) = 10 u(t - t_0) e^{-5 (t - t_0)} = 10 u(t - 4) e^{-5 (t - 4)}$.

So, basically, it's the output you'd expect from the "pure" transfer function, just delayed by 4 time units (presumably seconds).

Does that help?


Creating a matrix of Gaussian wavelets at dyadic scales

In this case the problem was that I was working in the frequency domain, but implementing an algorithm in the time domain.

I've removed the pre-multiplications of the Fourier matrix in the time domain and everything works fine.


This is likely due to a truncation of your sine wave in the time interval you have chosen, such that there is not an integer multiple number of sine wave cycles chosen. I am guessing that your frequency view that you provide above is a subsequent measurement (DFT) from the resulting time domain waveform you have created, and not a view of your original frequency function you are using to create the time domain waveform (as the original function in the frequency domain would not have those artifacts!).

That said, the specific answer will depend on how exactly you create and repeat/extend the final time domain waveform from the fixed length IFFT result that you get, but below are some insights on the effects of waveform truncation in the time domain that may help you to understand how you may be getting such artifacts.

A direct and concise mathematical answer of frequency effects from waveform truncation can be derived by describing your fixed length (sampled) waveform in the time domain as an infinite length waveform multiplied that is not truncated by a rectangular (boxcar) window. Multiplication in the time domain is convolution in the frequency domain, so the result of a boxcar window truncation would be convolving the frequency you would get without truncation with a sinc function (or when sampled, the approximation of a sinc function).

I like to view the resulting effect by picturing the time domain waveform extending from minus infinity to positive infinity in time, made up of your specific waveform snippet (IFFT result) continuously repeating in time. For a direct DFT this analogy this is mathematically valid (meaning the DFT over the finite time interval would match the DTFT of the same sequence repeating it time (just that the DTFT would be a continuous function and the DFT would be samples of this function- see below). Repetition in the time domain over an interval T results in a FT with discrete impulses in the frequency domain spaced at 1/T.).

From this view it may be clearer to see the effects of truncation: if we had 1.5 cycles of a sine wave in your time interval, the sine wave would suddenly snap every 1.5 cycles in the repeating waveform, requiring higher frequencies to achieve this "snap" that takes place in time. In contrast if you had an exact integer number of cycles (for a samples system, samples 0 to N-1, where sample 0 would be identical to the sample N for a pure undistorted sine wave), then the FT result (and DFT and DTFT) would only exist at one frequency location with 0 elsewhere; as expected for a pure tone.

FT, DTFT, DFT Background

I reference FT, DTFT, and DFT above, so I included this graphic and primer below to help clear up any confusion for those less familiar.

What is shown is the FT (or CTFT), DTFT, and DFT for a pure tone (that in "truth" extends from minus infinity to infinity in time), along with a modulated waveform of finite bandwidth.  In the Fourier Transforms, the tone is shown in blue while the modulated waveform is shown in red, and $\omega_s$ is the sampling frequency.

The CTFT is the Continuous Time Fourier Transform, or Fourier Transform (FT) for short.

$$ X(\omega) = \int_{t=-\infty}^\infty x(t)e^{-j\omega t}dt $$

Observe:

continuous in time -> aperiodic in frequency

aperiodic in time -> continuous in frequency

The DTFT is the Discrete Time Fourier Transform; the time interval still goes from minus infinity to infinity, but the system is now sampled. It is a sampled FT, which create repetition in frequency.

$$ X(\omega) = \sum_{n=-\infty}^\infty x[n]e^{-j\omega n} $$

Observe:

discrete in time -> periodic in frequency

aperiodic in time -> continuous in frequency

The DFT is the Discrete Fourier Transform; the time interval is defined over a finite interval. There is an implied periodicity in time, meaning if you made a waveform by repeating this one in time from minus infinity to infinity, the DTFT of this periodic waveform would match the DFT at the DFT sample locations (see plot above).

$$ X[k] = \sum_{n=0}^{N-1}x[n]  e^{-j\omega_o k/n} $$

Observe:

discrete in time -> periodic in frequency (implied)

periodic in time (implied) -> discrete in frequency

*note the time domain and frequency domain waveforms are not really periodic, as they only exist over a finite number of samples- but I like to make this interpretation, as if you did repeat the waveform in both domains you would get the same valued discrete samples (but they would be continuous waveforms with zero values in between). This interpretation has been extremely helpful for me in working with mixed signal (analog/digital) systems.


Your question is fairly broad, let me answer it step by step.

First of all $H(s)$ is indeed called the transfer function, and is the laplace transform of the impulse response. It's useful for finding the poles and zeroes, which is what the fourier transform can't do alone. $H(\omega)$ is the frequency response of the impulse response, not really the "transfer function".

If the frequency response is what you want to attain, then you can use the frequency response (fourier transform). It does not exist for unstable signals though.

The relation between Fourier Transform, discrete Fourier transform and Laplace transform/z-transform:

Fourier transform and DFT are in principle the same. One is for discrete signals, and the other for continuous signals. Because of that there are some different properties like frequency domain aliasing. What both of these transforms do is decompose the signal into complex exponentials (easily converted into sines and cosines (or sinusoids with magnitude and phase), which is the original idea).

Laplace and z-transform are again in principle the same, z-transform is the discrete equivalent to Laplace. It's simply the Fourier transform/DFT with the signal multiplied by a varying exponential $e^{\sigma t}$, or in discrete case just $r^t$. Then you have multiple fourier transforms/DFTs relating to each different value of the exponential ($r$ or $\sigma$). Point of interests are the poles and zeroes, IE. the points where the impulse response in question  multiplied by an exponential has integral that is exactly infinite or exactly zero.

Note that when $r=1$, or $\sigma = 0$, the Laplace/Z-transform is equal to the Fourier transform. In other words, the Laplace/Z-transform contains the Fourier transform in it.

EDIT: I wrote this fairly quickly so please point out any errors.


Your question is a bit broad. Since you probably have Matlab (I hope), Matlab help is where to start, for instance Measurements and Feature Extraction. From there, you can start with signals available in Matlab, get some tutorial code, etc.


Adjust the power of a digital signal to a given value


Why does fine CFO cause rotation of the constellation?

CFO is modeled like this $x(t) = s(t) e^{j2\pi f_{\text{CFO}}t} $. This is a phase rotation only and does not effect the magnitude. Convince yourself of it by breaking $e^{j2\pi f_{\text{CFO}}t}$ up as it magnitude ($= 1$) and its phase ($=2\pi f_{\text{CFO}}$). This is all happening in the time domain. If you plot the CFO-free signal and the CFO-corrupted signal magnitude over the top of each other, you should see little to no difference.
This is different from plotting the CFO-free spectrum and the CFO-corrupted spectrum. From the Fourier properties, multiplying by a complex exponential in time will result in a shift in frequency, so the CFO-corrupted spectrum will be shifted by $f_{\text{CFO}}$.
Edit
After seeing the comments, I think I understand the confusion better now. I'm using the term "sampling" to refer to taking a continuous time signal and sampling it. In a receiver, you do it at the symbol period to get a nice clean constellation back.
From what I can tell, the OP uses sampling to mean something else (not wrong just something else). In OFDM, the receiver does a FFT, which is also sampling but in frequency. The CFO introduces a shift in frequency so when the receiver goes to do the FFT, it is going to sample at points where you are getting a component from the main subcarrier and also other subcarriers (ICI).
To answer the question, CFO does not effect the magnitude of the signal (only the phase) and it does effect the magnitude of the spectrum at the subcarrier sampling points.


Discrete signals notation underlined value

An arrow or an underline is used to indicate the value at $n=0$. Example,

$$y(n) = \begin{array}{ccccc}[&9 &8 &\underline7 &6 &5 &]\end{array}$$ means that,
$$\begin{array}{ccccc}&y(-2) = 9, &y(-1) =8, &y(0)=7, &y(1) = 6 &\mathrm{and} &y(2) =5\end{array}$$

In your case, it says that $$\begin{array}{ccccc}&X(0) = 1, &X(1) =-1, &X(2)=2, &\ldots\end{array}$$

Note: No underline or arrow is equivalent to underlining the first value.


Eigenvalues and eigenvectors of a 3D image Laplacian

This answer comes a little late, but I think that it's necessary to clear up some of the confusion about what the eigenstructure of a Laplacian is and how it is calculated.

First of all, it's important to stress that this is not about properties of the local kernel used for calculating discrete derivatives. Instead you have to understand the Laplacian as linear operator on a vector space whose elements are in your case the volumetric data sets.

That means applying the Laplace operator is a linear map that maps one vector (data set) to another vector (data set). In this context, the eigenvectors of the Laplacian are again vectors (data sets) in that very same vector space. So the answer to your question would be a set of volumetric data sets, in fact as many as the dimension of your vector space (i.e the number of independent voxels).

Let's consider a very simple example. Take a 1-dimensional image, i.e. a single row of pixels, and let's also only use very few pixels, namely 4. Then the two center pixels have to direct neighbours, while the first and last pixel only have one neighbour each.

$$1 - 2 - 3 - 4$$

With this pixel geometry, we can give the discrete Laplace operator's result for the two center pixels 2 and 3 as linear functions of the pixels' values:

$$l[2] = p[1] - 2*p[2] + p[3]$$
$$l[3] = p[2] - 2*p[3] + p[4]$$

The other two pixels 1 and 4 don't have enough direct neighbours to calculate the discrete second derivative. We could fix that by assuming that pixels 1 and 4 are direct neighbours, closing the topology to a circle and imposing what is called a circular boundary condition. Or we just take the second derivatives at the boundaries to vanish. Let's do both, but start with the cyclic boundary condition. So we have:

$$l[1] = p[4] - 2*p[1] + p[2]$$
$$l[4] = p[3] - 2*p[4] + p[1]$$

This map is linear and we can write it as a matrix equation, mapping the column vector $p := (p[1],p[2],p[3],p[4])$ to $l = (L[1],L[2],L[3],L[4])$ by multiplication with the matrix $M$.

$$ L = M \cdot p $$

We call this matrix the discrete representation of the Laplace operator, and for our case it is
$$M=
\left(
\begin{array}{cccc}
 -2 & 1 & 0 & 1 \\
 1 & -2 & 1 & 0 \\
 0 & 1 & -2 & 1 \\
 1 & 0 & 1 & -2 \\
\end{array}
\right)
$$

The eigenvectors of this matrix are
$$ v_1 = (-1,1,-1,1) $$
$$ v_2 = (0,-1,0,1) $$
$$ v_3 = (-1,0,1,0) $$
$$ v_4 = (1,1,1,1) $$
with the associated eigenvalues
$$ \lambda_1 = 4
, \lambda_1 = 2
, \lambda_1 = 2
, \lambda_1 = 0 $$

You may recognise these vectors as the basis vector of the discrete Fourier transform on this vector, and the eigenvalues as their discrete frequencies. This is true in general, and in fact the decomposition of a vector (or more generally, a function) into the eigenspectrum of a Laplace operator generalises the idea of the Fourier transform.

Now let's investigate what happens if we use the alternative boundary conditions where $l[1] = 0$ and $l[4]=0$. The Matrix $M$ is then

$$M=\left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 1 & -2 & 1 & 0 \\
 0 & 1 & -2 & 1 \\
 0 & 0 & 0 & 0 \\
\end{array}
\right)$$

Cearly, this matrix has a different set of eigenvectors and eigenvalues. They're not as intuitive and insteresting, so I won't list them explicitly. It's worth noting however that we now get the eigenvalue $0$ twice, that means the eigensubspace where the laplacian vanishes is two dimensional.

So how do things change if we have a proper image instead of a single row of pixels? Not much. We only need to write down the Laplacian for each single pixel while considering the direct neighbour relationship, or topology, of the image. To make things a little trickier, let's go with an irregularly shaped two dimensional image.
$$
\begin{array}{ccccccccc}
    && 1 &-& 2 &-& 3 &&  \\
    && | & & | & & | \\
 4 &-& 5 &-& 6 &-& 7 &-& 8 \\
 | & & | & & | & & | & & | \\
 9 &-& 10 &-& 11 &-& 12 &-& 13 \\
    && |  & & |  &&  | & &  \\
    && 14 &-& 15 &-& 16 && \\
\end{array}
$$

Obviously now we have to take a 2-dimensional Laplacian by summing the second partial derivatives in horizontal and vertical direction. For that we require the two direct neighbours in each direction. The inner points $5,6,7,10,11,12$ therefore have a full 2-d laplace expansion. For point $5$ it looks like this for example:

$$l[5] = p[4] - 2 p[5] + p[6] + p[1] - 2 p[5] + p[10] \\= p[4] + p[6] + p[1] + p[10] - 4 p[5]$$

For the corner points $1,3,4,8,9,13,14,16$ we cannot construct a discrete second derivative so we use a boundary condition, e.g. $l[1]=0$

There are two points left, $2$ and $15$. Both have two direct neighbours in horizontal direction, but not in vertical direction. We can therefore apply a boundary condition that only affects the vertical direction, by setting the vertical second derivative to zero, while we evaluate the discrete second derivative horizontally and get $l[2]=p[1]-2p[2]+p[3]$ and likewise for $l[15]$.

Following this construction we get a linear equation for every point that relates it to the pixel values. Again, we write it as a matrix equation $L=Mp$, where the matrix in this case has $16\times 16$ entries. To be specific, it's
$$M=\left(
\begin{array}{cccccccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)$$

And again we can solve for the eigensystem of this matrix. A nice physical interpretation of the images that you get as eigenvectors is that they represent the vibrational modes of a membrane shaped like your image, with frequencies given by the eigenvalues.

You can easily step this game up to any number of dimensions, as long as you know the neighbourhood relations between your voxels. Simply formulate the individual linear equation like above, construct a matrix, find the eigensystem.

With the information gained from the eigenvectors of the Laplacian solving difference equations involving the discrete Laplacian can be greatly simplified. Once the eigenstructure is found, depending only on the geometry of the region, all data sets can be easily decomposed into the eigenbasis and the difference equations become trivial.


Low Pass Filter to Maintain Edge Information

There are many "Edge Preserving" filters in the image processing world.
2 very popular would be:

The Bilateral Filter.
Anisotropic Diffusion Filter.

Google search on each of the terms will give you plenty of data and code samples of each.

For any specific question about them, let me know, I'd be happy to assist.


$\tt fft$ vs. $\tt fftshift$ in MATLAB: Please explain

I believe you simply misunderstand what these functions are supposed to do. You get different results because they are not meant to do the same thing. I think it would be relevant that you read the documentation of these 2 functions.

Answers below:

  How do I decide which one to use?

fft : Use it when you need to compute the Discrete Fourier Transform of your time series
fftshift : Use it when you want to to swap the position of the negative and positive frequencies of a spectrum gotten with fft()

  Which is "correct"?

What do you consider correct ? fft will be correct if you need a frequency spectrum of your time domain signal. It doesn't really make sense to apply fftshift on a time domain series.

  Why are they giving such different answers?

Because they don't do the same thing.

  Why is the fftshift only real-valued while the fft is complex?

fft computes the discrete Fourier transform and by definition the output is complex.   fftshift doesn't compute anything except swaping the position of the samples, so if your input is real, you get real output.

Normally we use these two functions this way : spectrum = fftshift(fft(x))


Analyzing audio by looking at its time-frequency points in real- imaginary plane

This should be a comment but comments have limittaitns

I don't kno if this fixes your problem because I'm not sure what you are doing but using the handel file that comes with matlab

Modifying your code

figure(1);
ii=intersect(find(f>300 ),find(f<500 ));
for j=ii %I just want to look at these frequency indices
plot(real(Sl(:,j)),imag(Sl(:,j)),'Ob'); %grid on;

title(['F=', num2str(f(j)), ' hz']);

produced something.

Hope this helped


Continuous frequency grid to discrete grid

If I understood that right, you are trying to minimize the infinity norm over a discrete set of points. If that is the case you set over which you are trying to minimize is not convex and this is very much a NP hard problem, unless the set is very small, so that we can rigorously loop over it. Those are my thoughts on this.


GNU Radio | Async FEC encoder/decoder loopback | Issue during decoding?

The issue is that the FEC Async Decoder expects the input to be soft decisions with one decision per item instead of packed PDUs.

Here's a working example:

Encoder deployment:

Packed: yes
Rev. unpacking: yes
Rev. packing: yes

Repack Bits:

Packet alignment: Output
Endieness: LSB

Decoder deployment:

Packed: yes
Rev. packing: yes


Sound quality with really high sampling rate, but really low bit depth with dithering?

well, there is a sorta information-theory theorem regarding this (called the Gerzon-Craven limit) that was published long ago by Michael Gerzon and Peter Craven back in the late 80s.  i can't get a direct link to the paper reference, but type in "Optimal Noise Shaping and Dither of Digital Signals" into the AES Convention paper site.  this paper is a quarter-century old.  i just discovered that i have a camera copy of that old AES preprint and can send it to you if you send me an email to rbj@audioimagination.com .

i also have a copy of a paper making reference to the Gerzon-Craven limit On psychoacoustic noise shaping for audio requantization .

there is an interesting title regarding this called Approaching the Gerzon-Craven Noise Shaping Limit Using Semi-Infinite Programming Techniques. i don't have a copy of that and do not know how to get a free copy.

essentially, because of the Shannon Channel Capacity theorem, the answer to your question is "yes almost, theoretically", as long as you don't care about the magnitude of the quantization noise outside of your original 22.05 kHz baseband.  this is because the bit rate at 16-bit, 44.1 kHz is about the same (maybe just a teeny bit more) as the bit rate at 3-bit, 192 kHz and you can make an information-theory case that, because the information rate of the two cases are about the same, then the S/N ratio in the band of interest is also about the same.

but the real answer is "no, practically", because of the difficulty in constructing a 3-bit noise-shaped modulator with such extreme specs and that remains stable.  with a higher sampling rate than 192 kHz, the specs get less extreme.

if you have a working system (what's it coded in, MATLAB?), i would be interested in your noise-shaped modulator.  what order it is and how well it behaves.  sometimes these nasty high-order ΣΔ modulators behave badly in the presence of silence and we get, what are called, "idle tones".  if you think you have a working solution, i would be interested in hearing how it works for very low-level input.


This is the cross-correlation function of the two discrete-time random processes $x[n]$ and $y[n]$. It measures the similarity of the processes for a given time lag between the two. If $x[n]$ and $y[n]$ are jointly wide-sense stationary, then the cross-correlation only depends on the lag $k$, and not on the absolute time index $n$.

In terms of probabilities (or rather probability density functions), the expectation $E\{x[n]y[n-k]\}$ is given by

$$E\{x[n]y[n-k]\}=\int_{-\infty}^{\infty}xyf_{XY}(x,y;n,n-k)dxdy$$

where $f_{XY}(x,y;n_1,n_2)$ is the joint probability density function of $x[n_1]$ and $y[n_2]$.

In practice, we usually hope that the processes are ergodic, which allows us to estimate such an expectation by taking the appropriate time averages.


How to design very narrow, sharp lowpass filters - Only DC needed?

Your questions still leave me wondering as to what you're actually designing. For software implementation on modern x86 CPUs, CICs make almost no sense, but they are extremely elegant in hardware.

These filter definitions are ridiculous if you're planning to use a FIR – a transition width of 1mHz means that a minimum phase equiripple filter [1, (5.75), p. 142] would have

$$
N_e' \approx \frac 16 \log_{10}\left(\frac{1}{10\delta_1\delta_2^2}\right)\frac{f_s}{\Delta f}
$$
taps, with $f_s$ the sampling rate, $\Delta f$ the transition width, and $\delta_{1,2}$ the passband ripple and the stopband attenuation, respectively.

Now, your use case; I'm guessing from your the graphs from your question from 1 year ago, and also bearing in mind that you'd need incredible numerical accuracy to make use of > 80dB attenuation:

$$
\begin{align*}
f_s &= 6.4\cdot10^4\,\text{Hz,}\\
\Delta f &= 10^{-3}\,\text{Hz,}\\
\delta_2 &= -90\,\text{dB}\\
&= 10^{-9}\\
\delta_1 &= \delta_1\,\text{, everything else being wasteful either on suppression or precision side,}\\\\
N_e' & \approx \frac 16 \log_{10}\left(\frac{1}{10\cdot10^{-27}}\right)\frac{6.4\cdot10^4}{10^{-3}}\\
&=\frac16\cdot 26\cdot6.4\cdot10^7\\
&\approx 2.7 \cdot 10^8\,\text.
\end{align*}
$$

Obviously, Bellanger's approximation wasn't meant for such extreme cases, but the order of magnitude alone highlights the impossibility. Of course, you'd implement this as a decimating polyphase FIR, which means that you'd only need to apply
$$\frac{N_e'}{\frac{f_{s,input}}{f_{s,output}}}=N_e' \frac{f_{s,output}}{f_{s,input}}= N_e'\frac{\frac12}{6.4\cdot10^4}\approx2110$$
tapped filter every input sample clock, but the filter length puts impossible constraints on your machine accuracy if you're using floating point, and very hard constraints if you're building this in hardware on the bitwidth of your samples. As you notice, it does make a lot of difference for what we're designing a filter here – Matlab is not hardware :)

Yes, you can split your filter into multiple stages, but that will usually only reduce the computational load by running "parts" of the filter at lower rates, and will not solve the precision/bit width problems.

So we're ditching FIRs here, unless you're willing to reduce the requirements on stopband attenuation or transition width.

On to IIRs:

There's an approximate formula for the order necessary to achieve given specifications, and that is [1, (7.18), p. 184]:

$$
N_{IIR}\approx 1.076\log_{10}\left(\frac{2}{\sqrt{\delta_1}\delta_2}\right)
\log_{10}\left(\frac{f_s}{\Delta f}\frac4\pi\sin\left(2\pi\frac{f_1}{f_s}\right)\right)\text,
$$
$f_1=f_2-\Delta f$ and $f_2$ being your pass- and stopband edge. Since you're lowest output rate is 0.5 Hz, I'm guessing $f_2=0.5\,\text{Hz}$.

Applying the same specifications as above,
$$\begin{align*}
N_{IIR}&\approx 1.076\log_{10}\left(\frac{2}{\sqrt{10^{-9}}10^{-9}}\right)
\log_{10}\left(\frac{6.4\cdot10^4}{10^-3}\frac4\pi\sin\left(2\pi\frac{4.99\cdot10^{-1}}{6.4\cdot 10^{4}}\right)\right)\\
&\approx 1.076\log_{10}\left(2\cdot10^{13.5}\right)
\log_{10}\left(\frac{2.56\cdot10^8}\pi\cdot4.899\cdot10^{-5}\right)\\
&\approx 1.076\log_{10}\left(2\cdot10^{13.5}\right)
\log_{10}\left(3992\right)\\
&\approx 1.076\cdot13.80\cdot3.601\\
&\approx 53.48
\end{align*}
$$

Note that this is but an estimate! Also notice that this is far from a design approach, it just gives us an idea of how complicated the IIR will end up being.

Now, realistically, let's round up to $N_{IIR}=64$, and then do the classical approach of directly calculation:

We'll consider elliptic filters, which are optimal in the sense that for any given order, and accepting a certain ripple (which is your DC estimate accuracy, here), they have the most narrow transition width [1, p. 183].

Now, designing these taps isn't a trivial feat; scipy gives me something that roughly looks like that:

Notice the enormous dynamic range of each tap set ($\approx 10^{18}$), and of the overall filter ($\approx10^{65}$). That filter is nothing that reaches its accuracy on a double precision IEEE754 float.

I propose a hybrid approach: take cascaded halfband filters to get down to something like $\frac1{512}$ of the input sample rate, and then apply an elliptic filter.

1 Bellanger, Maurice: Digital Processing of Signals: theory and practice, second edition, Wiley 1988.


Kalman Filter with Accelerometer with DC Offset

Q1: Without actually implementing it and seeing that it works, your formulation looks pretty good... except that I do not understand your final equation:
$$
z_k=\left[\begin{array}{cccc}
1 & 0 & 0 & 0       \\
0 & 0 & 0 & 0       \\
0 & 0 & 1 & -1       \\
0 & 0 & 0 & 0
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
0 \\
a_k^{obs} \\
0
\end{array}\right]
$$
That doesn't make sense because a) the output matrix is too big and b) $a_k^{obs}$ does not form part of the state, $a_k^{true}$ does. It should look more like:
$$
z_k=\left[\begin{array}{cccc}
1 & 0 & 0 & 0       \\
0 & 0 & 1 & +1       \\
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
v_k \\
a_k^{true} \\
\Omega
\end{array}\right] = \left[\begin{array}{cccc}
d_k\\
a_k^{true} + \Omega
\end{array}\right]
$$
Notice the change in sign to get the expression for the acceleration measurement to come out correctly.

Q2: What noise is driving your system? Is it just the acceleration?

Your state update equation should look more like:

$$
\hat{x}_{k+1}=\left[\begin{array}{cccc}
1 & dt & dt^2/2 & 0       \\
0 & 1 & dt & 0       \\
0 & 0 & 1 & 0       \\
0 & 0 & 0 & 1
\end{array}\right]\cdot\left[\begin{array}{cccc}
d_k \\
v_k \\
a_k^{true} \\
\Omega
\end{array}\right]
+
Gw_k
$$
where $G$ determines which part of the state is driven by the noise process $w_k$.  The dimension of $w_k$ (is it a scalar or a vector?) will determine whether, if more than one part of the state is being driven, whether those drivers are different for each part.

If the driver is just acceleration, then choose:
$$
G = \left[\begin{array}{cccc}
0 \\
0 \\
1 \\
0
\end{array}\right]
$$
assuming that the offset in acceleration is not changing.

Then the question of what the covariance of $Q$ is boils down to: how quickly does your acceleration change?  That will really depend on your precise system. In this case, $w_k$ is just a scalar.

As to the covariance of $R$, that will depend on how different your measurements are form reality. I'd need more information to make a stab at that.


Maximum NCO frequency and filter design for down-conversion of complex input signal

Given the input signal is complex, then this can possibly be down-converted with a single real NCO as in the block diagram as follows as long as the frequency distance from DC or the sampling rate to the digital IF carrier is greater than the half the signal bandwidth (meaning none of the spectrum crosses those boundaries). By doing so increases the challenge of the low pass filter for isolating the images, especially when the signal is close to the DC or sampling rate boundary.

This is intuitively explained through observation of the "unrolled digital spectrum", such that we extend the frequency spectrum to plus and minus infinity and know through sampling that the spectrum from DC to the sampling rate will repeat in every additional integer multiple of the sampling rate. And further for real signals, we known that the positive and negative frequencies must be complex conjugate symmetric. We see all this in the diagrams below for an arbitrary complex IF waveform that is above half the sampling rate. These diagrams show the frequency domain spectrums for the down-conversion process with complex NCO's as typically done (and what happens if we rotate the wrong way), and finally what the above block diagram would result for a digital baseband spectrum.
The top diagram would be the most typical approach, using a complex NCO (not as in the block diagram above, but a full complex multiplier requiring four real multipliers and two adds, together with the complex NCO). The second diagram is what would occur if the NCO was wired up incorrectly (Sin/Cos swapped), and the bottom diagram is what we could do with the simpler diagram above, and shows the greater challenge in filtering out the image (but otherwise quite feasible). Note that because the real NCO signal actually creates two tones in the unique frequency span of plus or minus half the sampling rate (or equivalently  DC to the sampling rate), we can indeed use the real NCO to down-convert any signal in the span from DC to the sampling rate (as long as the signal's bandwidth doesn't cross those boundaries as already noted). We see this directly with Euler's formula: $2\cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$; a single sinusoid consists of two "tones", defining a "tone" as an impulse on the frequency domain.

The Complex IF is multiplied in time by the NCO, so in the frequency domains these spectrums would convolve resulting in the shifted spectrums at complex baseband. I find some grasp this quicker by knowing first that each impulse shown in the frequency domain represents a single spinning phasor in the time domain (given by $e^{j \omega_o n}$), so if you imagined an impulse at the center carrier for the complex IF given by $e^{j \omega_o n}$, you can see that it would be the one impulse out of many that is located at $e^{-j \omega_o n}$ that would move our signal to baseband, since in the time domain the product would be  $e^{j \omega_o n}e^{-j \omega_o n} = e^{j0}$, (and similarly you will see there are infinite combinations that do the same since the spectrums as shown repeat, but ultimately we know that everything from DC to fs repeats, so once we solve for all the signals in DC to Fs, we can simply populate the rest in this graphical approach). I highlight one example in the top diagram with the carrier at $f_c$ and the specific component in the complex NCO spectrum that is created at $-f_c$ that moves the complex IF to baseband through the time domain multiplication process. Through further observation we also see the same for the component in the fourth Nyquist zone just below $2f_s$, as we can find the component in the complex NCO spectrum that is the negative of this that would do the same (again this is just a visualization of the one digital spectrum that is unique over any frequency span of the sampling rate, such as $-f_s/2$ to $+f_s/2$ of $0$ to $f_s$, another approach is to view it as a cylinder but I find this clearer as can be represented easily on 2-dimensional plots). The math is very straight forward but given alone may mask this intuitive insight which helps in explaining many multi-rate and mixed-signal (analog/digital) principles.
A complex NCO will create a unique tone given by samples of $e^{j 2\pi f t}$ for any frequency from DC to the sampling rate (the sampling rate is an alias of DC, so to be accurate it would be unique from DC to one frequency spacing (given by $f_{\Delta} = f_{clk}/2^{acc}$ where $f_{clk}$ is the NCO clock rate and $acc$ is the accumulator bit size.). This would require an NCO with two outputs representing the real and imaginary components. For implementation this can be done with two real NCO's synchronized in quadrature given the relationship:
$$e^{j 2\pi f t} = \cos(2\pi f t) + j\sin(2\pi f t)$$
Implemented with "I" and "Q" where the $j$ is implied in the "Q" path. Ultimately the sinusoidal look-up-table would be shared with a simple offset in the phase accumulator of a quarter cycle (and only a quarter cycle needs to be stored to implement a full cycle with efficient up/down counting and sign change).
A complex NCO frequency translating a complex input is shown in the block diagram in this post:
Recovering signal for psk
Often with a DDC implementation, the input is real and the NCO is complex with a product of real and complex with complex output. The positive and negative half spectrums both get translated in the same direction (either left or right) such that one will end up at baseband (if shifted to the right, meaning the negative half spectrum, a spectral inversion will result which is easy to correct). The remaining higher frequency component needs to be filtered out after the multiplication.


Resize High Resolution Camera Capture to Low Resolution Camera Capture?

If you are using a “single-shot” sensor you could simply just not debayer the raw channels.  That would cut your sampling without adding the extra overhead that is already baked into the image acquisition.  A benefit of this approach is that you could probably boost your frame rate at the same time.


Adding the signals in time domain(sample by sample addition) while ensuring they don't add up to saturate will give you the desired effect. operation is adding and normalizing(like arithmetic mean). in question it is mentioned that taking mean is not giving desired results, this puzzles me. if i understand the desired effect properly then, adding element by element in time domain and normalizing it should do.


Calculate the harmonics in the power supply of $50\textrm{ Hz}$ using FFT

First of all, you should sample at rate greater than 2 times the highest digitized frequency, not equal to.

Second of all, do you have an anti-aliasing filter?


[EDITED FROM DISCUSSION]

On the first order, your data looks like a decay with a positive origin on a small-valued range $[0.7 \; 0.49]\times 10^{-7}$, and very tiny fluctuations with respect to the area under the curve.

So from afar, your data is much closer to an almost constant function than to some putative oscillations. So the the zero, or DC-component, is important. Indeed, the amplitude of your maximum FFT is around $5.9 \times 10^{-7}$, which happens to match the average value of your data, if I use the trapezoidal approximation $(0.7+0.49)/2= 0.595$ (I hope your FFT is normalized).

For second order observations, you can try to remove the exponential-looking decay trend (by an exponential fit, or a low-degree polynomial, of 2nd or 3rd order for instance). With the sufficiently dense number of points, you can also try again on the derivative (taking care of end-points) or second derivative, to see if periodic phenomena  appear. But this is likely to amplify noise.

At the third order, on the index interval $[0,\,2000]$, I see a kind of dropout around an order of magnitude of $0.05 \times 10^{-7}$. Once you remove the trend, this could appear like a single step down, which may cause a low-frequency effect in your signal.


Cascading conventional filter with matched filter

No, it is not “always” better.

One would  prefilter in circumstances where it makes sense, like Eckart filtering where the signal is random and the background has a known fixed component.

Filtering, generally introduces correlation in a time series.  A simple matched filter is derived under the assumption of a white background noise. If the prefilter is a whitening filter, it is better but each situation requires its own analysis.

Food for thought, if both filters are LTI, a prefilter can be a post filter.


$\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1} f(am,bn)e^{-j \cdot 2\pi(\frac {k}{M} \cdot m + \frac {l}{N} \cdot n)} $  (let $m'=am, n'=bn$)

$=\sum\limits_{m'=0}^{(M-1)a} \sum\limits_{n'=0}^{(N-1)b} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')}$

($m'$, $n'$ with the incremental interval of $a$ and $b$ respectively)

Note that
$\sum\limits_{m'=0}^{M-1} \sum\limits_{n'=0}^{N-1} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')}$
$=F(k/a,l/b)$

Here $m'$, $n'$ with the incremental interval of $1$, and the points $m$ and $n$ are evenly distributed on the circle, according to the symmetry, you have

$\sum\limits_{m'=0}^{(M-1)a} \sum\limits_{n'=0}^{(N-1)b} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')} = \frac {1}{|ab|}F(k/a,l/b)$


Why is the close looped model under control called plant?

I have heard that the word comes from "plant" as in "steam plant".  But I don't know of any historical research that's been done to pin down the term.
I did just check my copy of "On Governors" by James Clerk Maxwell himself (that guy got around).  At least in 1868, Maxwell wasn't using the word "plant" for "thing being controlled".


Double check z-transform ROC of $a^nu[n]$

You're right, the sum doesn't converge for values $z$ with $|z|\le \frac12$, that's why the equality only holds for $z$ inside the region of convergence ($|z|>\frac12$).


How do I get the phase angle from Cross Correlation?

There is no such thing as the phase angle between two signals unless they both consist of a single sinusoid at the same frequency, that is, $x(t) = A\cos(\omega t+\psi)$
and $y(t) = B\cos(\omega t + \phi)$.  If you have $N$ samples of these signals $x(t)$ and
$y(t)$, taken at times $0$, $T$, $2T, \ldots$, $(N-1)T$, so that
$$x[n] = x(nT), ~~ y[n] = y(nT), 0 \leq n < N,$$ and either $N\omega T$ is an
integer multiple of $2\pi$ or $N\omega T \gg 1$
then the phase angle between the two sinusoids is
$$\theta
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n]y[n]}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}\right).
\tag{1}$$
For complex-valued sinusoids $x(t) = Ae^{j(\omega t+\psi)}$ and
$y(t) = BAe^{j(\omega t+\phi)}$, $(1)$ should be replaced by
$$\theta
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|x[n]|^2\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}\right).
\tag{2}$$

Of course, these formula can be used for arbitrary signals, not just for pure sinusoids,
but then, what you get is the _angle between the two vectors
$\mathbf x = (x[0], x[1], \ldots, x[N-1])$ and
$\mathbf y =(y[0], y[1], \ldots, y[N-1])$ in the $N$-dimensional spaces
$\mathbb R^N$ or $\mathbb C^N$, and not a phase angle between two sinusoids
at the same frequency.
Note that the three points $\mathbf 0 = (0,0,\ldots,0)$, $\mathbf x$ and
$\mathbf y$ lie in a (two-dimensional) plane in $N$-dimensional space and
the $\theta$ that you get is the angle between the line segments with endpoints
$\mathbf 0$ and $\mathbf x$ and endpoints $\mathbf 0$ and $\mathbf y$ which
lie in this plane. Another way to think about this is that
$$\langle\mathbf x, \mathbf y\rangle = ||\mathbf x||\cdot||\mathbf y||\cdot\cos(\theta)$$
and thus $(1)$ is obtained from
$$\cos(\theta)
= \frac{\langle\mathbf x, \mathbf y\rangle}{||\mathbf x||\cdot||\mathbf y||}
= \frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|x[n]|^2}\sqrt{\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}.$$

If you must use FFTs because that's the way you have been told to do it, then you have
$$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} X[n](Y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|X[n]|^2\displaystyle \sum_{n=0}^{N-1}|Y[n]|^2}}\right)\tag{3}$$
so that you have the ineffable pleasure of not only needing to
compute two FFTs first, but also of using complex multiplications in $(3)$
instead of the real multiplications in $(1)$ (for real-valued
signals).  This is overkill in my estimation,
but as usual, YMMV, and what your boss insists on is always right, regardless
of what people write on Internet forums.


You can see for algorithms in ieee sites there are tons of algorithm which might suit your requirement and

An Algorithm for Real-Time Eye-Movement Tracking with Free Head Movement

Robust algorithm for video based eye tracing.
Are some examples.


I have created a C++ library that I haven't officially released yet that could be used for modulation.  It is called NimbleDSP.  It doesn't have a function or method called "modulation", but it does have interpolation and supports complex numbers, so it could easily be used to do modulation.

Liquid-dsp appears to be a mature SDR library, and it is used by the government framework RED HAWK that was recently open sourced.  My main problem with the library is that it uses the GPL license (may not be a problem for you) and the interface seems complicated.  Having said that, I looked at it fairly briefly so I could be wrong about that.

There is also SPUC (Signal Processing Using C++) and Aquila.  SPUC is also GPL I believe, but looks to be simpler than liquid-dsp while still being very comprehensive.  Aquila has a very permissive license but has some odd holes in its functionality.


K-means for 2D point clustering in python

It can be done very easily with the scikit-learn.  Examples are easy to find on their website, i.e. here. In my opinion it is the best way to go.

Modified code example from the above link:

import numpy as np
import matplotlib.pyplot as plt

from sklearn.cluster import KMeans
from sklearn.datasets.samples_generator import make_blobs

##############################################################################
# Generate sample data
np.random.seed(0)

batch_size = 45
centers = [[1, 1], [-1, -1], [1, -1]]
n_clusters = len(centers)
X, labels_true = make_blobs(n_samples=3000, centers=centers, cluster_std=0.7)

##############################################################################
# Compute clustering with Means

k_means = KMeans(init='k-means++', n_clusters=3, n_init=10)
k_means.fit(X)
k_means_labels = k_means.labels_
k_means_cluster_centers = k_means.cluster_centers_
k_means_labels_unique = np.unique(k_means_labels)

##############################################################################
# Plot result

colors = ['#4EACC5', '#FF9C34', '#4E9A06']
plt.figure()
plt.hold(True)
for k, col in zip(range(n_clusters), colors):
    my_members = k_means_labels == k
    cluster_center = k_means_cluster_centers[k]
    plt.plot(X[my_members, 0], X[my_members, 1], 'w',
            markerfacecolor=col, marker='.')
    plt.plot(cluster_center[0], cluster_center[1], 'o', markerfacecolor=col,
            markeredgecolor='k', markersize=6)
plt.title('KMeans')
plt.grid(True)
plt.show()

Yielding:


Which transformation in frequency domain equals a x-axis shift of a signal in time domain?

If $s(t)$ is your signal in time domain, you want to make the operation $s(t+\Delta t)$, which according to the Fourier's transform properties, is equivalent to the operation :
$$
S(f)e^{j2\pi f\Delta t}
$$
$S(f)$ being $s(t)$ Fourier's transform.

In case of a discrete Fourier's transform, which is defined as
$$
S[k] = \sum^{N-1}_{n=0}{s[n]e^{\frac{-i2\pi kn}{N}}}
$$
with $s$ the signal and N the number of points of both signal and DFT.

The operation to be applied in order to induce a shift by a number c of samples to $s[n]$ and obtain $s[n+c]$ is :
$$
e^{\frac{-i2\pi kc}{N}}S[k]
$$

Which in matlab code would be done with :

exp(-j*2*pi*(1:(N-1))*c/N).*S_fdt;


Peak location estimation for a triangle wave in AGWN?

You may be able to determine this in the frequency domain. The noise signal is broad band and has low spectral density at any one frequency. The triangle, however, has all energy concentrated in just a few frequencies, so these should stick out in a short term Fourier spectrum. You can match against your expected spectrum (odd harmonics, amplitude proportional to 1/n^2).

The harmonics of a triangle are all in phase. Since your triangle wave is arbitrarily located inside your analysis window, the harmonics will appear to have a linear phase. You can probably model your measured phase (after proper unwrapping) as a linear phase using a simple least squares error approach. This should give you the exact position of the peak of the triangular wave with respect to the beginning for your analysis window.


Simple, streaming, lossless image compression

You can consider using Huffyuv : http://neuron2.net/www.math.berkeley.edu/benrg/huffyuv.html

This is no great better than simple zip, but still slightly optimized for images.

Any Image related compression comes from the techniques like Vector quantization or Transform coding. In order to make use of transform such as DCT/Wavelet yet make it lossless you can think of JPEG-LS or JPEG2000 for compression. Only thing is, it is not streaming in your sense of definition.


For a given signal $g(t)$ with Fourier transform $G(\omega)$:

$$g(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} G(\omega)e^{i\omega t}\textrm{d}\omega$$

If you change variables such that $\omega=2\pi f$

$$g(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} G(2\pi f)e^{i2\pi ft}\textrm{d}(2\pi f)=\int_{-\infty}^{\infty} G(2\pi f)e^{i2\pi ft}\textrm{d}f$$

So if you present your Fourier transform as a function of $f$ instead of $\omega$, your Fourier pair of $g(t)$ will now be $\hat{G}(f)=\frac{G(2\pi f)}{2\pi}$ and:

$$g(t) =\int_{-\infty}^{\infty} \hat{G}(f)e^{i2\pi ft}\textrm{d}f$$

Note that with this definition, there is a scaling of $2\pi$ in comparison with the first equation I wrote. So now, for example, the Fourier transform of a cosine will consist of two deltas not with amplitude $\pi$, but with amplitude $\frac{\pi}{2\pi}=\frac12$, as correctly stated in your question.


How to quantize filter coefficients to 10 bits

Don't quantize the numerator and denominator coefficients of the total transfer function. Use the coefficients of the second-order sections - sos in Matlab - and quantize those. You will usually implement an IIR filter using second-order sections, especially in a fixed-point implementation.
The first thing to do is to scale the coefficients of the second-order sections (SOS). Note that you can scale each section individually. Scaling means multiplying the numerator coefficients of each SOS such that you minimize the chance of overflow in that section. The scaling depends on the chosen structure of the SOS, and on the filter coefficients. A detailed discussion on scaling can be found in the textbook Digital Filter Design by Parks and Burrus. If you're interested in the design and implementation of digital filter you should definitely read that book.
After having scaled the coefficients of each SOS, you quantize them by simply rounding them to the closest number on an equidistant grid consisting of $2^N$ numbers, where $N$ is the number of bits.
Example: The coefficients of one of the SOS are given by
b = [1, 1.621783996, 1]
a = [1, -4.030702997e-01, 2.332661953e-01]
You decide to scale the coefficients of that section by a factor of $\frac12$, which of course only affects the numerator coefficients:
bs = [0.5, 8.10891998e-01, 0.5]
You use $N=10$ bits to quantize the coefficients. Since the coefficients are in the range $[-1,1)$, we can multiply them by $2^{N-1}=2^9$, round the result to the nearest integer, and then multiply by $2^{-9}$ to scale back to the original range. In Matlab/Octave this is easily done as follows:

b = [1, 1.621783996, 1];
a = [1, -4.030702997e-01, 2.332661953e-01];
scale = 0.5;
b = b * scale;
N = 10;
range = 2^(N-1);
ac = round( a(2:3) * range ) / range;
bc = round( b * range ) / range;

The quantized coefficients are given by
ac = [-4.0234375e-01, 2.32421875e-01]
and
bc = [0.5, 8.10546875e-01, 0.5]
Or, in $10$ bit binary two's complement form
ac = [1100110010, 0001110111]
bc = [0100000000, 0110011111, 0100000000]


Only waveforms that have infinite support in one FT domain have finite support in the other.

This implies that only impossible (in the real world) or fictitious signals comply perfectly with the sampling theorem (at a finite sampling rate).  However given a non-zero noise floor (finite sized samples, quantum/Planck sampling time jitter, etc.), a common assumption is that the portion of the infinite spectrum sufficiently around or below the noise floor can be treated as zero, thus one can pretend that these signals are band-limited for practical purposes.


Butterworth filtering behaving unexpectedly -MATLAB

The effect you observed is mainly due to spectral leakage. It is a property of the DFT and it is only indirectly related to the Butterworth filter. Note that there are an integer number of periods of both sinusoidal components of y inside the DFT window. That's why you see two nice peaks in the frequency domain without any leakage. If you truncate y only by a few samples, you'd get spectral leakage even without the Butterworth filter, because now the period of the signal doesn't line up anymore with the DFT window:

plot(abs(fft(y(1:97))))

The output signal of the Butterworth filter shows some transient before it converges to its steady state. So it is not purely sinusoidal, and - more importantly - the number of periods of the strong sinusoidal components are not aligned with the length of the DFT window. If you remove the initial transient and if you make sure the DFT window aligns with the periods of the sinusoids you get two perfect spectral peaks:

plot(abs(fft(dataOut(11:100))))


How to detect "fast" changes in signal processing

The classic reference for this problem is Detection of Abrupt Changes - Theory and Application by Basseville and Nikiforov.  The whole book is available as a PDF download.

My recommendation is that you read Chapter 2.2 on the CUSUM (cumulative sum) algorithm.


Look up acoustic fingerprint or audio fingerprinting and dynamic-time warping for some basic technologies used to compare audio files.


Non-orthogonal basis which offers better energy compaction in image compression

Given a signal $x(t)$ of energy $\mathcal E$, suppose that $x(t)$
can be expressed as $\sum_i x_i\psi_i(t)$. If the $\psi_i(t)$ are a
complete set of orthonormal signals, then
$$\sum_i |x_i|^2 =\mathcal E$$
and one measure of energy compaction is the number of $x_i$ whose
squared magnitude is larger than some fraction of $\mathcal E$,
say $0.01\mathcal E$.  Another could be the smallest number of
$x_i$ such that $\sum_i |x_i|^2$ captures $99\%$, say, of the energy.

Some sets of orthonormal signals can
be much better for representing a specific $x(t)$ than other sets of orthonormal
signals. For example, sinusoids and complex exponentials are commonly
used to represent periodic signals, and do a good job with respect to
energy compaction when used with low pass or bandpass signals.  But they
require infinitely many components to represent signals such as
sawtooth waves that can be represented using just one or two coefficients
using orthonormal Legendre polynomials.
So why are generations of electrical engineering students forced to
calculate Fourier series of sawtooth waves at all?  Well, Fourier
series simplify life
in other calculations such as what happens when the signal passes through
a linear time-invariant system while Legendre polynomials do not
If signal-dependent orthonormal bases are allowed, then if we start
a Gram-Schmidt orthonormalization procedure beginning
with taking $\psi_1(t)$ to be the unit energy signal
$x(t)/\sqrt{\mathcal E}$ as the first member of a
set of orthonormal signals, then this orthonormal set can be used to represent $x(t)$ with
$x_1 = \sqrt{\mathcal E}, x_i = 0, i > 1$ for maximal compaction!  In fact,
we don't even need to bother to find $\psi_i(t)$ for $i > 1$.
So, how much compaction can be achieved with signal-dependent
bases depends on the degree of dependence allowed between the
signal and the basis to be used, etc.
The representation of $x(t)$ with respect to a set of non-orthonormal signals
is more complicated to determine.  With orthonormal signals, the coefficient
$x_i$ of $\psi(t)$ is easily determined since
$$\langle x(t), \psi_i(t)\rangle = \langle \sum_k x_k\psi_k(t), \psi_i(t)\rangle
= \sum_k x_k \langle \psi_k(t), \psi_i(t)\rangle = x_i,$$
since
$$ \langle \psi_k(t), \psi_i(t)\rangle =
\begin{cases}1, & \text{if}~k = i,\\0,  & \text{if}~k \neq i.,\end{cases}$$
while the simplification obtained via the orthonormality is not available
for non-orthonormal signals.  The question of how many of the $x_i$ are
needed to capture most of the energy is also difficult to determine
since the energy is not $\sum_i |x_i|^2$; it is necessary to take
the cross-products $x_ix_k$ into account as well.

So, unless there are specific
reasons why a particular non-orthonormal basis is important to use
(e.g. it helps with things one needs to do in image processing
or video processing or audio processing etc), simply using such a
basis because it provides better energy compaction might well be putting the cart
before the horse.  But, as always, your mileage may vary, and there may
be non-orthonormal signal bases that not only achieve better energy
compaction but also ease your life in other aspects.  It is just that
it is doubtful that a general theory exists saying that non-orthonormal
bases achieve better energy compaction than orthonormal bases.


Why does my sinusoid look "AM" in shape?

This visual phenomenon appears because the maximum frequency is close to the Nyquist frequency, or half the sampling frequency. Sampling begins to approach the limit of $2$ samples per period, and thus the linear interpolation performed by Matlab becomes highly inaccurate. However, samples are correctly located, as you can see from the code where an higher sampling ('Oversampled') is superimposed:

Fs=200e6;
Ts=1/Fs;
NFFT=2^14;
Runtime=(NFFT-1)*Ts;
t=0:Ts:Runtime;
f_in=90*1e6;
y_in=sin(2*pi *f_in *t);

Fs2=20*Fs;
Ts2=1/Fs2;
NFFT=2^14;
t2=0:Ts2:Runtime;
f_in=90*1e6;
y_in2=sin(2*pi *f_in *t2);

clf;hold on;
plot(t,y_in,'x')
plot(t2,y_in2,'-')
ylim([-1.5 1.5])
xlim([5.2 5.4]*1e-6)
legend('Original','Oversampled')


Assuming you have completed the basic course in DSP with a book like Oppenheim's Discrete-Time Signal Processing book then the next step could be anyone of the followings:

Multiresolution Signal Decomposition: Transforms, Subbands,Wavelets, A.Akansu

Multirate Digital Signal Processing, Rabiner

Multirate Systems and Filter Banks, P.P.Vaidyanathan

Wavelets and Filter Banks, Gilbert Strang

Multidimensional Digital Signal Processing, Dudgeon

Subband Adaptive Filtering: Theory and Implementation", by Woon-Seng Gan, Kong-Aik Lee, Sen M. Kuo

Wavelets and subband coding, Vetterli

Advanced DSP, John G Proakis

For a good book on multi rate DSP from communications perspective, please refer to the book below

Multi rate signal processing for communication systems, Friedrich J Harris

For a MATLAB specific handling, please refer to the book below

Multirate Filtering for Digital Signal Processing: MATLAB Applications, Ljiljana Milic

Personal recommendation the books by Vetterli and Fredrich Harris


Thanks to Dilip, I managed to figure out the answer to the question.
$y_1[n]$ will  be a normal random variable, with its variance given by the following equation:

\begin{equation}
\sigma_{y_1}^2=N\sigma_{th}^2
\end{equation}

This can be explained by using the fact that the variances will just add up. The mean will stay zero.$\\$
A similar approach can be used to find the variance of $y_2[n]$ ($\sigma_{y_2}$), which begins by writing the equation for $y_2[N-1]$ ($y_2[n]$ after $N$ cycles):

\begin{equation}
y_2[N-1]=(N-2)x[0] + (N-3)x[1] + (N-4)x[2]...+x[N-3]
\end{equation}

Taking variances on both sides:

\begin{equation}
\sigma_{y_2}^2=\sigma_{th}^2(1+2^2+3^2+...(N-2)^2)=\sigma_{th}^2\frac{(N-2)(N-1)(2N-3)}{6}\approx \sigma_{th}^2\frac{N^3}{3} (N>>1)
\end{equation}

Again, the mean will stay zero.


Need critical help: How to detect and distinguish two very similar looking signals?

As Conrad pointed out, a correlator is probably your best bet.

The correlation of a signal with itself (also known as its self-similarity) is larger than its correlation with any other signal (except for a constant factor related to the signals' energy).

In your case, you would implement two correlators, one for Signal 1 and one for Signal 2. Then, you'd look for one of the correlators' output to be larger than a certain threshold. You may have to run a few tests to determine the threshold.

One nice benefit of using correlation is that you don't need to synchronize your time window to the signals. Just be sure to accumulate the correlation calculations from one time window to the next.


First, you'll probably have better luck posting this on dsp.stackexchange.  That's a more specialized group that does stuff like this all the time.

In terms of your problem, here's a couple of options.

One is a machine learning approach.  e.g. create a training set by taking a bunch of data and hand marking the points that are good versus bad (like you've done above).  Then take an Artificial Neural Network, feed it the training data and teach it to recognize the pattern you want.

Another approach might be to take a Statistical process control approach.  Lots of tools there for detecting when something gets out of line.

Another way might a floating threshhold.  Compute the mean of the last X samples, and if the signal drops more than Y below that mean, call it bad.  could also compute the standard deviation of the last X samples, and if the signal drops more than Z std devs then call it back.

or, potentially just high-pass filter the signal to remove the slowly drifting mean and then you can use a fixed threshhold.


Are these Tx chain and Rx chains the correct way to implement this transmission?

No. When you want to transmit a real cosine in passband, then you need a single complex tone in baseband. Not a real one. Otherwise, you'll have two real tones on the air, one at the mixing frequency + cosine frequency, and one at the mixing frequency - cosine frequency.
But maybe, you wanted that, I'm not sure from your description. However, the fact that you're calling things "chirp" makes this unlikely.
Then, a real-valued sampled signal at a sampling rate of 60 MHz can't represent 56 MHz. The highest frequency it can represent is less than 30 MHz. Nyquist (The Nyquist-Shannon sampling theorem, to be exact)!
Also, you can't throw away the imaginary part in your receiver: you have no knowledge of your phase. So, you're unsynchronized in phase, and then you can't assume that the real signal you've sent is received as real signal. Especially if this, as I might guess, is about channel estimation or radar, you can't make that assumption.
So, I'm sorry but: you really need to refresh your education on sampling theory, and on what equivalent complex baseband and bandpass signals are, first!

I think I can't since the RF BW is 112 MHz (full BW) around the center frequency, but am I wrong?

Yes, the bandwidth is 56 MHz, not twice that.

Also, delete the file sink on the transmitter site: that signal you're transmitting is deterministic, and it's easier to just recreate it than read it from a file; in fact, that file sink hurts a lot, since you need to generate samples at a pretty high rate under hard real-time constraints. You're trying to write 60 MS/s · 32 bit/S = 1.92 Gb/s to a storage device; that takes computational power and bus bandwidth away from your PC, and acts as an artificial slow down.
I've seen few SSDs that can sustain that rate for a significant amount of time under these constraints.
On the RX side, as said, you need to be complex, so you need to write 60 MS/s · 2·32 bit/S = 2·1.92 Gb/s = 3.84 Gb/s to a storage device; that really is challenging. You might want to start with Vector sink first, which stores the samples in RAM, and use the "reserve memory" parameter of that to preallocate RAM for these samples; use a head block to cut off after you've collected enough samples and read the contents of the vector sink in Python. (That "reserve memory option" might only available in recent versions of GNU Radio – you're using an outdated one)


If you are applying the Short Time Fourier Transform (STFT) , then your data are more likely to look like a spectrogram, where the "height" is frequency and the "width" is time.

The..."legit" way to train an autoencoder type of network so that it learns to generate periodic-like sounds is to "teach" a bank of oscillators to "fire" in a predefined pattern. This "pattern" is what the VAE will "learn" from the training dataset.

By applying the STFT, you are segmenting your time into $N$ bins, where $N$ depends on the length of your STFT "window". Consequently, the STFT "window" also determines your frequency resolution. If you choose a low $N$, you will get very good temporal resolution but very poor frequency resolution too and therefore, this is a design choice.

Once you apply the STFT, you get the equivalent of a spectrogram. Taken a bit more liberally, the spectrogram gives you this information. Only that, in the case of the spectrogram, the "keyboard" is rotated 90 degrees clockwise (in a spectrogram the vertical axis is frequency and high pitch is upwards, while the horizontal axis is time with positive time extending to the right).

Every one of the keys is an oscillator. If you scan the spectrum at a given "time window" $t_n$ within a sound file and you find a high amplitude in that frequency then the key is "pressed" (the oscillator is active). Of course, the other parameter is how hard is the key pressed (amplitude).

So, your auto-encoder has to learn to predict the sequence by which it has to activate the oscillators and at what amplitude, in order to generate the sounds your are after.

This is the bigger picture. The "detail" to look for here is that the output of the STFT is $\mathcal{C}$omplex and you really need it to be to be able to take into account the phase that these oscillators are firing with respect to each other. Otherwise, you will be getting a baseline error below which your model would not be able to optimise any more.

More generally, I think that you would benefit from a little bit of background on speech (or music) generation by looking into what a vocoder is and how does it work, Prony's method by which signals can be synthesized from damped sinusoids (This is the closest to the piano analogy), speech synthesis and of course, the STFT and spectrogram links from above.

Hope this helps


Definitely don't use the whole file as a single frame.
As Alexander says, you want to use 25 millisecond frames with a 10 second overlap/hop (the vocal tract is stationary over this duration, so this is a common choice in speech recognition).
Each frame of audio makes a certain number of MFCCs (say 20). This number is constant for each frame of audio. But a file of different duration will, of course, have a different number of frames (and thus MFCCs). The size of your features (20) should always be uniform, but the number of observations (how many frames) you get from each file usually will not be, and that is okay.
Once you have extracted these acoustic-features, you need a machine learning algorithm that can "learn" what MFCC features are characteristic to a particular speaker. This is a very broad topic, but I can point to you this project: https://github.com/ppwwyyxx/speaker-recognition where you can see some of the machine learning algorithms used.

Good luck!


Mathematical expression of Multicarrier frequency equalizer

Hint: the convolution is not circular in that case as we used DCT instead of DFT.

Luckily, that's not true! At the receiver, you're doing a DFT, and you've prepared the transmit signal such that a convolution with the channel looks like a circular one by prepending the CP.
The fact that the data inside has something to do with the DCT is irrelevant there.

Then, the signal is recovered by multiplying by the DCT matrix

Yep, because that inverts $D'$.

I  think it's called frequency domain equalizer,

No, the frequency-domain equalizer is $\text{IFFT}(\text{FFT}/H)$.


Understanding Cb and Cr Components of YCbCr Color Space

YUV (or YCbCr) is like HSV, but in different coordinates. (The difference between YUV and YCbCr is marginal - mostly related to exact formulas).

The $V$ component is the same. $(S,H)$ can be thought of as polar coordinates, and $(U,V)$ as cartesian. $H$ is the angle and $S$ is the radius. A rough conversion would be:

$ U = S \cdot \cos(H) $

$ V = S \cdot \sin(H) $

You can see this link for more information.

Another thing to add to your intuition list:

Saturation is how pure the color is from spectral point of view. For example, a laser has a very narrow spectrum, which implies high saturation.


Principal Component Analysis (PCA) on Convolutional Neural Network (CNN) Features

I'm not into details of this specific case but I can see some logic.

A convolution layer can be reformulated as a Matrix Multiplication:

$$ y = W x $$

Let's say we trained on Data Set $ {x}^{1} $ which is big and general.
Namely we expect the trained weights $ {W}^{1} $ to be good enough for almost any other data set.

Let's assume we have another data set $ {x}_{2} $ with its matching weights $ {W}^{2} $.
By our assumption about how good $ {x}^{1} $ generalizes it means $ {W}^{1} \sim {W}^{2} $ with regard to some metric.

If this metric is the Frobenious Norm / Matrix Induced $ {L}_{2} $ Norm then using PCA (Which is basically SVD) makes sense as using the PCA we take the component of $ {W}^{1} $ which generalizes the most and use it for other case.


Estimation of Input Signal to Obtain the Desired Output Signal for an Unknown Filter

In the most simple case, just to give intuition about the problem, it is really easy.

In the Frequency Domain:

$$ {Y}^{\ast} \left( \omega \right) = H \left( \omega \right) {X}^{\ast} \left( \omega \right) \Rightarrow {X}^{\ast} \left( \omega \right) = \frac{ {Y}^{\ast} \left( \omega \right) }{ H \left( \omega \right) } $$

Since $ {Y}^{\ast} \left( \omega \right) $ is known all needed is $ H \left( \omega \right) $.

Yet since we have access to a black box of $ H \left( \omega \right) $ we can set input of a known signal $ X \left( \omega \right) $ and have the output $ Y \left( \omega \right) $ which will give us $ H \left( \omega \right) $:

$$ Y \left( \omega \right) = H \left( \omega \right) X \left( \omega \right) \Rightarrow H \left( \omega \right) = \frac{ Y \left( \omega \right) }{ X \left( \omega \right) } $$

As others pointed out it is called Deconvolution and as you mentioned it can be done using Least Mean Square (LMS) Filter.


Principal Component Analysis as a step between STFT and MFCC

One dimensionality reduction method (PCA) is linear, the other (MFCC) is  not. It  is thus unlikely that there is no difference at all since they probably do not commute. One serves more decorrelation, the other feature extraction.

The importance of the order thus depends a lot on your signal class, and what your purpose is.


Can every type of linear filter be modelled by a convolution?

No. It's only LTI (Linear and Time-Invariant) systems that can be modeled with convolution through a unique single impulse response.

For example the systems

$$ y(t) = g(t) x(t) $$ or
$$ y[n] = \sum_{k=0}^{k < n} x[n-k] $$

are both linear but not time-invariant and their output $y[n]$ cannot be computed with the convolution operation ( $\star$ denoting convolution)

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

as such an impulse response $h[n]$ does not exist.


How to measure the distance between two circle-like-lines over 360 degree

As @A_A said, it's a relatively straightforward problem once you've got good images.

Assuming you do have that:

Use the Hough transform to find the two main circles in the image
From the center of the inner circle, for each alpha, get a segment (by rotating your image and extracting a triangle) which should contain both your blue and black line
within that triangle, use appropriate signal processing (a simple "average per triangle column" or something, followed by a low pass filter, followed by a threshold, for example) to detect the exact position of the blue and black line


1) Group of points: $\{ (x_{i},y_{i})\}_{i=1..n}$

2) Centroid of the group of points: $\left(\frac{1}{n} \sum^{n}_{i=1}x_{i},\frac{1}{n}\sum^{n}_{i=1}y_{i} \right)$.

3) Center of the image: $(x_{c},y_{c})$

4) Translation of arbitrary point $(x,y)$ to the center $(x_{c},y_{c})$: $t(x,y)=(x_{c}-x,y_{c}-y)$

Note that $(x,y)+t(x,y) = (x_{c},y_{c})$ ...

Now just substitute 2) in 4) and you have your formula to translate the centroid to the image center. Right?

I am sure that I am missing something here - the translation distance obviously cannot be kept the same for all points. Didn't you mean another type of transform (maybe scaling have to be included as well) ?


Random Process at a particular time instance

Your equation (3) and (4) are wrong.

When you say

  Also, the sum of probabilities of a random variable should be 1

you actually mean the total sum of the values of probability mass function $P_{X_k}[x]$ equals one; i,e.

$$ \sum_{k} P_X[x_k] = 1 $$

(or the integral of the probability density funtion $f_{X_k}(x)$ equals one as
$ \int_{-\infty}^{\infty} f_X(x) dx = 1 $ )

But in your equation (3) and (4), you are not summing the values of probability mass function but the values of the sample function itself, which leads to te wrong conclusion and confussion.

To reach the correct conclusion, in your example RP, the pdf associated with each RV $X_{t_k}$ is
$$
f_{X_k}(x) = P(H)  \delta( x- \sin( \omega_1 t_k) ) + P(T) \delta(x- \sin( \omega_2 t_k))
$$

which integrastes (sums) to one as

$$
\begin{align}
\int f_X(x) dx &= \int \left( P(H)  \delta(x- \sin( \omega_1 t_k)) + P(T) \delta(z- \sin( \omega_2 t_k)) \right) dx \\
&= P(H) \int \delta(x- \sin( \omega_1 t_k)) dx + P(T) \int  \delta(x- \sin( \omega_2 t_k)) dx \\
&= P(H)+P(T) = 1
\end{align}
$$


You understand the basic operation of the phase detector. Reproducing some of the above, the detected phase error is:

$$
e_k = M \arctan\left(\frac{i_k}{q_k}\right) \mod{2\pi}
$$

As you noted, for an ideal QPSK constellation at the input to the phase detector, $e_k = \pi$ for all possible symbol values. However, in the presence of phase error $\theta$, the soft decisions at the phase detector input will have the form:

$$
\tilde{s_k} = s_k + \theta, \ \theta \in \left[-\frac{\pi}{M}, \frac{\pi}{M}\right]
$$

In this model, $s_k$ is the value of the QPSK symbol that is closest in phase to the observation. If you were making a hard decision on the symbol value based on observation $\tilde{s_k}$, this would be the one you would choose. Because the $\arctan$ function directly calculates the phase of the observed symbol, this phase error simply adds to the phase detector output: $e_k = \pi + \theta$. This allows the phase detector to detect the phase error in each symbol, modulo $\frac{2\pi}{M}$.

The job of the loop filter is to turn this stream of error outputs from the phase detector into a control signal that will adjust the frequency of the NCO that is used before the sampler. The integrator in the loop filter that you showed combats any frequency error (which manifests as a component to the phase error that increases linearly each symbol period). The proportional part of the filter will handle any constant component to the phase error as well.

Carrier tracking loop bandwidth can have multiple interpretations: sometimes it can refer to the pull-in bandwidth of the loop. For example, how large of a frequency error can I have and still expect to achieve lock? If the frequency error is too large, then the loop cannot detect and/or correct for it properly.

As a degenerate example, consider a frequency error $\omega = \frac{\pi}{2T}$ radians per second: this means that the carrier would spin exactly one quarter cycle between successive symbols. The phase detector described here would not detect any phase error at all (beyond the constant $e_k = \pi$ component discussed above).


The output process is clearly zero-mean because the LTI system cannot add a mean to the zero-mean input process. The variance of the filtered process is given by
$$E\{Y^2\}=\frac{N_0}{2}\int_{-\infty}^{\infty}|h(t)|^2dt\tag{1}$$
which results in $3A^2N_0/2$.


Extract hand written text from US census images

I couldn't open your pdf (antivirus said it was suspicious) so i downloaded an image myself, I don't know if its the right one, but its also from a 1920 census

You could perform image dilation on the document until the text is gone. You can do multiple passes if needed. I just did one. We know dilation will make shrink the black sections (or the text) this essentially gives us just the background.

We then remove the background from the image. This can be tricky depending on the language (I used matlab, and had to make sure to keep the right range and units for the image to display properly) then you have a cleaner image. You should be able to run a ocr engine on the result

and here is a zoomed in section just so you can see the text quality doesn't degrade much with this operation

and of course the code. As I said it is done in matlab but openCv(in many languages), python, or just about any other image library will be able to perform a similar operation

im = rgb2gray(color_im);
bg = imdilate(im,strel('disk',3));
%we do a complement becasue the bg subtraction makes the text white and
%the background black. so we invert it
clean = imcomplement(abs(bg-im));

subplot(1,3,1); imshow(im); title('original')
subplot(1,3,2); imshow(bg); title('after imdilate with 3 pixel circle')
subplot(1,3,3); imshow(clean); title('original minus background')

and here is the image I used. Its the first page of the 1920 census https://archive.org/details/fourteenthcensus02unit

edit

I forgot you said it was handwritten images. so I got another sample, same exact code, new image

notice the detail even in very dark obscured places, as the zoomed in view shows. You would need to run some other filters, but removing the main background noise is accomplished using dilation


In computed tomography (CT), why is 'Inverse Problem of Radon transform' studied?

Note that the original mathematical back-projection method assumes that the value of the Radon transform is known for all lines. This is an infinite amount of information, so it is pure mathematics at that level.
Since practical tomography works with a finite number of values of the Radon transform, the inversion is ill-posed. This cannot be ignored due to the finite precision of any numerical measurement and to the impact of noise.
In addition, since exposure to X-rays is harmful, it is imperative to continue to seek and develop algorithms that achieve results with fewer X-ray exposures.


How to find Cutoff frequencies (f1 &f2) for a narrow band pass filter for a particular central frequency?

The formulas for cutoff frequencies of a bandpass filter are different for different filter designs and can be calculated from the transfer function $H(ω)$. For a "basic" series RLC filter
$$
H(ω) = {\frac {ωRC} {\sqrt{(1-ω^2LC)^2 + (ωRC)^2}} }
$$
and this filter can be made narrowband if implemented with a high quality factor (see later in this text), but there are better designs for narrowband bandpass filters.
The 3dB level in amplitude is a ratio $1\over{\sqrt{2}}$; for an RLC filter the cutoff frequencies are calculated by solving the equation
$$
{\frac {ωRC} {\sqrt{(1-ω^2LC)^2 + (ωRC)^2}} } = {1\over{\sqrt{2}}}
$$
which is reduced to a quadratic equation having two roots
$$
ω_1 = -{R\over{2L}}+{\sqrt{\left({R\over{2L}}\right)^2+{ω_c^2}}} \\
ω_2 = {R\over{2L}}+{\sqrt{\left({R\over{2L}}\right)^2+{ω_c^2}}}
$$
$ω_c = {\frac {1} {\sqrt{LC}}}$ is a central frequency, $ω_1$ is a lower cutoff frequency, and $ω_2$ is a high cutoff frequency of the RLC filter. Unlike your guessed arithmetic mean, the central frequency of a RLC filter is a geometric mean
$$
ω_c = {\sqrt{ω_1 ω_2}}
$$
Another useful formula defines the filter bandwidth
$$
BW_{LRC} = ω_2 - ω_1 = {\frac {ω_c} {Q}}
$$
where $Q$ is a quality factor of the series RLC resonant tank:
$$
Q = {\frac {ω_cL} {R}} = \frac {1} {ω_cRC}
$$
For the other filter designs, the formulas are different.
Notice the use of an angular frequency ω in the formulas, $ω = 2\pi f$, which is customary for a signal processing math. Angular frequency is measured in radians per second.
UPDATE on question edit
The practical approach to constructing the filter is in the opposite direction to what you are seeking: given the design requirements, as the critical frequencies, ripple levels, and the sharpness of the cutoff, the priority of these requirements and the bill of materials, the developer selects the filter type, as Bessel, Butterworth, Chebyshev, elliptic or else, and adjusts the parameters to satisfy the requirements. See the filter example in your question's reference
sos = signal.iirfilter(17, [50, 200], rs=60, btype='band',
                       analog=False, ftype='cheby2', fs=2000,
                       output='sos')

the second parameter, [50, 200], sets the lower (50) and high (200) cutoff frequencies. Notice also in passing that in the Wikipedia article on the Chebyshev filter you can read

The common practice of defining the cutoff frequency at −3 dB is
usually not applied to Chebyshev filters; instead the cutoff is taken
as the point at which the gain falls to the value of the ripple for
the final time.

Usually you do not select these frequencies arbitrarily. Your decision depends on what is the purpose of your filter. For example, when you design a filter for audio signal, you can set the lower cutoff 16Hz, and the high cutoff 20KHz, because it is the hearing range of humans, and the transmission or amplification of frequencies outside this range would only enhance distortions and waste the electric power.


DFT of complex signal using Goertzel algorithm in C

The Goertzel algorithm (which is really just an efficient way of calculating what amounts to a single DFT bin at an arbitrary location) is defined for complex input, just like the DFT. A real input signal is really just a special case where the imaginary part is equal to zero. As far as your last question is concerned, if your input signal is real, then its spectrum is conjugate symmetric. That means that the spectrum value at -4000 Hz is the conjugate of the value at +4000 Hz. If you want to differentiate between the two frequencies, you'll need a complex input signal.


You asked "how" without first asking "if", and, if so, under what circumstances.

If you add two completely unknown numbers A + B and get 100, can you un-mix the sum to the original two numbers A and B?

But if you know B = 7, you might have a solution.

In the receiver case, if you know the 2 spectrums are disjoint (or other properties, such as redundancy in modulation/coding, etc.), judging might be possible.  Or the answer might be statistical.


What is the effect of a larger window size on "temporal resolution"?

You have probably heard of the uncertainty principles that arise in physics, with the most famous of all being the Heisenberg's uncertainty principle. But the uncertainty principle applies to signals as well. Just as Heisenberg's uncertainty is described by the following inequality:
$$
\Delta_x\, \Delta_p \ge \dfrac{h}{4\pi}
$$
The above imposes a limit in the level of resolution in the measurement of a particle's position and momentum. In fact their relationship is inverse meaning that as the particle's position measurements precision increases, then its momentum measurement precision decreases or in other words the uncertainty in its momentum measurement increases.
Now in signal processing the goal is to extract information about the signal's structure (if any) usually in the frequency or time domain (it could also be in the spatial domain).
But as you correctly described and analyzed in your graphs the more precise we get in the time domain the more uncertain our estimation of its frequency spectrum becomes. You can see it intuitively in the ripple. The smallest temporal window $N=16$ results in the biggest ripple at higher frequencies, taking also a sizeable proportion in the signal's power spectrum, but as $N$ increases so does the ripple increasing our certainty about the signal's frequency spectrum structure.
$$
\sigma_t\, \sigma_f \ge \dfrac{1}{4\pi}
$$
where $\sigma_t$ and $\sigma_f$ are the standard deviations of the time and frequency domain's estimations, respectively.
This minimum threshold is called the Gabor limit, and the principle is called Gabor's uncertainty.
Note:
  What is meant by "a larger window provides less temporal resolution" is that some of the properties of the signal in the time domain may be less prominent whereas in the frequency domain we would be more certain about its properties and structure.
From a practical perspective this limit is also important because in addition to other factors it could help us choose an appropriate window for our particular application where it could be desirable to extract features from both domains, sometimes we care only for the time domain properties, sometimes for the frequency domain properties and other times for both (also sometimes implications can be made about properties in one domain from properties measured in the other domain).
In conclusion, the more a signal spreads in the time domain, the less it spreads in the frequency domain and vice versa. In analog signals, a classic example of this phenomenon is the a rectangular pulse in the time domain, which Fourier Transform is the sinc function (exhibits spectral leakage).