Spectral density interpolation: Applications?

TL;DR:
The basis of this response is that matrix valued PSDs imply vector valued quantities that are expressed as Tensors. The Tensor is explored, briefly, to provide a basis to think (and possibly search further, beyond this post, for) practical applications. Practical applications explored here include DWTI-MRI and Quaternion RGBs.

  Let $\Phi_1$ and $\Phi_2$ be two matrix-valued power spectral densities.

How does a Power Spectral Density ends up being matrix-valued?

By applying it over a Tensor field or other sequence of Tensors.

What is a Tensor? And why do I need one?

The answer to "What is a Tensor" varies wildly depending on who one dares to ask. That's because different specialities are interested in different aspects and applications of Tensors. From my point of view, I thought I would provide a "necessity" based 2 min definition to cover how do we get to that "matrix valued" part:

A "scalar field" type of signal, associates a single quantity ($x \in \mathbb{R}^d$ ) with some point in space. That is, independently of the dimensionality ($d$) of the scalar:

One dimensional signals such as a discrete sampled sound recording ($x[i]$) are represented as a sequence of real ($x \in \mathbb{R}$) samples of sound pressure at distinct points in time. In this case, all we are interested in is how does the Sound Pressure Level vary over time. Notice here that $x$ is assumed to be one continuous recording without "gaps". In this case, successive samples $i, i+1$ are assumed to have been obtained at some $T_s$ time interval called "Sampling Period" (and of course, $i$ is an integer ($i \in \mathbb{Z}$)) . That is, we don't need to convey that $x[i]$ contains successive samples. It is implied.
Two dimensional signals, such as grayscale images, are still scalars. In visible photography, a sensor (or film surface) records how does light intensity vary within its Field of View. Again, we associate a single measurement ($x[i,j]$) with a point in space ($i,j$).
Three dimensional signals, such as volumetric data coming out of structural Magnetic Resonance Imaging scanner, are still scalars. In this case, all we are interested in is the concentration of a substance ($x[i,j,k]$) in three dimensional space ($i,j,k$). This concentration implies that $x$ might be more dense at $i_1,j_1,k_1$ but less dense around $i_2,j_2,k_2$ as it happens with fatty tissue around muscles for example.

Three dimensional signals such as video are still scalars. That is, we associate light intensity $x[i,j,n]$ with a point in 2D space $i,j$ and a point in time $n$.

Four dimensional signals such as those coming out of functional Magnetic Resonance Imaging are still scalars. In this case, we associate the concentration of a substance $x[i,j,k,n]$ with a point in 3D space $i,j,k$ and a point in time $n$.

But, what do scalar fields look like and why are they not enough?

This is what a three-dimensional scalar field looks like:

And here is a more realistic one. I am showing this one because the effect is more pronounced.

What do we observe? Well, it's not smooth is it? Why is it not smooth? Because it is only composed by little cubes that represent sampled space without any regard for direction. I only know that a particular point in space had some value but I don't know where was this value coming from.

Contrast this previous image with something like this:

This skull looks like it is composed of tiny little surfaces that tell us not only that some point in space $i,j,k$ is occupied by some bone substance but also what the orientation of this substance, in space, is.

To be able to represent this at the point of measurement, we cannot use a scalar field. We need something different. We need something that can encode direction. We need a vector.

Now, someone might say here: Hold on, you can calculate the local slope of a scalar field by obtaining its Gradient. Which, by the way, is a vector field.

And, for some applications, that could be enough but that would also mean that we make a particular assumption about the underlying phenomenon. That is, that "things" flow from high pressures to low pressures. In other words, we would still sample a scalar field of pressure (or density) values and from that infer flow. What if these regions of "high pressure" don't move? (in reality). Or what if only a small amount of their "density" flows?

What if what we need, genuinely, is to be able to not only record the concentration of a substance $x$ at some point in three dimensional space ($i,j,k$) but also its flow towards some direction $\rho, \theta, \alpha$?

Such is the case with Diffusion Tensor Magnetic Resonane Imaging. MRI is a beautiful "sensor" that can return a huge amount of information about a biological sample in its field of view.

It samples the human body directly in the frequency domain and to produce an "image" it transforms the acquired signals, from the frequency domain to the time domain.

Because of the way magnetic fields interact with magnetic objects (even as tiny as atoms), it is possible to selectively "excite" and "listen" for signals at a particular point in three dimensional space $i,j,k$.

This is how structural and functional MRI works which return scalar fields. You "excite" a particular region and then listen to the signal of atoms (from that particular region) "reacting" to the excitation (typically in the region of MHz, e.g. for H the frequency is ~43 MHz per Tesla). The more atoms you "listen" the higher the $x[i,j,k]$ at that point. Or more generally, $x$ is proportional (even if negatively) to the concentration of specific atoms at that $i,j,k$

But, by modifying the way atoms are excited you can also infer if and how did they move. How do you do that? You "excite" a particular region but "listen" from another region in quick succession. In other words, you "tag" some atoms, wait and then look where did the tag go (but very quickly because the "tag" is "wearing out" with time).

The geometrical orientation of these "excite" / "listen" regions now forms a vector. If you found "a lot of tagged atoms" towards a particular direction then this direction has higher magnitude of "flow", otherwise it does not.

So, in this case, we need to represent both the magnitude of the quantity ("How much substance moved?") and its direction ("In what direction did it move?").

We do this using a Tensor Field. It looks like this:

$$A[i,j,k] = \begin{bmatrix}
 A_{ii} & A_{ij} & A_{ik} \\
 A_{ji} & A_{jj} & A_{jk} \\
 A_{ki} & A_{kj} & A_{kk} \\
  \end{bmatrix}
$$

So, imagine a three dimensional grid such as this one:

This grid has a "node" (a point where three axes meet) at distinct locations $i,j,k$. But this time around, instead of having a single number associated with that $i,j,k$ we have a Tensor ($A$).

The MRI populates this matrix by obtaining six sets of measurements corresponding to "geometrical displacements" (i.e. directions) $A_{ij}, A_{ik}, A_{ji}, A_{jk}, A_{ki}, A_{kj}$. Pairs of $indexVar_1 indexVar_2$ denote directions. So $ij$ can be left to right, $ji$ right to left, $ik$ top to bottom, $ki$ bottom to top, $jk$ front to back $kj$ back to front WITH RESPECT TO $[i,j,k]$.

It is important to understand that $A$ is a matrix and it is one thing. Just as some vector is a "packet" (tuple) of two numbers ($\vec{a} = (b,c))$ but we treat it effectively as one thing, in exactly the same way, $A$ is one thing. One, indivisible, inseparable object. ONE thing.

This is important for what is about to follow.

Now, if you check back our $A$, you will notice that it contains a huge amount of information. So much in fact, that we can infer flow towards any direction.ANY direction. No, seriously, look back at $A$ again. $A$ (and of course, not $A$ but the MRI itself) captures "symmetric" flows. It captures flow from left to right and at the same time it captures a flow from right to left. Therefore, we can "see" exchange phenomena. Something travelling from left to right (perhaps mostly) but at the same time another flow (perhaps a low amount) travelling backwards in the same direction at the same site.

Most of the times however what we are interested in is just one vector out of all the possible directions that $A$ can encode (The vector of "most flow"). We extract this via the eigenvectors and eigenvalues of $A$.

The last step in this process is to visualise the flow. In this case, we simply "follow the arrows" to reconstruct an orbit. This orbit coincides with neuronal axons (in the brain) and this is what the DTI images depict. About this image, it is important to understand that this is a static image. This is a snapshot of the tracts, although it has been created by an extremely large amount of measurements that may have taken some time. If we did DTI over time, then we would have one $A$ associated with each $i,j,k,n$.

OK, that's fine, but how do you end up with PSDs that look like matrices?

Remember that $A$ that is one thing? Just like some $x[i]$ (e.g. some integer number) is one inseparable thing?

Find the periodogram of $A[i,j,k]$.

This is where it pays off to think about $A$ as one thing. Because you don't treat it differently than any other vector, although, yes, it certainly is not "just a vector".

Now, I will use Welch's method of deriving the periodogram here which specifies splitting the "signal" into overlapping segments, taking the DFT of those segments, deriving the spectrum and then averaging those segments to derive the final power spectrum.

BUT!!!

The most challenging thing here is NOT the computational part. The most challenging thing here is the connection with reality. What does it mean? What does it look like? In fact, the most challenging thing in the aforementioned process of deriving the periodogram is not the "indexing" (chopping the Tensor field down to overlapping segments) because that's just like slicing a scalar array. Only that now, each one of the entries of this Tensor field is a Tensor. So, this would look like an Array of Tensors, which we cut down into overlapping segments. It is not the averaging, because, averaging is averaging. Tensor addition works as you would expect, you just add the individual components.

The most challenging part, in terms of comprehension is The Discrete Fourier Transform. Again, it's not about calculating it, because, the Discrete Fourier Transform of a Tensor field is the Discrete Fourier Transform of the individual components of the Tensors that compose it. In terms of notation, you still have sinusoids but they are packed in a Tensor "package" this time.

The challenging bit here is the meaning of it and what does it look like given that it describes so much information.

One (useful) view would be to treat it like the DFT of complex data. Because in that case, it's not $x \in \mathbb{R}$ anymore but rather $x \in \mathbb{C}$ and therefore $x$ is now actually a vector (i.e., it has magnitude and direction).

When $x \in \mathbb{C}$ then what the Fourier transform actually does is decompose $x$ to the sum of a set of "squiggly" Lissajous symbols. That is, the complete repertoire of "doodles" that two interacting (orthogonal) sinusoids can create, including a perfect circle (at DC) but more importantly ellipsoids at various orientations.

This concept can be applied as is to the Fourier Transform of the Eigenvector of $A$ described above. In that case what we would be getting after $\mathcal{F}(\lambda \cdot v), A = \lambda \cdot v$ would be a decomposition of $v$ down to a set of "squiggly" three dimensional surfaces including a perfect sphere and also ellipsoids with various orientations (and everything in between).

I suppose that something similar can be obtained to interpret $\mathcal{F}(A)$ although it probably does not look exactly the same because each frequency of the PSD is now a complex matrix and returns a huge amount of directional information about that particular point (which we kind of "discard" when we focus only on the eigenvector).

  ...a practical application in which the considered spectra are matrix-valued functions

This has been partially covered by the DTI and discussion in the previous section. The only thing that I would like to add is that in terms of interpolating, you can interpolate incredibly complex curves in surprising ways with results similar to this one.

And also that in the signal / image processing space, there is a "Tensor" expression of the DFT that would return spectra you can work with but I am not really sure about interpolation.

Perhaps you can define some sort of "smoothing" filter. That is, given a Tensor field with missing spectral values, fill them in by interpolating the known values around it.

Other than this, I suppose that you can define a metric of "distance" between two spectra (?) to "match" tensor fields.

Beyond these applications in DSP, the next stop would be the Physics SE. They use Tensors a lot in cosmology. I don't quite understand this stuff,...I mean, what does it look like...and sometimes the "why". I guess they like tensors because they express that oneness required to describe space and time...but how do you go from that to, "maybe "time" does not exist", is impenetrable for me.

Anyway, this is why I asked "what is it that you do"? What is your "necessity" for ending up with Matrix valued PSDs? What is the "problem"?

Hope this helps.


This question seems to be based on several misconceptions.

In a QPSK modulator operating at a rate of $N$ baud, two bits enter the
modulator during each $T = N^{-1}$ second interval.  If the bits are
entering on a single wire (at a rate of $2N$ bits per second, the
QPSK modulator first converts this serial input
bit stream at $2N$ bits per second into two parallel bit streams at
$N$ bits per second. For example, alternate bits are steered to
alternate streams so that, say, all the odd-numbered bits
are in one stream and all the even-numbered bits are in the
other stream.

In some cases, the bits entering the QPSK modulator
are already formatted into two parallel bit streams at $N$ bits per
second. One important and commonly observed arrangement of this type
is when a rate-$\frac12$ convolutional encoder precedes the QPSK modulator,
so that information bits enter the encoder serially
at rate $N$ bits per second, and the convolutional encoder spits out two
parallel bit streams at $N$ bits per second.  However, regardless of
encoding or not, the bits that the QPSK
modulator uses are in two parallel streams that are synchronous
(that is, the bit transitions occur at the same times (once every
$T$ seconds) and in each
bit stream, the bit occupies the entire $T$ second interval.

As the OP's question says, the bits are modulated onto cosine and
sine carriers.  As discussed in more detail here, if the bits in the two streams are denoted by $b_I$ and
$b_Q$, then the corresponding signals, say, over the interval
$[0,T)$ can be taken to be
$$(-1)^{b_I}\cos(2\pi f_c t) ~~ \text{and}~~ -(-1)^{b_Q}\sin(2\pi f_c t),
0 \leq t < T. \tag{1}$$
(Ignore the negative sign on the $\sin$, it is there to simplify
notation in more detailed analyses).  The QPSK signal is the sum
$$(-1)^{b_I}\cos(2\pi f_c t) -(-1)^{b_Q}\sin(2\pi f_c t),
0 \leq t \leq T. \tag{2}$$
It is not the case, as the OP seems to think that the two
signals are separated into even and odd bits, so that
$(-1)^{b_I}\cos(2\pi f_c t)$ lasts from $t=0$ to $t=T/2$ and
$-(-1)^{b_Q}\sin(2\pi f_c t)$ from $t=T/2$ to $t=T$: the
cosine and sine signals last for the full $T$-second duration
and the bits $b_I$ and $b_Q$ both affect the signal phase for all
$T$ seconds.

Offset QPSK also converts the serial bit stream into
two parallel bit streams in which the bit intervals are
of duration $T$ on each wire, but the bit transitions
are not synchronous but occur with an offset of
$T/2$ seconds.  Thus, $b_I$ modulates the cosine carrier
to produce
$$(-1)^{b_I}\cos(2\pi f_c t), 0 \leq t < T$$
but $b_Q$ modulates the sine carrier to produce
$$-(-1)^{b_Q}\sin(2\pi f_c t), \frac{T}{2} \leq t < \frac{3T}{2}.$$
The modulator output is still the sum of the two sinusoids,
but the sinusoids change phase at different times. In fact,
the signal can be deemed to be a $2N$ baud signal since the
phase can change once every $T/2 = (2N)^{-1}$ seconds.
Each bit still affects the phase for $T$ seconds, but the
bit transitions are staggered.  In the receiver, integrations
(or matched filterings)
are carried out over (staggered) $T$-second intervals in the
two branches of the receiver, and the decision devices
spit out one bit every $T$ seconds but at staggered times.
It is easy to multiplex these decisions onto a single wire
that thus carries a bit stream at $2N$ bits per second,
just as the input to the QPSK modulator (if anyone still
remembers that far back in this answer) is a bit stream
on a single wire at $2N$ bits per second.


The reason is g1 is not a complete Gaussian pdf since you limited the range to not extend past 0.

The following will work for you:

W=[-2500:1:2500];
g1= normpdf(W, 500, 100);
g2= normpdf(W, 0, 100);
ConvG1G2 = conv(g1,g2, "same");
plot(g1)
hold on
plot(g2)
plot(ConvG1G2)

Results in:


Plotting a Time Frequency contour/colormesh plot of a Discrete Wavelet Transform

With wavelets, that rely on scaling, you can expect time-scale, more than time-frequency. Obtained 2D images of scalograms from continuous wavelets is relatively easy. Most of them are not exactly invertible.
Exactly invertible discretized wavelets come in different shapes. To preserVe images with equal-sized pixels, you could prefer redundant discrete wavelets over DWT.
So, maybe rwt.rdwt()


You can force denormal numbers to zero by setting the "flush to zero" and "denormals are zero" bits (15,6) in the MXCSR register.

_mm_setcsr( _mm_getcsr() | (1<<15) | (1<<6));


What is the distinction between ergodic and stationary?

A random process is a collection of random variables, one for each time instant under consideration.  Typically this may be continuous time ($-\infty < t < \infty$) or discrete time (all integers $n$, or all time instants $nT$ where $T$ is the sample interval).

Stationarity refers to the distributions of the random variables. Specifically, in a stationary process, all the random variables have the same distribution function, and more generally, for every positive integer $n$ and $n$ time instants $t_1, t_2, \ldots, t_n$, the joint distribution of the $n$ random variables $X(t_1), X(t_2), \cdots, X(t_n)$ is the same as the joint distribution of $X(t_1+\tau), X(t_2+\tau), \cdots, X(t_n+\tau)$.  That is, if we shift all time instants by $\tau$, the statistical description of the process does not change at all: the process is stationary.
Ergodicity, on the other hand, doesn't look at statistical properties of the random variables but at the sample paths, i.e. what you observe physically. Referring back to the random variables, recall that random variables are mappings from a sample space to the real numbers; each outcome is mapped onto a real number, and different random variables will typically map any given outcome to different numbers.  So, imagine that some higher being as performed the experiment which has resulted in an outcome $\omega$ in the sample space, and this outcome has been mapped onto (typically different) real numbers by all the random variables in the process: specifically, the random variable $X(t)$ has mapped $\omega$ to a real number we shall denote as $x(t)$. The numbers $x(t)$, regarded as a waveform, are the sample path corresponding to $\omega$, and different outcomes will give us different sample paths. Ergodicity then deals with properties of the sample paths and how these properties relate to the properties of the random variables comprising the random process.

Now, for a sample path $x(t)$ from a stationary process, we can compute the time average
$$\bar{x} = \frac{1}{2T} \int_{-T}^T x(t) \,\mathrm dt$$  but, what does $\bar{x}$ have to do with $\mu = E[X(t)]$, the mean of the random process? (Note that it doesn't matter which value of $t$ we use; all the random variables have the same distribution and so have the same mean (if the mean exists)).  As the OP says, the average value  or DC component of a sample path converges to the mean value of the process if the sample path is observed long enough, provided the process is ergodic and stationary, etc.  That is, ergodicity is what enables us to connect the results of the two calculations and to assert that
$$\lim_{T\to \infty}\bar{x} = \lim_{T\to \infty}\frac{1}{2T} \int_{-T}^T x(t) \,\mathrm dt$$
equals
$$\mu = E[X(t)] = \int_{-\infty}^\infty uf_X(u) \,\mathrm du.$$  A process for which such equality holds is said to be mean-ergodic, and a process is mean-ergodic if its autocovariance function $C_X(\tau)$ has the property:
$$\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^T C_X(\tau) \mathrm d\tau = 0.$$

Thus, not all stationary processes need be mean-ergodic. But there are other forms of ergodicity too. For example, for an autocovariance-ergodic process, the autocovariance function of a finite segment (say for $t\in (-T, T)$ of the sample path $x(t)$ converges to the autocovariance function $C_X(\tau)$ of the process as $T\to \infty$. A blanket statement that a process is ergodic might mean any of the various forms or it might mean a specific form; one just can't tell,

As an example of the difference between the two concepts, suppose that $X(t) = Y$ for all $t$ under consideration.  Here $Y$ is a random variable. This is a stationary process: each $X(t)$ has the same distribution (namely, the distribution of $Y$), same mean
$E[X(t)] = E[Y]$, same variance etc.;
each $X(t_1)$ and $X(t_2)$ have the same
joint distribution (though it is degenerate) and so on.  But the process is not
ergodic because each sample path is a constant.  Specifically, if a trial
of the experiment
(as performed by you, or by a superior being) results in $Y$ having
value $\alpha$, then the sample path of the random process that corresponds
to this experimental outcome has value $\alpha$ for all $t$, and the DC
value of the sample path is $\alpha$, not $E[X(t)] = E[Y]$, no matter how
long you observe the (rather boring) sample path.  In a parallel universe,
the trial would result in $Y = \beta$ and the sample path in that
universe would have value $\beta$ for all $t$.
It is not easy to write mathematical specifications to
exclude such trivialities from the class of stationary processes, and so this is a very minimal example of a stationary random process that is not ergodic.

Can there be a random process that is not stationary but is ergodic? Well, NO, not if by ergodic we mean ergodic in every possible way one can think of: for example, if we measure the fraction of time during which a long segment of the sample path $x(t)$ has value at most $\alpha$, this is a good estimate of $P(X(t) \leq \alpha) = F_X(\alpha)$, the value of the (common) CDF $F_X$ of the $X(t)$'s at $\alpha$ if the process is assumed to be ergodic with respect to the distribution  functions.  But, we can have random processes that are not stationary but are nonetheless mean-ergodic and autocovariance-ergodic.  For example, consider the process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$.
Note that each $X(t)$ is a discrete random variable that, in general, takes on four equally likely values $\cos(t), \cos(t+\pi/2)=-\sin(t), \cos(t+\pi) = -\cos(t)$ and $\cos(t+3\pi/2)=\sin(t)$,  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions, and so the process is not even first-order stationary. On the other hand,
$$E[X(t)] = \frac 14\cos(t)+
\frac 14(-\sin(t)) + \frac 14(-\cos(t))+\frac 14 \sin(t)  = 0$$ for every $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[\cos(t)\cos(s) +  (-\cos(t))(-\cos(s)) + \sin(t)\sin(s) + (-\sin(t))(-\sin(s))\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s).
\end{align}
In short, the process has zero mean and its autocorrelation (and autocovariance) function depends only on the time difference $t-s$, and so the process is wide sense stationary. But it is not first-order stationary and so cannot be
stationary to higher orders either. Now, when the experiment is performed and the value of $\Theta$ is known, we get the sample function which clearly must be one of $\pm \cos(t)$ and $\pm \sin(t)$ which have DC value $0$ which equals $0$, and whose autocorrelation function is $\frac 12 \cos(\tau)$, same as $R_X(\tau)$, and so this process is mean-ergodic and autocorrelation-ergodic even though it is not stationary at all. In closing, I remark that the process is  not ergodic with respect to the distribution function, that is, it cannot be said to be ergodic in all respects.


If I can make a call to another city, why does wifi only reach 50 meters?

This question, having a physical and network rather than mathematical nature, is close to being off-topic, but for your convenience let me put an answer for it.

Making a phone call is traditionally a communications network operation. But recently IP telephony over computer networks (internet) also emerged and furthermore GSM (mobile cellular) systems are so much intermingled with internet that the traditional medium seems to have blurred boundaries.

Making a traditional phone call using your GSM cellular phone, involves accessing an operator base station (BS) closest to your location. Base Stations, work at 900/1800 MHz, ranges and can be reached at as low as about -110 dB signal levels. Therefore, your phone can reach a BS about a 5/10 kilometers away in free space and less than a km in city areas. Once it's reached, remaning connection to other cities etc., is powered and handled by the BS network itself. Your phone is only required to reach the closest BS.

Wi-Fi is term used to denote wireless internet access point, that enables internet capable devices to wirelessly connect to internet. It may be at your home, office or airport, operating at 2.4 (or recently 5) GHz. It's signal output power level is low, therefore it can be reached within a hundred or less meters for fast / reliable connections. Once you access a wifi, the remaining communication is powered by Wi-Fi adapter / router device and your device is not involved.

As another option, you can also make an IP phone call using WhatsApp, Skype etc with your smart phone-computer combo. Then it's either accesing the closest (and allowed) Wi-Fi access-point at home/office or it may also be doing it through the GSM BS internet using its mobile data feature...


Positive slope of BER with a greater frequency deviation?

This could happen as discriminator gain is increased with a filter discriminator approach since in many of those approaches the gain would be maximum and linear for small signals only and then the slope of the discriminator slowly goes down coinciding with the results in your plot (such that you no longer get a perfect sine wave out for a sine wave in—- so you can also look at it as converting some of your signal components to other harmonics at the discriminator output due to that non-linearity). Ultimately as you increase the gain further your signal will start to saturate and then even invert depending on the wider bandwidth shape of your discriminator.

Plot the derivative of your frequency discriminator to see this more clearly as to its usable range for a given frequency deviation- the ideal for no distortion is to have a constant slope.

I see from the link you provided that the slope is indeed constant over a limited usable range. So specifically in your case I assume as $k_f$ was increased, it would cause a proportionate amount of the signal to go beyond the discriminator range (where your response is 0).


As far as the recovery time and results for the timing recovery, this is a loop implementation and you would need to review the complete loop for stability and gain parameters to balance loop bandwidth & convergence time, stability and noise performance.

Even with direct board to board communication you will have a static timing offset as the receiver needs to determine the optimum point in which to sample the recovered data in each symbol.

Step response testing of the timing error at your loop filter input and output can be insightful; the input to the loop filter will likely be quite noisy but you can filter your observation (as long as the filter bandwidth is much wider than your expected loop bandwidth this will be valid) to be able to see how it responds to a step change in time offset. Even better is if you have an accurate simulation model of your receiver so that you can observe the time offset more directly, and to observe what occurs as you change loop gain parameters.

Also through simulation and with testing on your hardware I would recommend creating an eye-diagram to better understand the core problems you may be having with data recovery. With this you can best observe the predicted sampling location for making a data decision versus the ideal.


First of all make sure that your $\delta_p$ value is linear (and not in dB). The function firpm needs a linear value. Don't worry about the passband error ripple having both signs; the value of $\delta_p$ used to compute the weighting function is just the maximum error magnitude, which is $0.2$ in your case. So you correctly defined the stopband weight by $0.2/0.05=4$.

The question is if you need to write an algorithm that automatically finds the necessary filter order, or if you can just try to find the order manually. If you need an automatic procedure you could use the second output argument of [h,err]=firpm(...), which is the maximum approximation error. You could design filters with increasing orders in a loop until the error returned as a second argument from firpm is smaller than the desired value for $\delta_p$ (because the approximation error equals weight times deviation, and since you chose the passband weight to be $1$, the error equals $\delta_p$).

Finally, note that the stopband edge equals the passband edge plus the transition bandwidth. So if the transition bandwidth is really $0.3$, as in your example, then the stopband edge would be the cut-off frequency (passband edge) plus $0.3$ (which equals $0.4$ in your example). These specifications are extremely mild and can be satisfied by a filter of very low order.


Does instability make an otherwise LTI system nonlinear (or time-variant)?

I don't know how the discussion arrived to this point but it is pretty convoluted to follow after this much of time. It is an ordinary difference equation with constant coefficients thus defines a linear, time invariant system. There is no need to pursue further from that if a unique solution exists. Here we have a problem with the emphasized part.

Let's first write the descriptor state space equations: the system is described by

$$
\begin{align}
\begin{bmatrix}1&-1\\0&0\end{bmatrix}
\begin{bmatrix}s_1[k+1]\\s_2[k+1]\end{bmatrix} &=
\begin{bmatrix}p&0\\0&0\end{bmatrix}
\begin{bmatrix}s_1[k]\\s_2[k]\end{bmatrix} +
\begin{bmatrix}0\\1\end{bmatrix}x[k]\\
y[k] &= \begin{bmatrix}1&0\end{bmatrix}
\begin{bmatrix}s_1[k]\\s_2[k]\end{bmatrix}
+[1]x[k]
\end{align}
$$

Now this is where I think the problematic part is. This system, though being LTI, is not regular (relevant buzz words are regular, impulse-free, index-1 descriptor systems). In other words, there exists no $\lambda$ for which the expression $\det(\lambda E-A)$ is nonzero and hence one of the modes is $\frac{\bullet}{0}$ and in fact this is $\frac{0}{0}$ for this example. This means that our system has solution uniqueness issues, unlike the causal LTI system, there is no guarantee on the existence of a unique solution. There is no guarantee of an admissible solution for that matter(buzzword impulsive systems). Hence the LTI reasoning of the other answers won't cut it.

And this what causes the trouble as far as I can tell from Matt L's argument is that he found two nontrivially distinct solutions for the same system and concluded that this cannot be a linear system. But this also assumes the uniqueness and existence of a solution and initial conditions.

It only differs from the regular systems in the way that the uniqueness and existence guarantees of the standard LTI systems cannot be assumed. The models no longer can be assumed to have admissible trajectories for all possible signals.


Why does sampling frequency affect the filter results?

Limited numerical precision. The higher the sample rate, the closer the poles move to the unit circle, the closer to the unit circle, the less stable the filter is.

There are different implementation methods that are better than others: design as poles, and zeros and not as transfer function, use cascaded second order sections, use correct section ordering, higher precision data types, etc.

In most cases very high order filters with poles close to the unit circle are a bad idea. You get very high group delay and phase distortions and massive time domain ringing.


Thank's to MattL's comment I found the 'answer' myself... I was mixing up equations (1) and (2).

As MattL. pointed out, it makes sense to take the square root of equation (1), resulting in something like $$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|\sigma_x\le
\sigma_y\le\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|
\sigma_x\tag{1}$$ if all terms are positive.

If one takes equation (2), applies the integral and the same scaling as in (3) the result is $$\frac{1}{2\pi}\int_{-\pi}^{\pi}S_y(e^{j\theta})d\theta=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta}) \left|H(e^{j\theta})\right|^2 d\theta \tag{2}$$

The left hand side equals the variance of y, so taking the square root gives the standard deviation. However the right hand side is different, as $\left|H(e^{j\theta})\right|$ is a function of frequency.

Hence in general  $\sigma_y\ne |H(e^{j\theta})|\sigma_x$.


BER of QPSK modulation in Rayleigh fading channel model

Chapter 14 of John Proakis' Digital Communications would be helpful, and Appendix C of the same might actually do that derivation. At the least, Proakis gives the expression for M-ary PSK in the Rayleigh channel which you could reduce for 4-PSK.

I remember it by remembering the basic AWGN channel. Since the Q function is the tail of a Gaussian distribution, BER of QPSK in AWGN is found using SNR and the Q function. Rayleigh distribution can be thought of as a complex valued AWGN. One of the easiest ways to derive Rayleigh is to start with two AWGN PDFs that are orthogonal. So I'm pretty sure you just use the SNR and the CDF called Rayleigh distribution.


Expressing 2D Gabor transform as 1D spectrum

Equation 11.1 shows the equation for $A$:
$$A(\tau,\omega) = A_0 \exp\left(-\frac{\omega \tau}{2Q}\right)$$

We see that $A$ really only depends on the product of $\tau$ and $\omega$, therefore by defining a new variable $\chi = \omega \tau$ we can rewrite this expression as

$$A(\chi) = A_0 \exp\left(-\frac{\chi}{2Q}\right)$$

Perhaps this helps.


Can we calculate moments (mean, var, skewness, kurtosis) of a signal in frequency domain?

If your time series $x(t)=s(t)+n(t)$ where $s(t)$ is a deterministic signal,and $n(t)$ is independent Gaussian noise,  it has a pdf of the form
$$
p(x(t))=\frac{1}{\sqrt{2 \pi} \sigma} \exp\left( -\frac{1}{2} \left(\frac{x(t) - s(t)}{\sigma}\right)^2\right)
$$
the first moment is
$$
E\left\{ x(t) \right\} = s(t)
$$
which is NOT the DC component of the Fourier Transform of a sample. The first moment is not a constant and the time series is not WSS, and the first moment is not equal to the time average.


Downsampling a PSD that was generated by zero-padding - worth it?

Consider these two cases where you operate on a signal $y[k], K = 0, \ldots K-1$:

You take the DFT of the signal $y[k]$ to obtain a length-K frequency-domain signal.
You first zero-pad $y[k]$ to obtain a longer signal $y_p[k], k = 0, \ldots, 10K-1$. You then take its DFT. Finally, you downsample the DFT by discarding all but the first of every 10 samples (i.e. you keep the first, tenth, twentieth, and so on).

The two processes will yield the exact1 same result.

1 In practice, that is, when using finite-precision arithmetic, it's likely that you won't arrive at results that are identical bit for bit. However, in theory, and with the use of infinite arithmetic precision, the two approaches are equivalent.


FFT requirements of real time signal from sensors

Q1: is this a correct setup?

I'm guessing "no" but that really depends on your application and goal. The process will produce numbers but if these numbers are useful or not depends on what you want to do with them.

Q2: is there a minimum requirement for the buffer (FFT window?).

Mathematically speaking : no. You can run a FFT on a single sample if you want to. t doesn't do anything useful but it can be done. The real requirements come from your application.

I have now set it up to process 10 samples. Q3. what is the maximum number of frequencies that can be detected in this scenario?

An FFT setup like this will you give  you ten numbers which are the complex amplitudes from -0.5Hz to +0.5Hz in 0.1 Hz steps. That's not the only frequencies you "detect" but that's a more complicated topic. If you feed a 0.15Hz sine wave in the process all 10 numbers will be non-zero.
You may also get significant aliasing with this setup, if you don't use a decent anti-aliasing filter before the current sensor.


Hint:
Write the relationship between $z(t)$ and $x(t)$ and then simplify it a bit.


Estimating pole radius of formant filter from bandwidth

I have here a very pretty answer supplied by @Jazzmaniac, which I will transcribe from our conversation (on IRC Freenode ##dsp).

The key insight is that the gain is inversely proportional to the distance from the pole.

So, double the distance, halve the gain.

(This can be seen by imagining a single pole on the origin with the transfer function $H(z) = 1/z$)

The small circle represents the peak gain (at frequency $e^{jθ}$).

The bigger circle represents a locus of half that gain.

Hence, it is clear that Φ = bandwidth/2


use chirp z transform for spectral resolution

If you have no prior knowledge about the approximate locations of the frequencies, the Chirp Z-transform is of no immediate use to you. The Chirp Z-transform functions like a magnifying glass, so you need to know where you want to look and the Chirp Z-transform will show you the details. I would suggest you use an FFT to get an idea where the frequencies are, and if you need a very high resolution in a certain area of the spectrum, then the Chirp Z-transform can be useful. However, if computational complexity and memory are no issue, you can simply use a long FFT to achieve the same. If you want to know more about the Chirp Z-transform I would recommend this paper to you.


Quantization and Sampling - putting it all together

I believe that $|\Delta_l|$ in that equation is a typo and should actually be $|\Lambda_l|$, because from the first equation in $(2)$ we have
$$|\Lambda_l|=\frac{\psi_H-\psi_l}{k}$$


Firstly, you do not need to specify an input $x(n)$ to know whether the system is stable or causal, these are properties of the system itslef, I will interpret $h(n)$ you mention as a array describing the system response.

For causality in time domain perspectivefor an LTI system check if $h(n) =0 $ for all $n < 0$, that is the $h(n)$ is essentially a right sided sequence.

For causality in Z domain perspective equivalently, the Region of Convergence is the right side of the outermost pole (outermost in terms of the highest magnitude pole in z domain)

For stability in time domain perspective, the system response $h(n)$ should be absolutely summable, since you have only finite coefficients, this system is stable, any FIR filter is always stable.

For stability from Z domain perspective, the region of Convergence must contain the unit circle.

For stable and causal system all these conditions hold at the same time.


Given the phase coherence over non-blanked intervals, a simple and effective approach would be to additionally window the non-blank intervals with a matching window corresponding to the length of the time samples in between blanked intervals. Since the window only tapers the amplitude and does not modify the phase, there will be no issue with different window durations. This will significantly decrease if not eliminate the grating while provide phase coherent processing gain in the DFT from interval to interval over all usable samples. The Kaiser window would be a good window choice for this application. Consider how the DFT for any given frequency bin, as a correlation, de-rotates a signal that may be at that particular frequency such that each sample is aligned in phase, and then the result is summed; our best strategy would be to maximize the information from the samples that are there alone, rather than make our own attempts as to what we think is there to fill in the blanks, to then further convince our assumptions (unless we had other additional external knowledge, in which case such an approach could make sense).
The windowing approach would result with multiple various sized windows over the complete data set corresponding to intervals of non-zero data, and then an FFT would be performed on the entire data set.
To optimize this process further, Sparse FFT algorithms could be employed to take advantage of the gaps in time, but the windowing approach to modify the non-zero time domain data first would be functionally similar.


Iterative Kalman filters and system parameters estimation

Regarding the miss modeling, I am not sure how you can handle that because a good mathematical model is important inorder to get kalman filter getting to work. usually the more complex your model is, the better results you can obtain.

Regarding the process and measurement noise, that noise will not effect much of your result but even then if you want to estimate that noise covariance matrix, you can use iterative algorithms like LMS or RLS to fine tune your Noise parameters for some observations and use it further.


Read about tinnitus masking on Wiki. Yes, you can mask it with various kinds of noise and/or natural sounds that are "noisy" (ocean, rainfall, etc)

For the synthesis of masking noise, there is plenty of information via google & StackOverflow. Keep in mind that broadband noise (versus bandlimited noise) is present across the spectrum and thus, you only need to be approximately in the same bark frequency band of the tone you want to mask.

Here is a plot showing one measurement of Tinnitus spectrum (from http://www.thelancet.com/cms/attachment/2021743353/2041580240/gr1.jpg). I wasn't able to access the full thelancet article:


Determining Power Spectral Density from Autocorrelation Function and Time-Series Data

Taking the FFT of a time-series does not give you the power spectral density (PSD). What you can do is take the squared magnitude of the FFT to get an estimate of the PSD. And the latter is equivalent to taking the Fourier transform of the (deterministic) autocorrelation function of the data. See also this answer to a related question.

So the inverse Fourier transform of the PSD is the autocorrelation function.


Synchronize a software waveform generator

It is not clear to me how you are determining your phase error so there may be issues there, and it also looks like  you are adding the phase correction to the accumulator directly. Consider adding the corrections to the frequency control word instead which would be analogous to a traditional analog PLL such as that shown in the graphic below (with the NCO being the equivalent to a VCO); this would allow the NCO to be an additional integrator in the loop which is necessary to implement a desired second order type 2 PI loop with your use of a PI Loop filter.  If you feed the corrections to the phase accumulator directly, you are treating the NCO as a phase control device instead of a frequency control device.

A type 2 implementation will track a frequency offset with zero error, while a type 1 implementation (which is what would result if the control was fed to the phase accumulator directly) would have an offset that is proportional to the inverse of the loop gain.

For the phase detector implementation, consider using a simple multiplier followed by a low pass filter to remove the double frequency output. (This multiplier can be analog or digital-- in your case since you are dealing with an analog signal an analog multiplier may be simplest. If you can be ensured the duty cycle is 50%, an XOR gate is perfect for this purpose. Otherwise a passive mixer is often a good choice).   The average value will be proportional the phase error with a quadrature offset (so the Loop would therefore lock in quadrature to the input signal).  Care should be taken that the low pass filter simply removes the double frequency component and is not too tight in bandwidth so as to compete with the loop bandwidth (it is preferred to set the loop bandwidth with the loop filter alone). Depending on implementation (see diagrams below) a "Phase/Frequency Detector" may be a preferred choice for the phase detector given that it won't lock to alias frequency offsets. If a multiplier type phase detector is used, the VCO/NCO output must be constrained to a narrow frequency range to avoid such false locks.

Three variants of a PLL approach are shown below. Depending on your need for digital and analog signal types in the rest of your system and programmability of loop dynamics one may be preferred over the other.

For the last implementation shown (which I think would be simplest) there are plenty of integrated chip solutions out there to make this even simpler. For example this chip: https://www.analog.com/en/products/hmc984.html#product-overview incorporates a phase/frequency detector together with a high range divider. The divider is 14 bits but 18 bits are required to get to 240,000, so the approach in this case would be to "pre-scale" the 12 MHz VCXO output with a divide by 16 to create a 750KHz clock, which would then be divided by 15,000 with this chip. There are many other integrated solutions to consider if you went down this path.


The plot as shown makes sense as nothing has limited the time axis to be going from $-\infty$ to +$\infty$. So we are showing the finite duration sequence as it would appear on such a time axis.

We can limit the time axis to a finite number of samples as used in the Discrete Fourier Transform (DFT) or the an unlimited time axis as used in the Discrete Time Fourier Transform (DTFT). Neither approach would change the number of samples in a finite duration signal.

So the plot and text is implying the time axis extends from  $-\infty$ to +$\infty$. The number of samples in the finite duration sequence would be the non-zero samples.

My confusion when first answering this question was the case of the DFT and how adding zeros is significant to the result, but really this is just approximating the DTFT and not changing the duration of the sequence itself.


How can I specify the distance between values in such row

Just a little bit of python:

import numpy as np

def splitVector(v):
    thres = (np.max(v)-np.min(v))/2
    return v[v<=thres], v[v>thres]

This will yield two vectors with just the lower or higher elements. If you want to preserve overall size and position of arguments, just a little modification will do it:

import numpy as np

def splitVector(v):
    thres = (np.max(v)-np.min(v))/2
    return np.where(v<=thres,v,0), np.where(v>thres,v,0)

This will fill in zeros for all elements not fulfilling the condition.


Discrete-time sampling of filtered white noise

Hmm. For standard deviation, I see https://dsp.stackexchange.com/a/8632/829 which states for uniform power spectral density $N_0/2$, the standard deviation is

$$\sigma^2  = \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,\mathrm df$$

which in my case would yield

$$\begin{align}
\sigma^2  &= \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,df \\
          &= \int_{-\infty}^\infty \frac{N_0}{2}|H(\omega/2\pi)|^2\,d(\omega/2\pi) \\
          &= \frac{1}{2\pi}\int_{-\infty}^\infty \frac{N_0}{2}|H(\omega)|^2\,d\omega \\
          &= \frac{1}{2\pi}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{(\omega\tau)^2+1}\,d\omega \\
          &= \frac{1}{2\pi\tau}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{(\omega\tau)^2+1}\,d(\omega\tau) \\
          &= \frac{1}{2\pi\tau}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{u^2+1}\, du \\
          &= \frac{N_0}{4\pi\tau}\left[\tan^{-1} u\right]_{-\infty}^{\infty} \\
          &= \frac{N_0}{4\pi\tau}\left[\pi/2 - (-\pi/2)\right] \\
          &= \frac{N_0}{4\tau},
\end{align}
$$

and with unit WGN $N_0/2 = 1$ so $\sigma = \sqrt{1/2\tau}$.

Not sure how to compute the PSD of a sampled random process, however.


Downsampling, shifting, high pass and low pass filter commutativity

I am suspicious due to the claim that $\mathcal{R}_0$ is the inverse of $\mathcal{D}_0$. Decimation is not invertible. Once samples/entries are deleted, those values are forgotten. They cannot be reconstructed from what remains.

Let's say the input finite sequences are in $\mathbb{R}^{2N}$, numbered as $x = (x[0],x[1],\ldots,x[2N-1])$. Then $\mathcal{D}_0:\mathbb{R}^{2N}\to\mathbb{R}^N$, and $\mathcal{S}$ is overloaded, so that it maps from $\mathbb{R}^{2N}$ to $\mathbb{R}^{2N}$ and from $\mathbb{R}^{N}$ to $\mathbb{R}^{N}$:
\begin{eqnarray}
(\mathcal{S}x)[n] &=& x[n+1~\textrm{mod}~2N]~~\textrm{for}~~x\in\mathbb{R}^{2N},\\
(\mathcal{S}y)[n] &=& y[n+1~\textrm{mod}~N]~~\textrm{for}~~y\in\mathbb{R}^{N}.
\end{eqnarray}

Let $u = \mathcal{D}_0x\in\mathbb{R}^{N}$:

\begin{equation}
u[n] ~=~ (\mathcal{D}_0x)[n] ~=~ x[2n],
\end{equation}
where $n$ runs from 0 to $N-1$.

Then $\mathcal{S}u = \mathcal{S}\mathcal{D}_0x\in\mathbb{R}^{N}$:

\begin{equation}
\begin{split}
(\mathcal{S}u)[n] &=~ u[n+1 ~\textrm{mod}~N]\\
 &=~ x[2\times((n+1)~\textrm{mod}~N)]\\
 &=~ x[2n+2~\textrm{mod}~2N]
\end{split}
\end{equation}

Now let $v = \mathcal{S}^2x\in\mathbb{R}^{2N}$:

\begin{equation}
v[n] ~=~ (\mathcal{S}^2x)[n] ~=~ x[n+2~\textrm{mod}~2N]
\end{equation}

Then $\mathcal{D}_0v = \mathcal{D}_0\mathcal{S}^2x\in\mathbb{R}^{N}$:

\begin{equation}
\begin{split}
(\mathcal{D}_0v)[n] &=~ v[2n]\\
 &=~ x[(2n)+2~\textrm{mod}~2N]\\
&=~ x[2n+2~\textrm{mod}~2N]
\end{split}
\end{equation}

Hence
\begin{equation}
(\mathcal{D}_0\mathcal{S}^2x)[n] ~=~ x[2n+2~\textrm{mod}~2N] ~=~ (\mathcal{S}\mathcal{D}_0x)[n]
\end{equation}
for $0\leq n < N$ (because all results are in $\mathbb{R}^{N}$), so $\mathcal{D}_0\mathcal{S}^2$ and $\mathcal{S}\mathcal{D}_0$ are the same.


The advantage to using FIR-filters here is that you can (relatively easily) account for the delay.
Assuming filters with constant group delay (see examples here,) - which seems likely given that your filters have an odd number of taps - all you have to do is to remove (n-1)/2 samples (where n is the filter order) from the beginning of your data after it has gone through the filter.

For example:
Filter your data with the Alpha filter.
It is now longer by 133 samples.
Remove (133-1)/2 samples from the start of the filtered data, and the delay is gone.
Do the same for all of your filters, and you can keep the output of all of your filters in sync.

With an IIR, the delay is not that simple to account for.  Getting the outputs of you filter back in sync would be challenging if you used IIR filters.

If keeping the filtered signals in sync isn't important, you could switch to using IIR filters for all of your filters.

This would greatly reduce the computational load as well as the amount of memory needed.

As with the FIR, there are design methods for IIR that use the desired passband, ripple, and attenuation to calculate the needed coefficients so that you should be able to generate coefficients for the IIR filters that match your FIR filters in everything except having a constant delay.

As for the coefficients, I wouldn't worry about the actual values unless you have trouble representing them on your selected hardware.  Use what your design tool gives you for the coefficients.

There are design tools which can help improve your FIR filters.  Check some of the online tools (like this one.)  Given the separation of the various bands (wikipedia entry for EEG,)  it looks to me like your current filters can't be provoding very good performance, any way.  I get only 22dB of attenuation using 133 taps for the alpha filter.

IIR filters that sharp may be difficult to attain, and may be unstable.

It all depends on how much attenuation you need.  Since the bands are defined with only a 1Hz transition, any high degree of attenuation will require high order IIR filters or very long FIR filters.

As with any design project, it comes down to making the best compromise between performance and cost - and you need to know what performance is required.


Follows from the DFT definition. It's defined as

\begin{equation}
X(k) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi \frac{kn}{N}}
\end{equation}

So $X(0)$ is

\begin{equation}
X(0) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi \frac{0 \cdot n}{N}}
\end{equation}

Having $k=0$ gives $e^0=1$ all the time so that

\begin{equation}
X(0) = \sum_{n=0}^{N-1} x(n) 1
\end{equation}

Comparing this to the average

\begin{equation}
\overline{x} = \frac{1}{N} \sum_{n=0}^{N-1} x(n)
\end{equation}

shows that $X(0) = N \overline{x}$


For any LTI (Linear Time Invariant) system (considering the steady state response i.e. ignoring the transient response which generally becomes negligible in a very short period), if the input is a sinusoidal signal, the output is always a sinusoidal signal of the same frequency.  But the amplitude of the output sinusoidal can be different from the input amplitude.  And the ratio of the output amplitude to input amplitude is a constant for that particular frequency - whatever the input amplitude and phase are. Similarly, the difference between the input and output sinusoidal phases is also constant for the particular frequency - irrespective of the amplitude. The transfer function can be used to give this input-output amplitudes ratio and input-output phase difference as a function of frequency.  Hope this helps.  The frequency dependency is derivable from the input/output differential equation in analog case or difference equation in discrete case.


How to shift FFT spectrum by half a bin?

Well, if you're willing to perform multiplications in the time domain, you can cyclically shift the whole spectrum in the frequency domain using:

$$x'[n] = e^{j\pi \frac{f}{F_s/2}n} x[n]$$

Where $x[n]$ is your original input sequence, $F_s$ is your sample rate, and $f$ is the amount of frequency shift you desire.

The DFT of $x'[n]$ will have a spectrum that is cyclically shifted by $f$ Hertz, compared to the DFT of $x[n]$. And there will be no "extra products within the bandwidth of the overall FFT", if by that you meant mixing products, from say multiplying by a cosine to effect a shift.

I suppose you could merge that complex exponential factor into your twiddle factors of a custom FFT implementation, but I would advise against that.  Page 19 of these notes shows how that merge would happen: http://eeweb.poly.edu/iselesni/EL713/zoom/dftprop.pdf if you are dead set on modifying the twiddle factors instead of a separate multiplication.


Like Matt L. pointed out, alpha and beta is unclear as is 'type 1'.

You can use the command fdatool to interactively design filters in MATLAB, although you may need the DSP system toolbox installed.

From there, you can apply the filter using the filter() function.

As an example, if you design the filter using FDAtool, and the filter is an FIR, you can use filter(Hd.numerator,1,X,Fs) if Hd is your filter object, X is the input sequence, and Fs is your sampling rate. This is not the only way to do this, but just the way I would go about it.

X could be generated as a sin(2*pinf*t) where t = 1/fs and n is an index in a for loop or shortened to 1:ns for the number ns you wish to have your input sequence have as a length.


This question is well-addressed in other posts; I would refer you elsewhere, but since you probably want the answer, here it is.

A couple of notes:

There are three knobs you need to play around with to get the right FFT. The first knob is sampling enough (so the 1/700 parameter); in particular, if you don't do it at least twice the maximum frequency, you run into problems of aliasing (by far the best and easiest to understand resource is here: https://blogs.mathworks.com/steve/category/fourier-transforms/). The second knob is that you need to sample long enough (the 50 number below) - you could do 1/1000 sampling but if you do up to 10 - 1/1000 and do the FFT, you will notice the amplitude of the smallest frequency (1/(2 * pi)) is always a little too high.
The linked blog explains VERY well why this third knob is important, which is the zero-padding, which allows you to get a better DFT with the same data.
Finally, yes I am plotting the two-sided FFT, so usually people will take half the amplitude I am using: $$\text{plot}(\text{fshift},\text{abs}(\text{yshift})/\text{length}(x))$$ instead of $$\text{plot}(\text{fshift},2 * \text{abs}(\text{yshift})/\text{length}(x))$$.

But this recreates the right amplitudes. The idea is you need to normalize by the number of data points you are using (I didn't come up with this code myself; I stole it from here: https://www.mathworks.com/matlabcentral/answers/162846-amplitude-of-signal-after-fft-operation).

Without further ado, here's the code you want:

\begin{equation}
t = 0:1/700:50-1/700; \\
x = cos(5.*t)+cos(25.*t)+cos(10.*t)+cos(1.*t)+cos(100.*t)+cos(200.*t)+cos(300.*t); \\
y = \text{fft}(x, 2^{\text{nextpow2}(\text{length}(t))+2)}); \\
n = \text{length}(y);    \\
\text{fshift} = (-n/2:n/2-1)*(700/n); \\
\text{yshift} = \text{fftshift}(y); \\
\text{figure}; \text{plot}(\text{fshift},2 * \text{abs}(\text{yshift})/\text{length}(x)) \\
\end{equation}


Python has a lot of import libraries for image processing, and they are running on platforms that support Linux, such as Raspberry Pi.

Pi runs with ARM based processors, and every limitation of that CPU (especially on SIMD floating point ops) is also a limitation for image processing.

Pi also comes with a dedicated hardware GPU which can be used to accelerate certain image processing (computer graphics) functionality (in addition to its image/video encode and decode functionality)

Pi also comes with built in camera interface that would enable raw processing of camera data, using available Python libraries.


I naively assumed any (N,N/2+1) Fourier image had the right symmetry, but this goes wrong at the 0- and Nyquist frequency because they don't have a symmetric counterpart. Here, you have to impose the symmetry explicitly.

So, if we generate an image:

N = 100
fourier_image = np.random.uniform(-1,1,(N,1+N/2))

we impose the required symmetry as follows:

fourier_image[1:N/2,0] = np.conj(fourier_image[N:N/2:-1,0])
fourier_image[1:N/2,-1] = np.conj(fourier_image[N:N/2:-1,-1])

And doing the same as in the question above, we find the input - and recovered power:

print 'Original power: ', sum(sum(abs(fourier_image)**2.))
print 'Recovered power: ', sum(sum(abs(fourier_recovered)**2.))

which are indeed the same:

Original power: 1664.3177945
Recovered power: 1664.3177945

For an odd-shaped image this will be slightly different but the idea is the same.


If I want to apply a filter to an Image should I use the IIR or FIR filters?

Some important vocabulary first to get the terminology right:

DFT - Discrete Fourier Transform. This is a transform. It is quite computationally intensive. ($O(N^2)$).
FFT - This is an algorithm that makes computing the DFT very fast. ($O(N \log N)$)
Z-transform: This is just a mathematical generalization of the DTFT - it is not an algorithm. When $z = e^{jw}$, then we have the DTFT. Moreover, the DFT is a (uniform) sampling of the DTFT.

So from a practitioner's standpoint, what you will use in your filtering technique is the FFT - at least with a lot of classical techniques. This is because as stated above, the FFT is an algorithm.

In relation to implementing actual difference equation VS doing your filtering via FFT: Most of the time, implementing a convolution is done via the FFT, because of the computational simplicity. There is a point however at which it becomes most efficient to perform a convolution directly, but that is only if your image/signal sizes/lengths are really small.


Discrete signals in spectrum to Analog : DFT

For signal 1 ($x_1(t)$) :
The Fourier's transform ($F(f)$) of $Acos(2\pi ft)$ is $\frac{A}{2}(\delta(f-f_0)+\delta(f+f_0))$ and for $Asin(2\pi ft)$ it's $\frac{jA}{2}(\delta(f-f_0)-\delta(f+f_0))$

Replacing frequencies and amplitudes with the appropriate ones, the Fourier's transform of $x_1(t)$ is :
$$
F(f) = \delta(f-3000 Hz)+\delta(f+3000 Hz)+\frac{j}{2}\left(\delta(f-2000 Hz)-\delta(f+2000 Hz)\right)
$$

Hence the observed Diracs on fig. 1. You did not scale your axis according to the fft sampling while plotting it so the frequency does not appear clearly. The amplitudes are respectively 4 and 8 which is the right ratio so I guess you also did not scale the fft result (matlab doesnt do it itself), it should be 0.5 and 1 according to what I calculated.

For signal 2 it works the same way

edit : there is a mistake in the sampling frequency (1 kHz) you mention. To fit the graphs properly (and respect Shanon) it should be 8 KHz.


Comparing signals recorded using diffrent methods

Different flavors of cross correlation are probably the best choice here. It really depends on how exactly the recordings are different. Two microphones side by side behave very differently than two microphones in different spots in the room. Your best method will depend on a clear understanding of what your application needs to cover and what doesn't.


Ways to identify if music has vocal timbres in it?

Scalable audio separation with light Kernel Additive Modelling

http://www.loria.fr/~aliutkus/kaml/


Understanding the Conditions for Recovering a Discrete Time Signal Through Sampling

Lets have a higher level of the idea of signal reconstruction from samples.

When you try to reconstruct something from partial information it is important to know what you know about the result required.

For instance, if you know the signal is a Line if $ \mathbb{R}^{2} $ then all you need is 2 samples, any samples of it.
You can even generalize that to any polynomial with having enough points samples from it.

You can even find more interesting rules (For instance, can you think about the case of having single Sine Function with Zero Phase ad Unit Amplitude, how many point will you need?).

So when dealing with reconstruction we need assumption of the model and then derive the requirements from the samples.

If we assume the signal has limited bandwidth, namely from a given frequency and above it vanishes in the Fourier domain then the Nyquist Shannon Sampling Theorem states a simple rule, sample the signal in a uniform manner with sampling rate which is larger then twice the bandwidth of the signal.
Once you do that, you have enough information to reconstruct it.

So in your case, saying that the signal vanishes from certain frequency and above is just stating it obeys the assumption of the model.


Ah okay thank you. I think what irritated me was that I didn’t consider that the sampled function is a periodic or periodically continued function so that the Fourier transform delivers amplitude values for the base frequency and multiples of it only (so the Fourier transform is discrete), which are indexed as n*(base frequency), so that it is obvious why F „depends“ on n. The amplitude value for n can then be converted into an amplitude value of a frequency given by n*(base frequency) (and the base frequency can be obtained by the period since a signal is only sampled between its starting point and endpoint, whose distance is the period).


How can I learn image processing in a practical way?

The book you mention in your question is a highly practical one. In fact since digital image processing (DIP) is by definition an applied field, any book on DIP must be and will be practical. They should be showing results of algorithms on images. But if you mean actual code, very little number of DIP books would give away anything other than pseudo or matlab / octave code.

In my opinion the preferred language for research in DIP is Matlab / Octave or Python, and the standard tool for production environment is C/C++ with MSVS (licence) or a number of excelent free alternatives.

Note that pattern recognition and medical image processing are two sub-fields of signal processing and image processing respectively. The former subject requires a lot of statistics and probability theory background, which is not included in a typical DIP book.


What is the use of perception calculation in Speaker Identification?

Well, I don't know whether it'll actually help you – you just said it would!

Now, in any case, using an algorithm to extract features from a signal that mimics or resembles human perception should inherently give those feature vectors a higher "mathematical" resemblance when a human would find the original signals similar, too.

This is obviously what you want in speech recognition – after all, the computer is supposed to do something that is very human.

Now, the question is whether that helps in identifying speakers.

In my first, intuitive, reaction to that, I'd say, no, applying what is used for speech recognition to speaker identification might be a bad idea: A good speech recognizer should be as invariant as possible to the speaker – and thus, "lose" the information of who spoke somewhere on the way.

In reality, something less than that extreme is true: You still want to analyze speech (as opposed, to let's say, interpret engine sounds) – but under a different aspect. So, using the same  basic steps to extract feature vectors sounds pretty clever. You might want to tweak parameters; for example, though not overly important for speech recognition, higher-frequency sound might be helpful in telling high-pitched female teens from smokey-voiced old men.


Pitch detection algorithms "give up" after very little

You might just have the magnitude or gain on your test input signals turned down too low.  Or the threshold set too high.  Or poor microphone placement.  Or perhaps you need an AGC controlled gain block before the pitch estimation algorithms to stabilize the decay portion of any ADSR-"like" amplitude envelopes or sound evolutions.

The above assumes you have a decent S/N ratio (e.g. little background noise and no harmony or accompaniment).

Also, you use a singular "best algorithm", when multiple algorithms (plural) may more appropriate for different portions of a note's envelope, and/or different sized stringed instruments.  In addition, your list of algorithms does not include any of the newer machine learning (ML or Convolutional DNN) inference methods of pitch detection/estimation or octave resolution.  Perhaps even RNNs to use melodic history to improve prediction statistics.


I believe the issue is the OP is simulating a spectrum analyzer with a "Zero-IF" architecture so is seeing the image issues as the sweep signal approach the input signal (since both are real signals). The product is down-converted to baseband with a real LO and passed through a low-pass Gaussian filter, so the resulting signal as it approach the DC center as a very low frequency real sine-wave has increasingly fewer cycles (and ultimately partial cycles) which then contributes to high variability in the power estimate. By using a complex LO instead, the envelope of the resulting signal is constant eliminating that variability (or by observing the two real components of the baseband complex signal $e^{j\omega t} = \cos(\omega t) + j\sin(\omega t)$ we see that as one partial cycle is decreasing the other is increasing to maintain the constant power at the detector).
I confirmed this was the case by changing the LO to a complex signal in the OP's code which would then match the traditional architecture of an actual Zero-IF receiver:
product(:,i) = (signal.*exp(-j*2*pi*freqs(i)*time))';

and use a complex conjugate product for the power (which also could be done using abs(x).^2)
values(i) = 10*log10(sum(product((order+1):end,i).*(conj(product((order+1):end,i)))))-power;

Resulting in:

The other option (that I would NOT recommend for simulation! - unless you need to prove to yourself that it is identical) is to move the filter from being a low-pass filter to a band-pass filter and offset the LO such that an IF frequency passes through the filter (matching traditional spectrum analyzer architectures).
For simulation purposes I would always simulate the equivalent baseband signal using complex signals rather than actual carriers which would be significantly more time consuming to process.
Further, this process can be done in one operation using an FFT which would be significantly faster and for a stationary input such as this produce an identical result. The Gaussian filter could be implemented by multiplying the time domain waveform by the Gaussian window prior to taking the FFT. For non-stationary signals further additions could be made to properly emulate the effect of the sweep time.
As far as the bandwidth, the Gaussian window used in the OP's code with
$\alpha = 10$ has an equivalent noise BW of 5.64 bins as computed from:
$$\text{ENBW} = N\frac{\sum (w[n]^2)}{(\sum w[n])^2} \tag{1} \label{1}$$
That said for a sampling rate of 10e6, and a target resolution BW of 10KHz, this would require a time capture of
$$N = 5.64 \times (10e6/10e3) = 5640 \text{ samples}$$
This would then also match the duration for the Gaussian filter for a 10KHz resolution BW in the OP's code by selecting the coefficients directly from the window function of that length, as in:
b =  gaussian(5640, 10)

I use equivalent noise bandwidth (ENBW) as resolution BW instead of the 3 dB BW, as it can be directly computed from the window as in \ref{1} and accurately indicates the same power that would be in a brick-wall filter of that bandwidth for a white noise source. The ENBW is slightly higher than the 3 dB bandwidth (on the order of 5% but depends on the window).
As recently reiterated by fred harris at the DSP Online Conference (https://www.dsponlineconference.com/), the Gaussian window is NOT the best choice. As he detailed, the reason Gaussian is often selected for spectral analysis is that it is believed to have the minimum time-bandwidth product, so minimum localization in time and frequency. However this is only true for a true-Gaussian which requires an infinite time support and only in that case when $t \rightarrow \pm\infty$ results in the minimum achievable time-bandwidth product BT of $BT = 0.5$. When we truncate the Gaussian response for the time-limited window (as required for practical implementation) the truncated Gaussian $BT>0.5$ and is no longer the minimum BT window function available for the case of finite time duration.
A window that has minimum time BW product under finite time support is the DPSS (Digital Prolate-Spheroid Sequence) window and very close to that is the Kaiser Window, both of which are available in MATLAB, Octave and Python scipy.signal. These windows also have decreasing stop-band attenuation which is also attractive for many applications. For more info on this see: https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html , https://www.dsprelated.com/freebooks/sasp/Kaiser_DPSS_Windows_Compared.html#fig:dpsstest and https://ieeexplore.ieee.org/document/543677


It is difficult to advise about the noisy constellation.
If it is a channel with no added noise, the constellation looks bad.
Not easy to investigate but here are some suggestions to explore:

The signal power might be too low, too close to receiver sensitivity level, so the noise you see on the constellation could be the one of the receiver.
Try to increase the power of the signal to be out of the receiver noise zone.

Automatic gain control issue: Try to freeze the power and AGC or bypass it if possible check the result

Synchronization recovery issue: Synchronize the transceiver and the receiver with their 10 MHz (so you have no carrier and sample offset). Then you can bypass symbol sync module in GNU radio

RCC filter: Maybe your filters do not have the same response in the transceiver and receiver. You could check their roll-off


How do scaling of the coefficients for the poles of a biquad filter affect the gain?

Sorry, that doesn't work. You have to scale all numerator (the "b") or all denominator coefficients (the "a") at the same time. That doesn't work for the denominator since you all practical implementations require $a_0=1$. Scaling $a_1$ and $a_2$ without scaling $a_1$ will move the poles and change the overall shape of the transfer function and not just the gain.


Negative result from filtering a positive signal with a band-pass filter

What you see is what one would expect. As pointed out in Marcus Müller's answer, your band pass filter has a relatively poor stop band attenuation, and since DC is in the filter's stop band, DC is not sufficiently attenuated. You can easily predict what is going to happen:

compute the filter's DC gain, which is just the sum over all filter coefficients: $H(0)=\sum_nh[n]$ This will be a not so small negative number.
Take the DC offset of your input signal ($2^{15}$), and multiply it with $H(0)$. This must be the DC offset of your output signal.

So the DC offset of the output signal can be easily predicted, even without computing the output!

What can you do to get rid the DC offset of the output signal? Since you can perfectly predict it, you can just subtract the DC offset $=2^{15}\cdot H(0)$ from your output signal. Another fun solution is to increase the filter length by $1$ sample, and set that new last sample to the value $-H(0)$ (as computed above). This will make sure that $H(0)$ of the new filter is exactly $0$, and it will only insignificantly change the overall frequency response of the filter.

For the ones who wonder if it is normal that the first and last sample of the impulse response are much larger than all other samples: yes, this is normal, especially for narrow band pass filters, and it is a direct consequence of the equi-ripple optimality criterion implied by the Remez algorithm. These two impulses (first and last sample) are echos generating the sinusoidal (i.e., equi-ripple) stop band behavior in the frequency domain. The ripple in the impulse response generates the (impulse-like) narrow pass band in the frequency domain. So everything can be explained by the basic Fourier relationship "sinusoid $\Longleftrightarrow$ impulse".


Fourier transform of $\sum_{n=-\infty }^\infty e^{-|t-2n|}$

First of all, the signal is periodic, so its spectrum must be discrete. Hence, your final solution ("Original answer") cannot be correct.

I would go the following way:

Calculate $X(f)=\mathcal{F}\left\{\exp(-\lvert t\rvert)\right\}$
Apply the time-shifting property of the Fourier Transform to transform all time-shifts of the function.
See that you can factor out $X(f)$ from the resulting sum and then see that the terms in the sum correspond to a comb function where the diracs are $\frac{1}{2}$ apart from each other.

Edit: Below is some code to plot the function:

Fs = 100  # some sampling frequency for discretizing the signal to plot it
T = 10    # the time range
t = np.arange(-T, T, 1/Fs)  # the time samples

# define the basis function
f = lambda t: np.exp(-abs(t))

# plot the basis function and the periodic extension
plt.plot(t, f(t))
plt.plot(t, sum(f(t-2*n) for n in range(-10, 10)))

Edit: Regarding your comment:

$$
\begin{align}
\mathcal{F}\left\{e^{−\lvert t\rvert}\right\}&=X_1(j\omega)=\frac{2}{1+\omega^2}\\
\mathcal{F}\left\{\sum_{n}e^{-\lvert t-2n\rvert}\right\}&=\sum_nX(j\omega)\exp\left(j2n \omega\right)\\
&=X(j\omega)\sum_n\exp(j2n \omega)
\end{align}
$$

Here's the trick:
$$\sum_n\exp(j2n \omega)=\frac{1}{2}\sum_n\delta\left(\omega-\frac{2\pi}{2}n\right)$$

Have a look at the wikipedia article and also this article for example on why this holds.


What is the basic idea behind Fourier transform?

I think the discrete Fourier transform (DFT) is typically the easiest to use as a base of explanation:

A DFT is a linear transform that operates on a time domain signal. After the transform, the signal is in “the frequency domain”, i.e the signal has now been broken down into a vector/signal which contains information about what sinusoids are contained in the signal, and how strong they are. That’s pretty much it, but I’ll continue on to be more thorough.

So let’s say I have a time domain signal that has a 200Hz size wave and some noise: assuming my sampling rate is sufficient, after I do my DFT, I’ll easy be able to see the sine wave as a peak in the frequency domain, whereas the noise won’t really “add up” to much in the frequency domain.

Of course it’s important to mention that once you do a DFT, there is also an inverse DFT operation that “unapplies” the transform; doing this returns you to the time domain and gives you your original signal back.

For a continuous Fourier transform, the idea is essentially the same, except it is just expanded to calculus rather than discrete mathematics/linear algebra

EDIT:

Okay, we'll try to make this even simpler using some pictures via MATLAB. I've added a good amount of words to explain in detail, but there is almost zero math in the text; hopefully that will be sufficient for you.

We'll start off by showing two signals: one is a sinusoid (with a frequency of 50Hz, that is it oscillates 50 times per second. If you need more information than that, please consult a high school text book), the other is that same sinusoid, but with noise added.

Since you asked what noise was, you can think of noise as random numbers added into the signal which "corrupts" the signal, i.e. it is no longer a pure sinusoid.  A plot of such signals in the  time domain looks like this:

On the top plot, we can see that there is clearly a sinusoid present.  In the bottom plot, it pretty much just looks like a bunch of random numbers.  A DFT essentially breaks a signal down into sinusoidal components. It does this by "looking" for the each sinusoid in the signal (it's a bit more complicated than that, but since you wanted limited to no math we'll skip it).

So now that we know that a DFT can "look" for signals, we take the DFT of our two signals, and plot the results to find the following figure:

In the top plot, we have the DFT of our pure sinusoid.  Recalling that I assigned it a frequency of 50Hz, we would expect there to be a peak at 50Hz, and indeed there is.  You may notice there's a negative frequency component.  There is more mathematics required to explain that which involves unit circles, which you said in your original question that you "don't care about", so we'll just leave  that be.

Now, on the bottom plot, we have the DFT of the sinusoid plus our random noise.  We can see that those two peaks are still there at -50Hz and +50Hz.  When we just looked at the time domain signal, it looked like "nothing" was there, just a bunch of random numbers.  Through using the DFT, we were able to show that there is actually a sinusoid buried in that noise, and we were able to identify it easily in the frequency domain.  This is a very basic motivation for how/why one might want to use a DFT.

Naturally, the inverse DFT (IDFT) does this operation in reverse, i.e. it maps a signal in the frequency domain to the time domain. Knowing that, it's easy to see that if we take the DFT of a signal, perform no other operations, and then eventually take the IDFT of the output, we are given back the  signal we started with.  This is because the DFT is what is called a linear transform.  If you'd like to learn more about that, you'll have to learn more about linear algebra.

This is just one application of the DFT, there is a myriad of other applications, and it is a VERY important subject/tool in signal processing.  The Wikipedia article for the DFT has a list of a wide range of applications. If you're interested in learning more, I'd start there.

If someone else would like to provide details for the continuous time Fourier Transform, please do! I think you'll find its more or less going to be the same, but a bit more mathematical in its explanation.

For completeness, here's the MATLAB code for the figures I've shown:

f = 50;               % sinusoid frequency
fs = 500;             % sample rate
t = linspace(0,1,fs); % create a time vector
% create a sinusoid
x = sin(2*pi*f*t);
% create a noise signal
n = 1.5*randn(size(x));
% add in the noise to the signal x
xn = x+n;
% create a frequency vector that we'll use to label the plot
f = linspace(-fs/2, fs/2, numel(t));

% create a some plots
figure(1);
subplot(2,1,1);
plot(t,x);
xlabel('Time (s)'); ylabel('Magnitude');
title('Pure Sinusoid (No Noise)');

subplot(2,1,2);
plot(t,xn);
xlabel('Time (s)'); ylabel('Magnitude');
title('Sinusoid Buried in Noise');

figure(2);
subplot(2,1,1);
plot(f,20*log10(abs(fftshift(fft(x)))))
xlabel('Frequency'); ylabel('Magnitude (dB)');
xlim([-250,250]); ylim([0 60]);
title('DFT of Pure Sinusoid');

subplot(2,1,2);
plot(f,20*log10(abs(fftshift(fft(xn)))))
xlabel('Frequency'); ylabel('Magnitude (dB)');
xlim([-250,250]); ylim([0 60]);
title('DFT of Sinusoid Buried in Noise');


The order of processing the prerecorded samples (DC removal, Hilbert transform, frequency shifting)

Order of operations can be changed but will change where the processing is applied and that may be simplified in certain conditions which will be very obvious once the actual signal is processed if the goals of each step are understood. For example, if the OP means by DC a non-zero carrier (which will appear at DC after moving the signal completely to baseband), this is much easier to remove after moving to baseband by simply subtracting the mean. However in actual applications, the carrier may not initially be perfectly estimated, in which case a tone close to DC will occur and subtracting the mean will be ineffective (but the tone can be used to continue to refine the carrier estimate).  If the OP meant the waveform itself as received prior to any processing has a large DC offset, this too can be removed prior to any processing with a simple subtraction of the mean (since the samples are pre-recorded this would be simplest) assuming the DC offset is not desired for some reason, otherwise the tone it becomes will be filtered out regardless in the subsequent processing below. (There can of course be advantages to removing the DC signal first in terms of dynamic range available but that would all be clear in the specicic procssing done and precision of the operations.
No Hilbert transform is necessarily needed, if the frequency initially is at a carrier sufficiently higher than its bandwidth, simply frequency translate the real signal to baseband using using $e^{-j\omega_c t}$ which means feed the signal into to multipliers and multiply one with $\cos(\omega_c t)$ and multiply the other with $\sin(\omega_c t)$ to get a complex I and Q output and filter the remaining higher frequency component that will reside at $e^{-j 2\omega_c t}$ (this is a typical Digital Downconverter (DDC) architecture).
For further details of this approach, please see:
https://dsp.stackexchange.com/a/51280/21048


What factors determine synchronous signal transmission ability of a channel?

Not sure what exactly you mean by "synchronous". Most digital transport mechanism (such as Ethernet or 802.11) are "asynchronous", i.e. a continuous data streams  are chopped up in packets, the packets are transmitted whenever the connection allows it, the packets are assembled in the correct order at the receiving end and then turned into a continuous data stream again. Since it's not known when exactly the packet can be transmitted, there has to some buffering at the receiving end so there is some data available even if the next packet comes in late. This results in a certain amount of "latency". In the case of Bluetooth + VoIP you are probably looking at several 100 milliseconds to a few seconds of overall latency.

You may also be concerned about inter-channel sync, i.e. do both channels of the stereo stream stay time aligned. This can easily be achieved by interlacing the channels to create a single stream.

This all depends heavily on the specific implementation: Bluetooth supports many profiles. For your application you probably want the A2DP profile. This uses the SBC codec and there are dozens of ways to configure the codec ranging anywhere from high quality to totally unusable. Both parts of the interface need to agree on the codec implementation.

There are also many different VoiP protocols and codecs out there. Most operate in the 10kb/s-40kb/s range which seems way too low for you application, so you may need to create a custom protocol.


Determining the resulting pixel value before creating a compressed JPEG image

To give you an upfront answer your question you want to determine the $pixel-values$ from $DCT-Coefficients$ probably there is nothing on earth better than doing IDCT itself!

However, leaving aside some vagueness i can re-phrase your question:

If i have an algorithm x - that is defined to work on Pixels, here is what i do:
DCT-Coeffs ->  IDCT  -> pixels  Algorithm_X  -> pixels  DCT  -> DCT-Coeffs
(ignoring other stuff that is peripheral).

  The question is: is there a short cut to this brute force where i am
  obviously undoing quite some stuff done in first stage?

What you are referring to is what is generally known as Compressed Domain Processing.

DCT-Coeffs ->  Compressed_Domain_Algorithm_X  -> DCT-Coeffs

Basically, all compressed domain algorithm exploits the linearity of DCT. For example

$$
Y = \alpha * X_1 + \beta * X_2
$$
then,
$$
DCT(Y) = \alpha * DCT(X_1) + \beta * DCT(X_2)
$$
where $\alpha$ $\beta$ are scalars.

There are many things possible, and while many NOT. It is too vast to answer the set of all compressed domain algorithms which are possible, but i would recommend key readings that will set a great stage for you.

Compressed Domain Processing of JPEG-encoded images by Brian C. Smith , Lawrence A. Rowe Real-Time Imaging (1996)
A Survey of Compressed Domain Processing Techniques by Brian C. Smith
Proc. NSF Workshop on Reconnecting Science
Algorithms for Manipulating Compressed Images, Brian C. Smith, Lawrence A. Rowe,
IEEE Computer Graphics and Applications, Sept. 1993, vol.13, (no.5):34-42.
Manipulation and compositing of MC-DCT compressed video by Shih-Fu Chang, David G. Messerschmitt IEEE Journal on Selected Areas in Communications, 1995/1, Volume
13 Issue 1
Compressed-Domain Video Processing by Susie Wee, Bo Shen, John Apostolopoulo
HP Labs Tech Report 2002.

Specifics:
If you are using libjpeg you can actually get quantized DCT co-efficients directly. However, most above literature refers un-quantized DCT co-efficients.


Applying fourier transform twice (DSP course)

Just use the formula for the geometric sum:
$$\sum_{k=0}^{N-1}e^{-j\frac{2\pi}{N}kl}=\frac{1-e^{-j2\pi l}}{1-e^{-j\frac{2\pi}{N}l}},\qquad l\neq mN,\;m\in\mathbb{Z}$$
Now realize that $e^{-j2\pi l}=1$ for $l\in\mathbb{Z}$ and the result follows.


Algorithm for detecting the time where the signal is above a threshold

There is a book by Basseville and Nikiforov called "Detection of Abrupt Changes : Theory and Application" that they released to the public as a PDF several years ago (it's out of print, now, I believe).

That book looks at the basic CUSUM (cumulative sum) algorithm and how to choose appropriate thresholds for it.


For fractional delay filters see this article: https://www.dsprelated.com/showarticle/22.php.

Timing error delay is not used for carrier recovery. Carrier offset and symbol timing offset are not the same thing.

Carrier recovery can be done at 1 sample per symbol assuming you do timing recovery first. For this a simple phase detector is the cross product phase detector where you derive phase from the imaginary term of the complex conjugate product.

The Gardner Timing Error Detector is a good choice as it can determine timing offset over a relatively large carrier offset. If resolving timing offsets to within 1/10th of a symbol then this would be a very simple approach without further interpolation. However, I question why the OP precludes interpolation; combining a timing error detector with a polyphase interpolator is an excellent approach to resolve timing with fractional delay filters without having to increase the sampling rate- however interpolation techniques are still applied.

I have more information in other posts for both of these approaches listed below:

Gardner TED:

Isn't Gardner's algorithm and Early-Late gate the same thing?

Fractional spaced equalizer + timing (clock) recovery

Phase Detectors and Carrier Recovery Loops

Phase synchronization in BPSK

More Carrier Recovery and Timing Recovery

How to demodulation BPSK in real case(through a channel)?


How to Generate Band Limited Gaussian White Noise in MATLAB?

You would generate bandlimited Gaussian noise by first generating white noise, then filtering it to the bandwidth that you desire. As an example:

% design FIR filter to filter noise to half of Nyquist rate
b = fir1(64, 0.5);
% generate Gaussian (normally-distributed) white noise
n = randn(1e4, 1);
% apply to filter to yield bandlimited noise
nb = filter(b,1,n);


Basics

The amplitude of an IQ signal is just the vector magnitude, $\sqrt{I^2 + Q^2}$.

The power of an IQ signal is the squared magnitude, $I^2 + Q^2$.

When you see a logarithmic (dB) meter, it is usually measuring the log of the power, i.e. $10 \log_{10}(I^2 + Q^2)$. (This can also be calculated as $20 \log_{10}$ of the amplitude, but unless you already have the amplitude that wastes a square root operation.)

Units

Remember, dB is a relative figure. If you just take $10 \log_{10}(I^2 + Q^2)$, then 0 dB corresponds to an amplitude of exactly 1. If your hardware driver takes the usual floating-point convention of the absolute extreme sample values being from −1 to +1, then you can say that your dB power values are dBFS — decibels relative to full scale. Any signals stronger than that level will be clipped, distorting the signal.

dBm, decibels relative to one milliwatt of power, just uses a different reference level. You can convert to dBm just by adding or subtracting the proper calibration value from the dBFS value — but you need to know that calibration for your hardware at the frequency of interest, such as by measuring it (using a signal source of known power output); it is impossible to perform that calibration purely in software since the digital samples are just numbers with no inherent units.

(A mistake I've seen is to refer to the sample values, or parameters which scale according to them such as an amplitude threshold, as being in “volts”; this is complete nonsense unless your ADC (and other hardware) is actually calibrated to one volt. This is unreasonably large for a radio receiver.)

Practical application

If you are just looking to ignore signals that aren't sufficiently strong (this is known as carrier squelch or power squelch), it doesn't matter what units you use, or even if they're logarithmic or linear, because you're just doing a greater-than comparison. The only critical component is that you start with $I^2 + Q^2$ (as opposed to, say, $I + Q$, which would be just plain wrong).

Note on bandwidth you probably don't need to read

It may also be relevant to note that if you filter a signal, you are by definition removing some of the signal power, so the measurement will be smaller.

In particular, a FFT (such as is the primary visual display in tools like SDR#) can loosely be thought of as a large collection of extremely sharp filters; each output “bin” collects some fraction of the input power. Accordingly, the power in each bin depends on the width of the bin. If you divide by the width of the bin in hertz (that value being $\text{sample rate}/\text{FFT length}$) before taking the logarithm, then instead of dB power, you measure dB power spectral density, which has the advantage of being independent of the FFT bin width if the features you care about are wider than one bin (e.g. a wideband modulated signal); if they are narrower (e.g. pure tones) then the power value is more useful.


You want to resample down to alias frequencies higher than your maximum desired formant.
If you understand aliasing in the z-transform and the transfer function of a system then this should be obvious.


Power-law behaviors in frequency can be found in several unrelated observations and systems. This is apparently the case for $1/f$ or flicker noise. Note that an exact $\alpha=1$ exponent might be too stringent, and people can be interested in a wider range, like  $1/f^\alpha$ with $1/2\le\alpha\le3/2$.

The plot above is from 1/f noise.
Since only looking at the decay of a PSD does not provide you with much ideas about the morphology of signals, you might be interested in looking at generative processes, or accumulated actions that may generate such types of spectra, from physics phenomena (superposition of relaxation processes, infinitely large fluctuations, noise in diffusion processes etc.). And then maybe find analogies when the phenomenon you are studying.
Here are three papers of interest:

$1/f$  noise: a pedagogical review by Edoardo Milotti
$1/f$ noise’’ in music: Music from $1/f$ noise by Richard F. Voss and John Clarke

The spectral density of fluctuations in the audio power of many
musical selections and of English speech varies approximately as $1/f$
(f is the frequency) down to a frequency of 5×10−4 Hz. This result
implies that the audio‐power fluctuations are correlated over all
times in the same manner as ’’$1/f$ noise’’ in electronic components.
The frequency fluctuations of music also have a $1/f$ spectral density
at frequencies down to the inverse of the length of the piece of
music. The frequency fluctuations of English speech have a quite
different behavior, with a single characteristic time of about 0.1 s,
the average length of a syllable. The observations on music suggest
that $1/f$ noise is a good choice for stochastic composition.
Compositions in which the frequency and duration of each note were
determined by $1/f$ noise sources sounded pleasing. Those generated by
white-noise sources sounded too random, while those generated by
$1/f^2$ noise sounded too correlated.

Power laws, Pareto distributions and Zipf's law


An $N$th-order linear constant-coefficient difference equation (LCCDE) is of the form

$$y[n]=\sum_{k=0}^{M}b_kx[n-k]-\sum_{k=1}^{N}a_ky[n-k]\tag{1}$$

It is linear because the sequences $x[n]$ and $y[n]$ appear linearly in $(1)$, and it has constant coefficients because the coefficients $a_k$ and $b_k$ do not depend on the index $n$.

LCCDEs are important because they can be used to describe many practically useful discrete-time (or discrete-space) systems, such as linear time-invariant (LTI) filters (which are also called linear shift-invariant (LSI) filters, if $n$ does not represent a time index).

Given appropriate initial conditions, Eq. $(1)$ can be used to recursively compute the sequence $y[n]$ given its past values and the current and past values of the sequence $x[n]$. For causal LTI filters, $y[n]$ is interpreted as the filtered output, and $x[n]$ is the input. In that case, Eq. $(1)$ describes an infinite impulse response (IIR) filter.

If $a_k=0$, $k=1,\ldots,N$, Eq. $(1)$ reduces to

$$y[n]=\sum_{k=0}^{M}b_kx[n-k]\tag{2}$$

which is a convolution sum, and the coefficients $b_k$ represent the system's impulse response. Note that the impulse response has finite length $M+1$, i.e., the system described by $(2)$ is a finite impulse response (FIR) filter.


Laplacian of Gaussian Approximation and Gaussian Blur as the Solution of Heat Equation

I'm not sure I fully understood what's the issue you're having.
Yet I will show a simple property of the Gaussian filter which might make things clearer.

For simplicity, I will use 1D Signal.
Yet it is easy to extend it into 2D.

The Problem

For $ u \in \mathbb{R} $ (1D, unbounded domain), show that a Gaussian convolution with the initial condition solves the linear diffusion equation:

$$ {u}_{t} = {u}_{xx}, \quad u \left( t = 0; x \right) = f \left( x \right) $$

That is:

$$ u \left( t; x \right) = f \left( x \right) \ast {g}_{ \sigma \left( t \right)} \left( x \right) $$

Where $ \ast $ denotes convolution.

Find the function $ F \left( \cdot \right) $ relating the
evolution time $ t $ to the Gaussian standard deviation $ \sigma \: : \: \sigma = F \left( t \right) $.

Solution

Looking at:

$$ \frac{\partial u \left( x, t \right ) }{\partial t } = \frac{\partial^2 u \left( x, t \right ) }{\partial x^2} , \; u \left( x, t = 0 \right) = f \left( x \right) $$

Defining $ U \left( s, t \right) $ as the Fourier Transform of $ u \left( x, t \right) $ with respect to the parameter $ x $.
This leads to:

$$ \begin{align}
\mathfrak{F} \left( \frac{\partial u }{\partial t} \right) & =  \frac{\partial U \left( s, t \right ) }{\partial t}  \\
\mathfrak{F} \left( \frac{\partial^2 u}{\partial x^2} \right) & = -4 {\pi}^{2} {s}^{2} U \left( s, t \right ) \\
\mathfrak{F} \left( f \left( x \right) \right) & =  F \left( s \right)
\end{align} $$

Hence the equation in the Fourier domain is given by:

$$ \frac{\partial U \left( s, t \right ) }{\partial t} = -4 {\pi}^{2} {s}^{2} U \left( s, t \right ) , \; U \left( s, t = 0 \right)  = F \left( s \right) $$

The solution of this Linear Ordinary Differential Equation is given by:

$$ U \left( s, t \right) = F \left( s \right) {e}^{-4 {\pi}^{2} {S}^{2} t } $$

Hence the solution in the time domain is given by (Using the identity $ \mathfrak{F} \left( {e}^{- \alpha {x}^{2}} \right) = \sqrt{\frac{\pi}{\alpha}} {e}^{- \frac{{\left( \pi f \right)}^{2}}{\alpha}} $):

$$ u \left( x, t \right) = \frac{1}{\sqrt{2 \pi 2t}} {e}^{- \frac{{x}^{2}}{4 t}} = \frac{1}{\sqrt{2 \pi {\sigma \left( t \right)}^{2}}} {e}^{- \frac{{x}^{2}}{2 {\sigma \left( t \right)}^{2}}} $$

Comparing it to the classic Gaussian Function yields that $ \sigma \left( t \right) = \sqrt{\left( 2 t \right)} $.

Remark
This is a solution of mine to a problem I was given in a course named Variational Methods in Image Processing.