Conceptually, how does real time convolution reverb work?

It comes down to latency vs. complexity. If your filter is 10 seconds long, you need to store the audio data of the last ten seconds and then you are able to calculate the current output audio sample with a latency of basically zero (ignoring the time required for calculations here) simply by doing:

$$y[0] = \sum_{k=0}^{l} x[-k] \dot h[k]$$

where $l$ is the length of the impulse response $h$ and $x[0]$ is always assumed to be the current sample. While this gives you basically zero lag, it boils down to the performance of naive convolution which can easily be too much especially if the impulse response is very long.

Now you could use techniques such as overlap-save which already perform convolution of very long signals by splitting the signal up and performing the convolution in chunks. However, to be able to calculate the output you first need to capture all the input samples for one chunk so this defines your latency (again ignoring the calculation time). The efficiency of these algorithms is usually in between the direct convolution method and the pure FFT method and the amount of calculations per sample will be less for longer chunks at the cost of latency.

The problem with overlap-save is that the chunk needs to be at least the size of the impulse response so for a 10 second response, the latency will be at least 10 seconds, hardly real-time. For this case so-called partitioned convolution algorithms have been developed that split both, the impulse response and the incoming signal. I've found a very good paper describing these techniques here (Click me) so I will just link it. I can hardly describe it better than they do.

From what I heard, there should be no more than 3-5ms latency or a good musician will notice the delay.


I think I found an answer to this question. So, I will post it here.

For power transfer function of $1-z^{-1}$ we know:
\begin{eqnarray*}
1-z^{-1}&=&1-e^{-jw}\\
&=& 1-\cos(\omega)+j\sin(w)
\end{eqnarray*}
Then by using complex conjugate properties we have:
\begin{eqnarray*}
|1-z^{-1}|^2&=&1-2\cos(\omega)+\cos(\omega)^2+\sin(\omega)^2\\
&=&2(1-\cos(\omega)).
\end{eqnarray*}
And simillarly, for $1+z^{-1}$, we have:
\begin{eqnarray*}
|1+z^{-1}|^2&=& 1+2\cos(\omega)+\cos(\omega)^2+\sin(\omega)^2\\
&=& 2(1+\cos(\omega)).
\end{eqnarray*}

-

Note that the integrator in this question is a special case (see the expression for $S_i$) and is different than the conventional integrator $y[n]=y[n-1]+x[n]$.


A binary input - binary output channel which flips bits with equal probabilities is known as a Binary Symmetric Channel.

(Image Credit: Wikipedia - Binary Symmetric Channel)

The channel capacity $C$ is the maximum (actually, a supremum) mutual information between the channel inputs and the channel outputs over input distributions. You can't get information through a channel more than its capacity. Intuitively, this says that by designing the way you put inputs into the channel well, you can maximize the amount of stuff coming out of the channel that you can recover (and this limit is the capacity).

Mathematically, this is
$$C = \sup_{p(x)} I(X;Y)$$
where $p(x)$ is probability distributions on the input $X$ and $Y$ is the output of the channel.

In this case, you can write the input of the channel as $P(X=0) = q$ and $P(X=1) = 1-q$. With this distribution, you can see that $P(Y=0) = P(X=0 \text{ and no crossover}) + P(X = 1 \text { and crossover}) = q(1-p) + (1-q)p$ and $P(Y=1) = 1-P(Y=0)$. The mutual information between $X$ and $Y$ is given by $I(X;Y)=\sum_{x,y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)}$ where $p(x,y)$ is the joint distribution of $x$ and $y$, and $p(x)$, $p(y)$ are the distributions of $x$ and $y$ respectively. By writing down $I(X;Y)$ and then maximizing over $q$, you will see $C = 1 - H(p)$ where $H(p)$ is the binary entropy function $-p \log_2 p - (1-p) \log_2 (1-p)$. Then, in your case, $p = \frac{1}{2}$ where $H(p) = 1$, so the capacity is $0$.

If you're interested in learning more about this, a few good books to look (at very different levels of difficulty):

Elements of Information Theory (Second Edition), by Tom Cover and Joy Thomas (Easy, Probably Best Starting Point)
Information Theory and Reliable Communication, by Robert Gallegar (A bit harder)
Information Theory: Coding Theorems for Discrete Memoryless Channels by Czisar and Korner (Pretty hard)
Information Theory by Robert B. Ash (Dover book, Easy)

And an unrelated book which I love, but couldn't keep myself from mentioning, is David Mackay's Information Theory, Inference and Learning Algorithms (free on his website), even though it doesn't go over channel capacity.


How/why are the $\mathcal Z$-transform and unit delays related?

Mathematically we can easily show that the z-transform of a unit cycle delay is $z^{-1}$, Just like the s-transform of a time delay $\tau$ is $e^{-\tau s}$. I would like to add some additional insights that may lead to a more intuitive understanding:
Why z??
NOTE: To really understand how the z-transform and unit delays are related (intuitively), it really helps to understand the Laplace domain, what it is used for, and how it maps to the z-domain. For brevity I am including the first graphic covering a main point needed and the following paragraph giving a high level perspective, without going into further details of how we use either plane, but the more you know about that, the more the rest will make sense of what I will show after that.-- For further details on s and z along with a great graphic by Julius Smith that inspired my plots below, see  https://ccrma.stanford.edu/~jos/mdft/Comparing_Analog_Digital_Complex.html)
To help understand $z^{-1}$ it is useful to understand why we use the transform into the s and z plane. This is an entire topic on its own but I will at least make note of this very interesting observation: The process of multiplying (complex conjugate multiplication to be specific) and then integrating the product for two functions in the continuous time domain, or multiplying and then adding the product in the discrete time domain can be viewed as "correlation". (In actuality it is impressive how many things we know in DSP are a form of this description of correlation!).  Take a close look at the equations for the Fourier Series Expansion (FSE), the Fourier Transform (FT), Laplace Transform (LT) and Z-Tranform (ZT) equations and consider how we are doing correlations to each possible frequency for the FSE ($e^{jn\omega_nt}$), each point on the frequency axis for the FT ($e^{j\omega t}$), and each point on the complex surface for S ($e^{s t}$) and Z ($z^n$):
Fourier Series Coefficients
(scaled correlation of x(t) with $e^{jn\omega_0t}$, x(t) is only defined from 0 to T):
$$c_n = \frac{1}{T}\int_0^T x(t)e^{-jn\omega_0t}dt $$
Fourier Transform (correlation of x(t) with $e^{j\omega t}$) :
$$X(\omega) = \int_{t=0}^\infty x(t)e^{-j\omega t}dt$$
Laplace Transform (correlation of x(t) with $e^{st}$):
$$X(s) = \int_{t=0}^\infty x(t)e^{-st}dt$$
Z-Transform (correlation of x[n] with $z^n$):
$$X(z) = \sum_{n=0}^\infty x[n]z^{-n}$$
(x(t) and x[n] assumed to be causal, hence the integrations and summations start at 0).
(Note that for the Laplace and Z Transform it is not strictly a correlation the way I defined it since s and z have real components, so $e^{st}$ is not simply the complex conjugate of $e^{-st}$; both the real and imaginary exponent signs change.)
In each case, we have mapped our time domain function by describing it in terms of other variables, with the primary benefit of revealing its extrema (the poles and zeros) that can uniquely and compactly describe a (linear time invariant) system, opening up a whole world of elegant problem solving (notably solving integro-differential equations with simple algebra, and describing/ manipulating systems based on their pole and zero locations).
So we convert our discrete time domain sequence into the continuous z-domain by correlating it with every possible value of $z^n$ where z is any complex number, thereby decomposing our waveform into an infinite number of z's each with its own complex weight (magnitude and phase). This is the surface of our function on the z domain, which we usually just view the pole and zero locations (as that is sufficient to determine the magnitude and phase of every other point on that surface!).
The unit circle in the z-plane is the frequency axis:
Below demonstrates one way, out of many, to map from the s-plane to the z-plane. For purposes here, the only point I want to make is the reminder that the vertical axis in s which represents frequency $\omega$ (the Fourier transform is the function H(s) when you set $s=j\omega$),  maps to the unit circle $|z|=1$ on the z-plane.
(This mapping is specifically the matched-z transform, other common approaches are the bilinear transform, which I use most often, and the method of impulse invariance.  In all these mappings, what is consistent is the unit circle in z is the frequency axis (imaginary axis in s) and the left half plane in s maps to the inside of the unit circle in z. The answer on this link was very good for more details into the mapping techniques and why: Are there alternatives to the bilinear transform?)
Additionally this answer by Matt L describes additional mapping techniques where the $j\omega$ axis in the s-plane does NOT map to the unit circle in z-plane: z-Transform Methods: Definition vs. Integration Rule

Time Delay of One Sample = $z^{-1}$:
The math to show that the z-transform of a unit cylce delay is $z^{-1}$ is extremely simple (which is why we love the z-transform; we could certainly compute everything in discrete time systems using Laplace if we wanted the punishment).
$X(z) = Z\left\{x(n)\right\} =\sum_{n=0}^\infty x(n)z^{-n} $
$Z\left\{x(n-1)\right\} = X(z)z^{-1}$
So knowing this I provide the following graphics below to make the main point: as we increase powers of z, we are moving forward in time. The first graphic is of a shift register, where we often see each register denoted as $z^{-1}$, indicating that that element provides a 1 sample delay. Therefore the contents of the shift register at any given moment contain the current sample in the left most position, and all the previous samples as we move forward to the right (so we are moving backward in time reviewing the contents of our waveform in the shift register at that moment).  The second graphic shows the same waveform, but as we would be used to displaying it, on an axis that is moving forward in time, and that axis has increasing powers of z.

What is z???
Next, let's talk about z and what it is. For purposes here,the important point to know is that z is any complex number.

Now, knowing that z can be any complex number, and increasing powers of z represents samples as we move forward in time (and with the reminder when you multiply two complex numbers, you add the phases and multiply the magnitudes), let's look at $z^n$ as $n=0,1,2,3...$ for various values of z :

Notice in these plots how when z as any complex number is limited to the unit circle, the waveform that results for $z{^n}$ is a constant amplitude complex frequency, just as we introduced in the first plot how the unit circle represents the frequency axis.  When z=1, the frequency is DC, and when z is magnitude 1 with a very small angle, the frequency is small, and as the angle of z increases, the frequency increases. A delay in time of a complex number with magnitude 1 results in a rotation (the angles add), so if we continuously add delay with time, we will get a continuous rotation (frequency is change in phase vs change in time). This hopefully will give more insight in tying together an interpretive meaning for z, and $z^{-1}$ and how it relates to a unit sample delay.


FIR butterworth also possible or just IIR?

Filters according to those optimality criteria only exist as IIR filters. They are derived from the corresponding analog prototype filters via the bilinear transform, and this naturally results in IIR filters, i.e., filters with zeros and poles (away from the unit circle).

Of course, you can approximate these filters by FIR filters. The most straightforward way would be to truncate the infinite impulse response using some window function. But such an FIR filter is not optimal anymore since it can only approximate the optimal response.


Signal bandwidth of an 802.15.4 signal at 2.4 GHz

I found this link that details the spectrum for 802.15.4: https://www.semanticscholar.org/paper/IEEE-802.15.4-PHY-analysis%3A-Power-spectrum-and-Gupta-Wilson/2892da39a7ecec1945595d46950438a5bb8777af

And pasted the plot below.

This is consistent with what I would expect for a half-sine pulse-shaped pulse for a symbol rate of 1 MHz (and that would apply to BPSK, QPSK, OQPSK etc where the power spectrum for a random data sequence is the Fourier Tranform squared of the base pulse that is transmitted.).  In Offset-QPSK the I and Q channels can be viewed as two independent (in data content) BPSK transmissions, in quadrature, offset in time by half a symbol. Thus each produces an independent but identical power spectrum that adds in power.

Below I compare the Fourier Transforms for a rectangular pulse of 1 us duration vs the half-sine windowed pulse of the same duration. I normalized the power of the two to be the same. Although the pulse shape increases the width of the main lobe, there is a clear spectral improvement by doing the half-sine pulse shape as demonstrated in this plot (Notably, the power spectrum of the rectangular pulse rolls of at a rate of $1/f^2$, or -20dB/decade, while the power spectrum of the half-sine pulse rolls of at a rate of $1/f^4$ in spectrum).

Fourier Transform of rectangular pulse (Red) vs half-sine pulse (Blue):

This clearly is not as good as that which can be achieved with more spectral efficient pulse shaping such as raised-cosine shapes but offers the significant advantage to low power systems of being a constant envelope modulation when used in an O-QPSK modulation. In this case we are able to drive the power amplifier (PA, which typically dominates power dissipation in the transceiver) into saturation without introducing AM/PM distortion. In saturation the PA can operate at a much higher efficiency (Tx power radiated at a given DC power in).


Here's Octave code for a (digital) RRC filter:
pkg load signal;  % Octave needs this; MatLab doesn't
Fs = 16000;  % sample rate
Rs = 400;    % symbol rate
sps = Fs/Rs; % samples per symbol

%
% Root raised cosine pulse filter
% https://www.michael-joost.de/rrcfilter.pdf
%
r = 0.22; % bandwidth factor
ntaps = 8 * sps + 1; % Pulse duration is 8 symbols

st = [-floor(ntaps/2):floor(ntaps/2)] / sps; % symbol time = t/Ts values
hpulse = 1/sqrt(sps) * (sin ((1-r)*pi*st) + 4*r*st.*cos((1+r)*pi*st)) ./ (pi*st.*(1-(4*r*st).^2));

% fix the removable singularities
hpulse(ceil(ntaps/2)) = 1/sqrt(sps) * (1 - r + 4*r/pi); % t = 0 singulatiry
sing_idx = find(abs(1-(4*r*st).^2) < 0.000001);
for k = [1:length(sing_idx)]
    hpulse(sing_idx) = 1/sqrt(sps) * r/sqrt(2) * ((1+2/pi)*sin(pi/(4*r))+(1-2/pi)*cos(pi/(4*r)));
endfor

% normalize to 0 dB gain
hpulse = hpulse / sum(hpulse);

For practical reasons, the filter impulse response is truncated to 8 symbols in duration, as there isn't much benefit in making it any longer.
I'm not sure what you mean by "analog" RRC filter in Matlab.  The continuous frequency response of the RRC filter is well understood.  It results in an infinitely long time domain impulse response, that is the equation in your question with a scaling factor out front.
In the above code, you can make the sample rate as high as you like, to get more closely spaced samples, to have the RRC filter taps get an ever closer approximation to the continuous function in MatLab.
Matlab itself may have facilities for describing continuous functions in the system toolbox, but I am not familiar with them.
As a technical reality, the infinitely long, continuous RRC filter is not realizable in real-life analog components, as it extends infinitely into the past.  Could you implement a truncated continuous RRC filter with analog components?  Well maybe, but why?  Its use is exclusive to digital modulations and it is easily implemented digitally.


How to do addition of sound pressure RMS in a DFT

The power spectrum of the DFT is the square of the magnitudes and is related to the power spectrum of the discrete time signal by Parsevel's theorem:

$$E_x=\sum_{k=0}^{N-1}\lvert x[k]\rvert ^2=\frac{1}{N}\sum_{k}^{N-1}\lvert X[k]\rvert^2.$$

The RMS is just the square root of the power. (I assume you've already taken care of the 1/N since you mentioned your DFT is already in the correct units).

$$\textrm{RMS}_{x}=\sqrt{\sum_{k=0}^{N-1}\lvert X[k]\rvert^2}.$$

You can compute the RMS of a range of bins $[i,j]$ by constraining $k$.

$$\textrm{RMS}_{x_{ij}}=\sqrt{\sum_{k=i}^{j}\lvert X[k]\rvert^2}.$$


Max is right that the difference between Hamming and Hann windows are small. (BTW, I am a proponent of the movement to totally do away with them term "Hanning".  There is no Dr. Hanning nor Mr./Ms. Hanning that the window is named after.)

The Hamming window is 92% Hann window and 8% rectangular window.  Hamming found out that he was able to reduce the height of the maximum side lobe by doing that.

Another issue that I just thought of:  The Hann window (and some others) is a complementary window; that is that the latter half of the window of one frame adds to the first half of the following frame to 1.  For analysis, the property of being complementary is usually not important, but for synthesis (or reconstruction), this property is salient.

So if you like symmetry about 0, the Hann window having even length $N$, is

$$ w\left( \tfrac{n}{N} \right) = \begin{cases}
\tfrac12 + \tfrac12 \cos\left( 2 \pi \tfrac{n}{N} \right) \qquad & |n| \le \frac{N}{2} \\
0 \qquad & |n| \ge \frac{N}{2} \\
\end{cases} $$

The property of being complementary is:

$$ \sum\limits_{m=-\infty}^{+\infty} w\left( \tfrac{n-mH}{N} \right) = 1 $$

then

$$\begin{align}
x[n] &= x[n] \times 1 \\
     &= x[n] \times \sum\limits_{m=-\infty}^{+\infty} w\left( \tfrac{n-mH}{N} \right) \\
     &= \sum\limits_{m=-\infty}^{+\infty} x[n] w\left( \tfrac{n-mH}{N} \right) \\
     &= \sum\limits_{m=-\infty}^{+\infty} x_m[n-mH] \\
\end{align}$$

where

$$\begin{align}
x_m[n-mH] & \triangleq x[n] w\left( \tfrac{n-mH}{N} \right) \\
x_m[n] & = x[n+mH] w\left( \tfrac{n}{N} \right) \\
\end{align}$$

So $x_m[n]$ is a finite-length frame of signal of length $N$ and represents the audio in the neighborhood of $x[n+mH]$ or $m$ hops from the beginning.  $H = \frac{N}{2}$ is the frame hop displacement.

The Short-Time Fourier Transform (STFT) is:
$$\begin{align}
 X_m[k] &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} x_m[n] e^{-j 2 \pi nk/N}
\\
\\
   &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} x[n+mH] w\left( \tfrac{n}{L} \right) e^{-j 2 \pi nk/N} \\
\\
&= \sum\limits_{n=0}^{N-1} \hat{x}_m[n] e^{-j 2 \pi nk/N} \\\end{align}$$

where

$$ \hat{x}_m[n] = \begin{cases}
x_m[n] \qquad  & 0 \le n < \tfrac{N}{2} \\
\\
x_m[n-N] \qquad & \tfrac{N}{2} \le n < N \\
\end{cases}$$


Why make the phase of an FIR filter linear?

Because linear phase of a filter results in constant group delay according to the relation between group delay $\tau$ and phase response $\phi(\omega)$ of the filter:
$$ \tau = - \frac{ d\phi(\omega) }{ d\omega} $$

But then you will ask: why then constant group delay? And the short answer will be that constant group delay will preserve the relative positions (in time or space) of different frequency (narrow band) components of the signal to each other.

Many applications, such as image processing, may require such a fidelity on the processed signals. But not all, and to some degree, non-constant group delay such as provided by IIR filters are well tolerated in applications like audio...


I can recommend the following two books (which are free!):

Introduction to Signal Processing by S.J. Orfanidis
Signal Processing for Communications by P. Prandoni and M. Vetterli

Also take a look at the answers to this question.


How do I determine if a discrete signal is periodic or not?

I would do a normalized autocorrelation to determine periodicity.  If it is periodic with period $P$ you should see peaks at every $P$ samples in the result.  A normalized result of "1" implies perfect periodicity, "0" implies no periodicity at all at that period, and values in between imply imperfect periodicity.  Subtract the data sequence's mean from the data sequence before doing the autocorrelation because it will bias the results.

The peaks will tend to diminish the farther away from the center they get simply because of having fewer overlapping samples.  You can mitigate that effect by multiplying the results by the inverse of the percentage of overlapping samples.

$$
U(n) = A(n) * \frac{N}{\lvert N-n\rvert}
$$
where $U(n)$ is the un-biased autocorrelation, $A(n)$ is the normalized autocorrelation, $n$ is the offset, and $N$ is the number of samples in the data sequence that you are checking for periodicity.

EDIT: This is an example of how to tell if the sequences are periodic.  The following is Matlab code.

s1 = [1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 0 0 0 1 0 1];
s1n = s1 - mean(s1);
plot(xcorr(s1n, 'unbiased'))

The "unbiased" parameter to the xcorr function tells it to do the scaling described in my equation above.  The auto-correlation is not normalized, though, which is why the peak at the center is around 0.25 instead of 1.  That doesn't matter, though, as long as we keep in mind that the center peak is perfect correlation.  We see that there are no other corresponding peaks except at the outermost edges.  Those do not matter because there there is only one sample overlapping, so that is not meaningful.

s2 = [1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0];
s2n = s2 - mean(s2);
plot(xcorr(s2n, 'unbiased'))

Here we see that the sequence is periodic because there are multiple unbiased autocorrelation peaks with the same magnitude as the center peak.


Effect of Resizing (Up Scaling) an Image on the Ability to Detect / Recognize Objects

It is theoretically impossible to add information by increasing the resolution of an existing image.

So from that point of view -- no.

However, I certainly find it conceivable that if you're using some canned algorithms they may be sensitive to pixelation.  Unless they're really lame, I would only expect this to be an issue if your original images are very sharp -- if the un-differenced objects tend to have blurry edges, either because of the limitations of your microscope or the properties of the objects themselves, then I'd think about trying what you suggest.

Otherwise -- you can't get something for nothing!


Decimation Question about using a low pass filter before downsampling

Sample rate conversation is easy in theory but tricky in practice.
Assuming you want to convert to the standard rate of 44.1 kHz (not 44 kHz), you have an awkward conversion ratio. $3800 =2^3 \cdot 5^2 \cdot 19$ and $441 = 3^2 \cdot 7^2$  are mutually prime that means that rational sample rate conversion is impractical,so you need irrational sample rate conversion. If you do not need 44.1 kHz exactly you can use $69/8  \rightarrow 44058Hz$ or $112/13 \rightarrow 44107Hz $
In either case you need a low pass filter. The choice of the filter depends A LOT on your specific application. Specifically

How much energy is above 20 kHz. There really shouldn't be a lot but you don't know before you look. If there is a lot above 40 kHz you can perhaps pre-filter this out before doing the actual conversion
What amount of aliasing and signal to noise can you tolerate.
Do you care more about preserving transients in the time domain or an accurate phase spectrum
Up to which frequency do you have to meet your requirements. Designing a converter that's good up to 20 kHz is a LOT harder than designing one that's good up to 18 kHz. Don't even try to go up to 22 kHz.
Do you care about latency, memory consumption and/or execution time.

The "art" of building a good sample rate converter is to really understand your requirements and trade-offs and then optimize the filter design for your specific case.
A bunch of "generic" hints for audio

Polyphase FIR filter is typically a good choice and a good trade off between time smearing and phase preservation
32 taps, 64 phases gets you around 90-100 dB SNR at 18 kHz and much better below.
If you need to do irrational sample rate conversion, do linear interpolation between the phases.
You have a conversation factor that is "far away" from unity. In this case, it's probably best to convert this in multiple steps and not in one go. Maybe from 380->95 (down by 4) and then from 95->44.1


Why Convolution between Signal and a Channel Is Not Done Directly in SIMO Channel

I've said, in another comment, convolution using the function conv "i.e in MATLAB" and convolution using the Toeplitx matrix must give the same results. That's ok.

Now, according to your code, the received signal $r$ is the results of convolution between channel $h$ and emitted signal $x$, which means $r = h*x + n$
* indicate to the convolution (which is circular convolution in case of using $CP$ with OFDM or with any other system).

So, in that conventional known case, we are using SISO system, where 1 antenna is used as transmitter and 1 as receiver, The length of our parameters should be:

$x = N$ x 1 ; $h = L$ x 1;

$N$ is the length of our signal, in your case $Q$

$L$ the length of channel or we sometimes call it IR, in your case $M$

Till here, it's clear, it's the normal process which can be read anywhere. Now suppose that you are using $SIMO$ system with $P$ antennas at reciever instead of 1 (by the way fractional sampling is also equivalent to SIMO system).  in that case you suppose to have $P$ copies for your signal instead of one. It's like you are doing $P$ times convolutions for your emitted signal with $P$ different channels. let's say $P$ = 4; means you have 4 receiver's antennas equivalent in our case into 4 different channels.

As mentioned, you suppose to have 4 copies for your signal, and 4 different channels compared with conventional case, so let's say parameter $H$ is toeplitz matrix which will represent those four channels in SIMO system. So we will have

$R = HX + N$

$R$ is received signal in SIMO, $H$ is the toeplitx matrix, $X$ emitted signal in SIMO system too and also $N$ represent the noise.

Now, your question how to build $H$ and what's dimension of $H$ when using SIMO system. (by they way, your code is correct)

The emitted signal $X$, represents 4 copies of $x$, so its dimension is $(N +L)$ x 1;

The toeplitz matrix $H$ has a dimension $NP$ x $(N+L)$ where $P$ = 4 in your case.  And the noise $N$ should be $NP$ x 1. (N as noise is different about N of emitted signal in this case).

So, Now, you can conclude the dimension of received signal $r$ and $R$ in every case easily.

and regarding your question, why using $x1$, it's to have dimension of emitted signal $X$ of $(N +L)$ x 1 instead of dimension $N$ in conventional case

thnx


Differential Pulse Code Modulation (DPCM)

Eq. $(5)$ in your question is really just a consequence of Eqs $(1)$, $(3)$, and $(4)$, so you have four independent equations. But that's all you need because there's a simple one-to-one correspondence between $\hat{d}[n]$ and $c[n]$ (given by the encoder). Note that the quantization error $e[n]$ is usually modeled as a zero-mean white random process with a uniform PDF.

You generally cannot compute an analytical solution of all the sequences in the coder due to the quantizer, which is of course a non-linear device. However, you can compute the sequences sample by sample by simply simulating the DPCM system as shown in the diagram in your question.

Assuming that the prediction filter has an initial state of zero (i.e., $\tilde{x}[0]=0$), and that the quantization step $\Delta=1$ for simplicity, we can easily compute all sequences given the input sequence. With an input sequence of x=[2,3,5,7] we obtain

$\begin{array}{cccrcr}
n & x[n] & \tilde{x}[n] & d[n] & \hat{d}[n] & \hat{x}[n] \\
0 & 2 & 0.00 & 2.00 & 1 & 1.00 \\
1 & 3 & 0.80 & 2.20 & 1 & 1.80 \\
2 & 5 & 1.44 & 3.56 & 1 & 2.44 \\
3 & 7 & 1.95 & 5.05 & 1 & 2.95
\end{array}$


Amplitude response of Lynn's low pass filter at the 3dB cutoff frequency

Here's your mistake, you forgot that $ \omega = 2\pi f$


Formula to compute output of a linear system

It's correct. Here's the general way how you can solve such a problem – in this specific example it would actually be easier to just solve the integral, but that seems indeed like homework you should be doing yourself.

I'll call the step function $\Theta$, that's how we always call it in physics – hope it doesn't confuse, and the convolution
$$
  (f\star g)(a) := \int_{\{(t,t')\colon t+t'=a\}}\!\!\!\!\!\!\!\mathrm{d}(t,t')\;f(t)\cdot g(t') = \int_\mathbb{R}\!\mathrm{d}\tau\; f(\tau)\cdot g(t-\tau)
$$
which is just the definition that you use. We know the convolution theorem
$$
  \mathrm{FT}(f\star g)(\omega) = \frac1{\sqrt{2\pi}}\mathrm{FT}(f)(\omega)\cdot\mathrm{FT}(g)(\omega)
$$
(the factor $\frac1{\sqrt{2\pi}}$ depends on your definition of the Fourier transform). To apply this to your problem first calculate
$$
  \mathrm{FT}(h)(\omega) = \frac1{\sqrt{2\pi}}\int_\mathbb{R}\!\mathrm{d}t\; e^{-t}\Theta(t)\cdot e^{-i\omega t}
 = \frac1{\sqrt{2\pi}}\int\limits_0^\infty\!\mathrm{d}t\; e^{-t}\cdot e^{-i\omega t}
$$
since the Heaviside function cancels everthing for $t<0$, then this is
$$
 = \frac1{\sqrt{2\pi}}\int\limits_0^\infty\!\mathrm{d}t\; e^{-(1+i\omega) t}
 = \frac1{\sqrt{2\pi}}\left.\left[\frac{e^{-(1+i\omega) t}}{-(1+i\omega)}\right]\!\!\right|_{t=0}^{\infty}
 = \frac1{\sqrt{2\pi}}\Bigl(\frac{0}1 - \frac{e^0}{-(1+i\omega)}\Bigr)
$$
which is just $=1/(\sqrt{2\pi}\cdot(1+i\omega))$. Meanwhile,
$$
  \mathrm{FT}(x)(\omega) = \mathrm{FT}(\lambda t.1)(\omega) + \tfrac12\mathrm{FT}(\lambda t.\cos(400\pi t))(\omega).
$$
These Fourier transforms only exist in a distributional sense; it is known (can easily be shown by inverse $\mathrm{FT}$) that
$$
  \mathrm{FT}(\lambda t.1) = \sqrt{2\pi}\;\delta
$$
and
$$
  \mathrm{FT}(\cos)(\omega) = \frac{\sqrt{2\pi}}{2}\bigl(\delta(\omega-1)+\delta(\omega+1)\bigr).
$$
So we have
$$\begin{aligned}
  \mathrm{FT}(y)(t) =& \frac{\sqrt{2\pi}\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{\sqrt{2\pi}(1+i\omega)}
  \\=& \frac{\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{1+i\omega}
\end{aligned}$$
and therefore
$$\begin{aligned}
  y(t) =& \mathrm{IFT}(\ldots)(t)
   = \frac1{\sqrt{2\pi}}\int_\mathbb{R}\!\mathrm{d}\omega\; e^{i\omega t}\frac{\bigl(\delta(\omega)+\tfrac14\delta(\omega/400\pi-1)+\tfrac14\delta(\omega/400\pi+1)\bigr)}{1+i\omega}
 \\=& \frac{e^0}{\sqrt{2\pi}\cdot 1} + \frac{e^{-400 i\pi t}}{\sqrt{2\pi}(1+\tfrac1{400\pi}i)} + \frac{e^{400 i\pi t}}{\sqrt{2\pi}(1-\tfrac1{400\pi}i)}
\end{aligned}$$
We could know calculate this as an exact expression of trigonometric functions by writing out the $1\pm\tfrac{i}{400\pi}$ as a polar complex number.


This happens frequently if your poles are reasonably close to the unit circle. Consider the following example

%% TF2ZP is problematic
fs = 44100;
% 6th order lowpass, fc = 50Hz, sampled at 44.1kHz
[z,p,k] = cheby2(6,80,50*2/fs);
% to transfer function
[b,a] = zp2tf(z,p,k);
% back to zpk
[z1,p1,k1] = tf2zp(b,a);
display([p p1]);

Displaying the poles side by side, you get

         Correct Poles                  TF2ZPK
      0.99729 + 0.00078062i       1.0007 +          0i
      0.99729 - 0.00078062i      0.99983 +  0.0033549i
      0.99814 +  0.0019985i      0.99983 -  0.0033549i
      0.99814 -  0.0019985i      0.99691 +  0.0030633i
      0.99936 +  0.0025696i      0.99691 -  0.0030633i
      0.99936 -  0.0025696i      0.99535 +          0i

The poles are very different and one pole has actually moved out of the unit circle and the resulting filter is unstable.

This is caused by limited numerical precision. $tf2zp()$ requires calculating the roots of a polynomial which is a numerically tricky problem especially if the roots are very close together. $zp2tf()$ works perfectly fine but  $tf2zp()$ doesn't.

Morale of the story: If you have poles close to the unit circle, keep your filters in ZPK or SOS representation and avoid transfer functions. Avoid $tf2zp()$ if possible.

  Are MATLAB function zp2tf() and tf2zp() are complementary or not?

Sort of. They are intended to be complimentary but like any complimentary pair this is subject to numerical noise of the implementation. In this particular instance there are reasonable use cases, where the numerical noise can be substantial to the point where the function becomes useless.


Difference Between Correlation and Convolution in the Context of Image Processing

Convolution is correlation with the filter rotated 180 degrees. This makes no difference, if the filter is symmetric, like a Gaussian, or a Laplacian.  But it makes a whole lot of difference, when the filter is not symmetric, like a derivative.

The reason we need convolution is that it is associative, while correlation, in general, is not. To see why this is true, remember that convolution is multiplication in the frequency domain, which is obviously associative. On the other hand, correlation in the frequency domain is multiplication by the complex conjugate, which is not associative.

The associativity of convolution is what allows you to "pre-convolve" the filters, so that you only need to convolve the image with a single filter. For example, let's say you have an image $f$, which you need to convolve with $g$ and then with $h$.  $f * g * h = f * (g * h)$. That means you can convolve $g$ and $h$ first into a single filter, and then convolve $f$ with it. This is useful, if you need to convolve many images with $g$ and $h$. You can pre-compute $k = g * h$, and then reuse $k$ multple times.

So if you are doing template matching, i. e. looking for a single template, correlation is sufficient. But if you need to use multiple filters in succession, and you need to perform this operation on multiple images, it makes sense to convolve the multiple filters into a single filter ahead of time.


Digital 16-bit PCM audio signal scaling and STFT spectogram amplitude

1: Scaling

Audio corresponds to a pressure wave, so it is mostly symmetrical in amplitude. Therefore audio is normally represented as a signed value (with 0 average). [-1,1) is the standard representation.

You may offset your signal by 1/2 and represent the wave in the range [0,1), but this adds nothing but a huge component at 0 Hz in your spectrogram.

2: Warining when log(0)

You cannot calculate log(0). If that bothers you, create a function that produces a minimum value. Something like mylog(x) = log10(max(abs(x),1E-5)).

3: Small values in spectrogram

Values in spectrogram indicate energy per frequency bin. The more bins, the smaller values you may get there. And you may well have a 0 as result (see (2) above).

4: Signs again

Stick to signed values. See above.


First, it is not true that an error in the parity portion of a Reed-Solomon code word renders the word useless.  RS codes can correct errors in their parity and their data.

Having said that, 10% errors is very high and I would not use a RS code in that situation.  In fact, in a situation with that many errors you will occasionally have > 50% errors just due to random distribution (random distribution is often "clumpy") of the errors.  No code you use will correct all of the errors all of the time.  Your system should be designed with that fact in mind.

I would use a big turbo code or LDPC code.  Either one will get you near optimal performance.  The code needs to be "big" in the sense of using large blocks to avoid the clumpiness problem noted above.  You will still get more errors on occasion, but (assuming a gaussian distribution) your tails will be MUCH smaller if you make your blocks large.

If your errors are bursty you can use interleaving to spread them out among multiple blocks.  If they aren't bursty, interleaving won't help any.


To explain the supposed contradiction between the number of samples calculated by symbol length and sampling rate on the one hand and the IFFT size on the other hand, you have to take into account that the IFFT output is upsampled by a factor of 4. This means that for every output sample of the IFFT block there are four samples after upsampling. There's a cyclic prefix of 16 samples so that we arrive at a total number of samples of $N=4\cdot(64+16)=320$ as stated on page 53.

On the other hand we can calculate the number of samples by $N=Tf_\mathrm{S}=32\mathrm{ms}\cdot8\mathrm{kHz}=256$. The symbol duration $T$ has been taken from page 108. It think the guard interval has not been taken into account there and thus $N=64\cdot4=256$.

The process that is described as "upconversion" is exactly the frequency shift that you are observing (cf. figure 3-8). Citing from the report:

  Once the baseband signal has been interpolated it is multiplied with
  the audio band carrier frequency to obtain a audioband representation.

To make this more clear using a formula, let $x(t)$ be the complex, analog OFDM signal (after digital-to-analog conversion (DAC)) and let $f_0$ be the desired carrier frequency. Then the radio frequency, upconverted signal $y(t)$ is
$$
y(t) = x(t)\operatorname{exp}(j2\pi f_0 t)
$$
This is a basic concept required in almost every radio system. For more information you might want to do some reserach on radio frequency upconversion and I/Q modulation.

The baseband bandwidth of an OFDM signal is approximately $f_\mathrm{S}/2$. But again, here we have to take the upsampling into account. If it is done correctly, it won't change the bandwidth of the signal. I.e. the sampling frequency before upsampling is actually $f_\mathrm{S}/4 = 2\mathrm{kHz}$ and we have to use this frequency for calcuating the OFDM signal bandwidth wich here is $B=1\mathrm{kHz}$. Finally, the band occupied by the radio frequency signal is $f_0 - B\ldots f_0+B$ which corresponds to 1kHz to 3kHz.

Concerning your suggestion to use a greater IFFT size and to left some subcarriers unmodulated to achieve upconversion - this is a valid approach in this specific example. Note however, that in general it requires a larger sampling frequency which is often undesirable.


there is a relationship between filterbank analysis and the STFT analysis, but it is not an identity.

perhaps it can be identical, where for filterbank, you are slicing the frequency domain into band slices and for STFT you are slicing the time domain into windowed segments (in audio we might call them "grains" or maybe sometimes "wavelets") in time.  but to make the two identical will require some thinking.  the BPF can't be simple regular-old BPFs.


Implementing a bandpass filter in R - this code's logic confuses me

Without knowing what values are being passed in, I will assume the +2 is an adjustment for R using 1 based indices.
The band in the upper half of the frequency data represent complex frequencies symmetric about the fft length/2. The motivation is using a property of the DFT that will result in a real value output from the ifft: http://en.wikipedia.org/wiki/Discrete_Fourier_transform#The_real-input_DFT


Norm MPEG-1 Layer III (Mp3) Filter banks distortion

I have tested your code.

There are problems with the decimation and upsample stages. When you run the analysis and synthesis filters at the undecimated rate, then you get the perfect reconstruction (plus a filter shift)

First of all in case of applications which are sensitive to group delay, you should better use the conv / filter function and then compansate for the group delay.

Also, the decimate function applies a lowpass anti-aliasing filter, which kills the signal content at the $k$-th channel considering the fact that its output is a bandpass signal rather than lowpass. Hence you should avoid the anti-aliased decimation but instead employ a bandpass decimation strategy.
A direct subsampling by 32 would provide a convincing result.

Same goes with interpolation part. After you expand your signal, you should select the proper band of image spectrum by a suitable bandpass filter, rather than a lowpass filter from a conventional interpolator. A direct expansion may convince as its output will already be bandpass filtered later by the polyphase filterbank.

That being said, your main problem is the document that you based your code. That's far too out of signal processing context. I suggest you read the MPEG audio standard book of BOSI. But it's up to you of course.

The following is a working corrected MATLAB / OCTAVEcode:

% S0 - Load the prototype lowpass filter impulse response:
% --------------------------------------------------------
load h2.mat;         % h[n] is the prototype lowpass filter... of length 512
L = length(h);

figure,subplot(2,1,1)
stem([0:L-1],h);title('The Prototype Lowpass Filter');
subplot(2,1,2)
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(h,4*L)))));
grid on;

% S1 - Create the 32 band filter bank by cosine modulation
% --------------------------------------------------------
numbands = 32;                  % number of bands
n=0:L-1;

hha=zeros(numbands,L);         % band of filters
for k=0:1:numbands-1
   hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end

figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
title('ANALYSIS FILTER BANK');
hold on
for k=1:numbands
    plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end

% S2 - Design the synthesis filterbank:
% -------------------------------------
numbands = 32;                     % number of banks
n=0:L-1;

hhs = zeros(numbands,L);           % band of filters
for k=0:1:numbands-1
   hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end

% S3 - Generate a test input signal:
% ----------------------------------
N = 2*1024;      % NOTE: According to the MPEG specification a determined LENGTH must be applied.
wav_in = sin(0.05*pi*[0:N-1]);        % pure sine tone

figure,subplot(2,1,1)
plot(wav_in);title('input signal')
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(wav_in,4*N)))));

% S4 - Apply a wav file to it,  get every output of the every sub-band:
% --------------------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands));   % decimated outputs..
yyu = zeros( numbands , N);

for k=1:1:numbands
    temp = conv( wav_in,hha(k,:));
    yyu(k,:) = temp(L/2+1:L/2+N);
    yyd(k,:) = yyu(k,1:numbands:end);
end

% S6 - Apply synthesis on the decimated signal
% --------------------------------------------
ys = zeros(1, N);
yi = zeros(1,N);

for k=1:1:numbands
    yi(1:numbands:end) = yyd(k,:);
    temp = conv(yi,hhs(k,:));
    ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;

% Display the result:
figure,subplot(2,1,1)
plot(ys);title('Synthesized Back');
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(ys,4*N)))));

With the outputs being:

1- The prototype filter:

2- 32 channel filter bank:

3- A pure tone input:

4- An reconstructed back output:


How to eliminate noise that has the same frequency of my signal

Whenever the desired signals and undesired signals are at the same
frequency band, you can't use the bandpass filter because it will remove both the desired and undesired signals.

In such cases statistical characterization for filter design is needed. Means you need to find out those statistical characteristics of background voice (noise) which differs your required signal from the noise. Examples include variance, expectation / mean or there can be others depending upon background noise, you might need to design adaptive filter if the noise characteristics changes with time.


MATLAB produces too short output for frequency modulation

The function fmmod cannot know the sampling frequency of your input signal. It assumes that the input signal has the sampling frequency Fs that you specify in the function call. Check out the example in the fmmod documentation.

By the way, in a software simulation, there's usually no need to choose such a high carrier frequency.


minimum frequency in periodic components

No, there is no lower bound to the sampling theorem. Any signal $x(t)$ sampled at rate $f_\mathrm s$ can be reconstructed from its sampled version if it contains only spectral components lower than (1) $f_\mathrm{s}/2$. I.e. the frequencies contained in $x(t)$ are allowed to be as low as 0 Hz. Note however, that perfect reconstruction of $x(t)$ is only possible if an infinite number of samples $N$ has been drawn from $x(t)$.

When it comes to spectral analysis of $x(t)$ using its sampled version $x_n = x(n/f_\mathrm s)$ the number of acquired samples $N$ plays an important role because it determines the frequency resolution of the DFT. In order to be still distinguishable the lowest frequency cannot be arbitrarily low. This is intuitively understandable as the lower the frequency the longer the period and thus the signal should be measured for a sufficently long time. Here is some further reading about DFT frequency resolution and spectral leakage

Number of FFT points required for a specific frequency resolution for an oversampled signal

FFT of sine wave not coming as expected i.e single point

(1) Actually $x(t)$ must not contain a sine wave with exact frequency $f_\mathrm s$. For an interesting discussion here on dsp.se see Nyquist Frequency Phase Shift


Different representations of frequency space of 2D image FFT

The second image you find online as you mentioned is magnitude and the first image must be either real part or the imaginary part (but not both). If you calculate complex value magnitude for each point, i.e. $Mag_{i,j} = R_{i,j}^2 + I_{i,j}^2$ you get something similar to the following image.

I = imread('cameraman.tif');
F = fft2(I);
F_Centered = fftshift(F);
Mag = abs(F);
imshow(Mag);

but taking log10 harnesses the large magnititde of some frequency bins so for visulasation the normalization does not round them to zero.

imshow(log10(Mag),[])

Thus applying the filter on the above image which is only processed for visualization is not correct and inverse FFT of the values would not have much meaning.

To apply your filter you get the fft of the image, F(image) , then you need to make a same size mask image (all 1 and the frequency section you need to be removed must be zero, also remember symmetry feature of Fourier) and then apply do a element wise multiplication of fft values with this mask.

Also, I would suggest to watch this short lecture by Hoff on this subject:

https://youtu.be/02c6ohQV2TA


How to get a Rician fading channel matrix from a Rayleigh fading channel matrix?

A Rician fading  channel is just a Rayleigh fading channel with a specular (that is, unfaded) signal component.  Depending on what exactly you mean by the  Rayleigh fading channel matrix $H$, the answer could be as simple as changing the distribution of each element of $H$ from a complex zero-mean proper Gaussian random variable to a complex Gaussian random variable with mean of the form $A_{m,n} \exp(j\theta_{m,n})$.


Channel capacity with additive Gaussian noise with changing variance

I think you are looking for the information rate,

$$ I(x_{1},x_{2},\cdots,x_{n}; y_{1}, y_{2}, \cdots, y_{n}) = \frac{1}{n} \sum_{i=1}^{n} I( x_{i};y_{i} ) $$

when input achieves the capacity of the channel $ x \sim N(0,P) $. Then,

$$ C = \frac{1}{2n} \sum_{i=1}^{n} \log \left( 1 + \text{SNR}_{i} \right) $$

I think every time we use the channel we have a new random variable $ y $ since we have a memoryless channel and $ w $ changes. Then, we have a sequence of related random variables $ \left\lbrace x_{i} \right\rbrace $ and $ \left\lbrace y_{i} \right\rbrace $. With the information rate we compute the mutual information between this two sequences.

A good reference about this Elements of Information Theory.

I hope this could help you.


Python MNE - reading EEG data from array

mne.EpochsArray is for 3-D data (epochs * channels * times). mne.RawArray is for 2-D data. Use EpochsArray.

events is an n * 3 integer array. The 3 columns are: time (in sampling points), length (you can put a dummy here - it is almost never checked - but you still need 3 columns), value (e.g. condition). You fed it a 1 * 2 array.

Try it:

import numpy as np
from mne import create_info, EpochsArray
n_epochs = 100
channels = ["a", "b", "c", "d"]
n_channels = len(channels)
events = np.array([np.arange(n_epochs), np.ones(n_epochs), np.ones(n_epochs)]).T.astype(int)
d = np.random.random((n_epochs, n_channels, 256))

info = create_info(channels, 256, "eeg")
epochs = EpochsArray(d, info, events)

epochs.plot(show=True, block=True)


historical context of impulse response-convolution

Stanley: Thanks for your input. Your answer led me to the "brilliant" idea of googling myself. Based on the link below, I think one can define the birth of convolution as sometime between 1760-1900 depending on how one defines birth. https://pulse.embs.org/january-2015/history-convolution-operation/.

Either way, it seems that Koyck was unaware of it since  he came  up with his model in the 1950's. He would have had a much easier time using the DSP framework in his derivation but he came up with the model without impulse response or convolution which is quite impressive in itself. I'm gonna check this answer but your mention of Mr. Heavside is appreciated.


They're essentially the same, with the steering vector typically representing phase shifts (==delays) using complex numbers from the unity circle, and the weight vector often containing real values to scale the samples. In practice, both can be unified into one and hence, the terms are sometimes used interchangeably. Also note that for many algorithms, the individual amplitude of the samples doesn't matter, and hence, the steering vector's magnitudes can be chosen arbitrarily (as long as non-zero).


It means that in a multipath fading environment, the signal will have several dominant spatial paths from transmitter to receiver.

From the perspective of a receiver, this means that the received signal will have components from several angles of arrival (AoA).

MIMO techniques exploit this, and use joint beamforming between arrays at both the transmitter and receiver to send different data streams along these spatial paths.  To get the full multiplexing gain, however, these signal paths must be statistically independent.

If the environment does not have any multipath characteristics (i.e. free space), then it is not possible to provide spatial multiplexing, and the best you can do with multiple antennas is conventional beamforming (CBF).


Does sum of random variables follows gaussian distribution although each have different mean and variances?

You can add random variables with different means and variances. For example, look at two continuous uncorrelated normal-distributed random variables $x_1$ and $x_2$ with means and variances $\mu_1, \sigma_1^2$ and $\mu_2, \sigma_2^3$
Then $x_3 = x_1 + x_2 $ is also a normal-distributed random variable with $\mu_3 = \mu_1 + \mu_2$ and $ \sigma_3^2 = \sigma_1^2 + \sigma_2^2 $. For uncorrelated random variables the means and variances just sum.
The distribution of a sum of two continuous random variables is the convolution of the individual  distributions. That's actually the core of the central limit theorem. If you convolve a Gaussian with a Gaussian you get another Gaussian (see for example https://jeremy9959.net/Math-5800-Spring-2020/notebooks/convolution_of_gaussians.html) . If you keep convolving reasonably looking functions with each other eventually it will approach a Gaussian.
A simple example is  summing a few uniformly distributed random variables. The first one has a rectangular distribution. The second one is triangular. The third one already shows a nice round belly in the middle and once you are up to summing 10 or 12 it looks very much Gaussian. That's in fact one method how gaussian random number generators are implemented.
For the Central Limit Theorem, there is no requirement for random variables to have the same mean or variance.


What are the advantages of barker sequence over marker sequnce?

A big advantage of Barker Codes over a generic marker is strong correlation when the codes are aligned and very low correlation for all other shifts, even by one sample.  This offers increased resistance to bit errors that would otherwise cause a false detection.
The presence of a Barker sequence can be optimally detected using correlation, which is the process of multiplying bit by bit the received sequence to the reference Barker code sequence at a given alignment and summing the result of the multiplications (multiply and accumulate). With a Barker sequence, you will get maximum correlation when the received sequence and reference sequence are aligned, and near 0 for all other shifts. This is an ideal characteristic for synchronization, and not a property of generic marker sequences.
Depending on the noise resistance desired, longer codes can be chosen, and/or code sequences can be repeated. Like the use of a generic marker sequence, the Barker code itself would not be an allowed pattern in the data payload, or properly escaped.
As an example, 10110111000 is an 11 chip Barker Code. As an easy way to explain the correlation as a multiply and accumulate, replace all the 0's in the sequence with -1, and all 1's in the sequence with +1. Therefore when correlating the multiplications before summing are an x-nor: 0x0 =1, 0x1=0, 1x0=0, 1x1=0. For demonstration, consider a sequence of three repeating Barker codes for a header:
received sequence
101101110001011011100010110111000
As the sequence is received, we will correlate it by shifting it against our reference sequence and performing the multiply and accumulate for each shift.
Case 1 aligned:
When a reference code is aligned with the received code, the correlation (multiply and sum) is:
101101110001011...
10110111000
1x1=+1
0x0=+1
1x1=+1
1x1=+1
0x0=+1
1x1=+1
1x1=+1
1x1=+1
0x0=+1
0x0=+1
0x0=+1
Sum=+11

Case 2 misaligned by 1 sample:
01011011100010...
10110111000
1x0=-1
0x1=-1
1x0=-1
1x1=+1
0x1=-1
1x0=-1
1x1=+1
1x1=+1
0x1=-1
0x0=+1
0x0=+1
Sum = -1

The same result of -1 occurs for all shifts, except when the code aligns at which point the correlation is maximum (the length of the code).  Since the length of the code is odd, the correlation cannot sum to the ideal value of 0, but gets quite close! Also note that the code can be inverted, and in that case it would sum to -11 when aligned (for a maximum distance of 22 between positive and negative copies of the code).
Maximum Length Sequences (generated with Linear Feedback Shift Registers) also have this shift-correlation property and include significantly longer sequences than Barker codes, which finds application in GPS and Direct Sequence Spread Spectrum communications where signals are received below the noise floor and recovered using these correlation techniques.
Another view from coding theory is the concept of Hamming distance.
(see [https://en.wikipedia.org/wiki/Hamming_distance] 1) To be consistent for that approach we would need to repeat what I did above with symbols 0 and 1 and use an x-or to do the multiply. If the symbols in two sequences are different a 1 is added, if they are the same 0 is added (x-or). The Hamming distance is the total addition of the x-or results between the two sequences.  If every symbol didn't match (maximum negative correlation) the Hamming distance would be 11, and if every symbol matched, the Hamming distance would be 0 (for the 11 chip Barker sequence I used).  Therefore the Hamming distance between 01001000111 and 10110111000 is 11 (maximum distance possible), and the Hamming distance between the code and cyclically shifted versions of the code would be approximately half way in between: For example, the Hamming distance between 01001000111 and 10010001110 = 5. In fact, since the Barker code has noise-like properties, the likelihood is highest for the distance to be near the middle for all random patterns such as what may be represented by our data payload.  It is much more likely that the Barker code will have distances between 4 and 8 to patterns in our data but much less so from 0-3 and 9-11; but that still can an will happen, especially with such a short sequence, unless we encode our data to ensure those sequences do not occur unless bit errors existed. The concept of Hamming distance, and choosing code words of maximum distance is a foundation of error detection and correction.  The bigger the distance, the more errors we can correct (see the reference for specific details on all of that).


In 1D, you are likely to obtain, more precisely, $f_0, f_1, f_2,\ldots, f_{N/2}, f_{-N/2+1},\ldots, f_{-1}$. In other words, the $N/2+1$ ones correspond to positive frequencies. In higher dimensions, the same applies. In 3D, you ought to get a "top-left-front" sub-cube of positive frequencies.

To reduce periodic artifacts, you could 3D window the data before, and take care of the values of the imaginary parts after IFFTing back. If the input data is real, they should be small, and you can keep the real part only. If not, something gone wrong.


Gain at critical frequency for Butterworth and Chebishev filters

If I understood your question correctly, Butterworth filter is maximally flat (no ripples) in the passband frequencies, and at the cutoff frequency, the gain is down by 3dB (3.01 to be exact). But disadvantage is slow roll off from passband. See https://en.wikipedia.org/wiki/Butterworth_filter especially the figure showing magnitude roll off.

For Chebyshev filter type 1 filter, the roll off is more sharp but result is ripples in pass band. You can control this ripple, which will in turn vary the Gp that you mentioned. See https://en.wikipedia.org/wiki/Chebyshev_filter. If you specify $\delta$ as the ripple factor in dB scale (like 1dB), $\epsilon = \sqrt{10^{\delta /10}-1}$ in linear scale. The gain Gp = $\frac{1}{\sqrt{1 + \epsilon^2}}$ in linear scale. The sharpness of your passband and passband ripple specification will affect $\epsilon$ which will inturn affect Gp.


Why non gaussian variables are independent

The model ICA uses says that there exist some unknown, statistically independent sources, $s_i$ that are non-normally distributed (their distributions are something other than Gaussian):

$$
s_i \sim S(\mu_{s_i}, \sigma^2_{s_i})
$$

where $S$ is some (possibly) known but non-Gaussian distribution with mean $\mu_{s_i}$ and variance $\sigma^2_{s_i}$.

Then it is assumed that what you can actually measure is the addition of these:
$$
\mathbf{x} = \mathbf{A}\mathbf{s}
$$
where $\mathbf{s}$ is the vector of the $s_i, i=1,\ldots,N$, $\mathbf{A}$ is an $N\times N$ matrix (usually assumed to be invertible) and $\mathbf{x}$ is the vector of the actual measurements.

Because the $x_i$ will just be weighted sums of the $s_i$:
$$
x_i = a_{i1} s_1 + a_{i2} s_2 + \cdots + a_{iN} s_N
$$
the central limit theorem says that the distribution of $x_i$ will be closer to Gaussian than the distribution of the $s_i$ (provided certain non-restrictive conditions on the true distribution of the $s_i$ are  met).

Then all ICA really tries to do is to find a $\mathbf{W} = \mathbf{A}^{-1}$ the inverse of $\mathbf{A}$.

There are any number of measures for "non-Gaussianity".  One of the simplest (?) is to use kurtosis as the measure of how far the sample of a random variable is from Gaussianity.

So, to attempt to answer:

  So my question is what is non-Gaussianity here and why its necessary to maximize it to extract the original sources.

There are several different ways one could measure non-Gaussianity. For example, kurtosis can be used as it is known that the kurtosis of a Gaussian is 3. Any distribution with a kurtosis different from 3 is therefore "non-Gaussian" to some extent.

The reason we want non-Gaussianity is because we assumed that the original sources are non-Gaussian.


Extremely large reconstruction errors in NMF

When you decompose $X \approx WH$, the columns of $W$ form the dictionary and the columns of $H$ are the time varying weights. In your example, if you choose 150 components, then $W$ should be $33k \times 150$ and $H$ should be $150 \times 1024$.

The error "in thousands" could mean one of two things:

It's practically "zero" if you normalize it by the Frobenius norm of the spectrogram $X$, in which case all is well, or,
Perhaps your sound signal is too rich to be completely representable by 150 components and you should experiment by increasing that number.

It is also worth noting that squared error is often not a good measure for audio perception. Two signals with a large squared difference may sound almost the same.


The frequency response of a filter (or, in general: of a linear time-invariant system) is a complex-valued function

$$H(j\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$

where $M(\omega)$ is the magnitude response, and $\phi(\omega)$ is the phase response. The most straightforward way to interpret those two functions is to look at the way such a system responds to a sinusoidal input signal. For an input $x(t)=A\cos(\omega_0t)$, the response is given by

$$y(t)=AM(\omega_0)\cos(\omega_0t+\phi(\omega_0))\tag{1}$$

So, as you've correctly pointed out, the magnitude defines the filter's frequency dependent gain. The phase response determines the phase shift experienced by a sinusoidal input signal.

Assuming $\omega_0\neq 0$, Eq. $(1)$ can be rewritten as

$$y(t)=AM(\omega_0)\cos\left(\omega_0\left(t+\frac{\phi(\omega_0)}{\omega_0}\right)\right)=AM(\omega_0)\cos\left(\omega_0\left(t-\tau_p(\omega_0)\right)\right)\tag{2}$$

with the phase delay $\tau_p(\omega)$ defined by

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega},\quad\omega\neq 0\tag{3}$$

From $(2)$ we see that a sinusoidal signal of frequency $\omega_0$ experiences a delay of $\tau_p(\omega_0)$. Note that the phase delay $\tau_p(\omega)$ is different from the group delay (cf. this wikipedia article).


Can we do cyclic prefix with zeros in OFDM rather than with samples from end of the ofdm symbol?

There are two major approaches for cyclic extension in OFDM systems - CP (cyclic prefix) and ZP (zero padding, also called Trailing Zeros, TZ). Generally they show the same performance, I mean in AWGN or Fading channel. CP method is the simplest one so it is preferred in the most cases. ZP approach leads to slightly less transmission power level for obvious reason. Doppler performance is the same. The main idea of using ZP as I understood is the existence of equalization techniques for combating spectral nulls that is the major problem of OFDM equalizers. I can't explain this technique because I haven't modeled or implemented it, but you can search some articles in the Web. In one of the article it stands that for CP the simple equalizer exists but symbol recovery isn't guaranteed. But for ZP symbol recovery is guaranteed with complexity of equalizer increasing. I think it is the ckear idea. But I suppose it leads to some matrix subspace algorithms e.g. SVD, so it can be too difficult to implement. I've heard ZP is popular in hydroacoustic OFDM techniques of equalization so you can also refer to this topic.

Hope it helps.


1 - How do I analyze the stability of the last case? AFAIK, I cannot use the margin and allmargin Matlab functions as the delay is not in the open loop transfer function ($G_{comp}G_{process}$).

In both cases the open-loop transfer function is $1+G_{comp}G_{process}e^{-Ts}$

The open-loop transfer function is the transfer function from the output of the summing junction, back to the input of the summing junction.  It includes the effect of all of the blocks -- controller, actuator, plant, and measurement.

You're being mislead by the fact that the closed loop transfer function doesn't include the delay in the numerator -- but it still does.

As to how to make that work with the Matlab functions -- I dunno.  Can they take fixed delays, or must they be presented with rational transfer functions?  If they can only be given rational transfer functions, then you need to calculate the frequency response of the open-loop transfer function and make a Bode or a Nyquist plot and assess stability margins yourself.

  2 - Now assume I'm stuck with a total delay of 100 µs. Say I can split the delay, 50 µs in the measurement delay and 50 µs and the actuator delay. Could I improve somewhat the performance of my controller compared to having an actuator delay of 100 µs?

In this regard, having less actuator delay is fundamentally better, because you have less delay from command to output.  Actually designing a controller that can take advantage of that may be a challenge, however.  Certainly the response to a disturbance cannot be made better -- only the response to a command.  In general pure delays like that just make life hard, no matter what.


The Effect of the Standard Deviation ($ \sigma $) of a Gaussian Kernel when Smoothing a Gradients Image

Let's analyze it in 1D as the intuition is the same.
First, let's have a look on a few different Gaussian Kernels:

As expected, they are wider as the Standard Deviation (STD) increase.
It means that when the kernel is applied using the convolution, more information is aggregates from farther samples. On the other side it means data is spread.
Now, in your images a gradient is a bump. So we can approximate it by a step function:

Let's have a look at the output of the convolution of the kernels from above with the step function:

As we can see, the higher the STD the data is more spread, Though the amount of energy is the same.
In image we expect it to be seen as data is more blurred, wider, while being darker (As the values are lower).
The full code is available on my StackExchange Signal Processing Q70725 GitHub Repository.


Is Fourier series a sampled version of Fourier transform?

There are 4 versions of Fourier transforms that are all close cousins.

It's all due to the basic property that "sampling in one domain corresponds to periodicity in the other domain". If a time signal is periodic, the frequency domain signal is discrete. If a frequency domain signal is periodic, the time domain signal is discrete.

So you have four variants

> `Time                       Frequency               Transform
Continuos & aperidoc          Continuos & aperidoc    Fourier Transform
Discrete  & aperidoc          Continuos &  peridoc    Discrete Time Fourier Transform
Continuos &  peridoc          Discrete  & aperidoc    Fourier Series
Discrete  &  peridoc          Discrete  &  peridoc    Discrete Fourier Transform

The flavors are different since a continuous signal requires integration and a discrete signal needs summation.

The naming conventions are non-intuitive: it would be easier if the DFT would have been called Discrete Fourier Series and the DTFT called DFT, but it is the way it is.

Which flavor you need to use depends on the nature of your signals and depending on what assumptions around the signals you can make, you can use one flavor to calculate a different one.


How to make pop-free, crisp two-beep ringtone?

The "click" is caused by the discontinuity in the waveform and its derivatives - even if the waveform stops at a zero-crossing you might still hear a pop if there's a discontinuity in the higher order derivatives! Its loudness depends on the amplitude at the discontinuity - and it is thus influenced by the length/frequency/phase of the sinusoidal tone. Note also that if the tone is contained within a clip in an audio editor, a short fade-in/fade-out might be added automatically by it, attenuating the problem. I know a few audio editors that dynamically apply a short fade (say 32 samples) at the boundary of each audio clip during mixing during playback.
You hear the click clearly with the sine wave because the tone has very few harmonics (ideally: none). You might not hear it with the other waveforms (square) because they already have high frequency content that will cover part of it (frequency masking phenomenon). For example, you won't hear a click when a burst of white noise stops.
The solution is to apply a short fade-in / fade-out. 1ms is a good start. Above 5ms it'll start feeling "slow". There is a trade-off to find between "clicky and snappy" and "slow" - this is a problem even well-established synth manufacturers have to educate the public about. With a sine-wave, an envelope curve like a raised cosine function ensures that the resulting signal has continuous derivatives throughout the fade.
Aliasing is the bad thing to get rid off, so it is normal that the "no-aliasing" option produces a better quality output.
A sine wave with an amplitude of 1V has twice as less energy as a square wave with the same amplitude. Assuming that the phone normalizes amplitude, it is normal that with a given dynamic range/amplitude budget, the sine wave will sound quieter. Energy/RMS is a rather bad proxy for the perceptual concept of loudness, but it's much more accurate at this task than peak-to-peak amplitude!


How did they color old black and white movies?

There is no way of recovering the original color information from a black and white photo, so whether Einstein (resp. Waheeda Rehman) was wearing a pink or green sweater (resp. Dupatta) is up to your imagination.

Historically, this has been done by hand, by painting over the film. The first digital techniques to automate the process consisted in "painting" a few dots of color on each frame, at the center of each uniformly colored region, and using something like a voronoi partition + some blurring to get a color map for each frame (see for example US patent 4606625).

Today, this can be done relatively easily (though manually) with video editing software, by using vector masks to indicate regions of uniform color on a few keyframes, and interpolating between them. Then a color transform is applied to each mask. See it in action here.

Standard image segmentation and region tracking techniques can be used to automate the task of segmentation and the marking of regions on each keyframe - for example by propagating manual annotations to similar/adjacent pixels in space/time, or by detecting uniformly textured regions. Texture and gray level similarity can be used to propagate color cues from a color image to a greyscale image depicting a similar subject - in this case the manual process only consists in finding a template color image - this later task can be itself automated using content-based image retrieval techniques.


A system is stable if $ \int_{-\infty}^{\infty}  |h(t) |dt < \infty$

In your case $ \int_{-\infty}^{\infty}  |h(t) |dt = \int_{-\infty}^{\infty}  |e^t \sin(-3t)u(-t)| \ dt$

Since $\sin(-3t) = -\sin(3t)$,

$ = \int_{-\infty}^{\infty} | e^t \sin(3t) u(-t)| \ dt$

$ = \int_{-\infty}^{0}  e^t |\sin(3t)| \ dt$ , $u(-t) = 1$ if $ -\infty < t < 0$ else $0$

$ < \int_{-\infty}^{0}  e^t \ dt < \infty$

So the system is stable.


What if they did all use the same carrier?

What you're describing is a single-frequency network (SFN). These are in common use for things that are not stupid FM broadcast.

The whole truth is that reception from different senders in an SFN just look like heavy multipath propagation, where the sent signal just takes multiple paths of different length to the receiver, and thus the receiver sees the signal as what would come out of a tapped delay-line.

If you consider it like that, the effect is clear: the different delays with which the signal can reach the receiver lead to something that for all practical purposes is an FIR filter. If you look at it from the viewpoint of digital communications, you get Inter-Symbol Interference (ISI) as delayed versions the last sent symbol still reach the receiver while the shorter paths already brought on the next symbol.

If you look at the same thing from the perspective of spectrum analysis, it is a frequency-selective channel.

Either way you look at it, you counter that with an equalizer that requires knowledge of how the channel looks (i.e. the different propagation delays, or the frequency-dependent behaviour of your channel, which are just two descriptions of the same thing, linked by the Fourier transform) and can then do its best to reverse that effect.

Now, sadly, FM radio is not digital communications. You don't get symbols, and by the very nature of FM  being a modulated carrier, you don't even get any kind of "known, not very repetitive" input with which you can estimate what your channel looks like.

This needn't be a problem – if the resulting delays don't make the channel too harsh, you don't lose much SNR at the receiver. We can measure that by comparing the delay spread (difference between shortest and longest path) to the inverse of the channel bandwidth. If the channel bandwidth is "sufficiently" smaller than the inverse delay spread, we're probably in luck and can still consider the channel "sufficiently" flat. "Sufficiency" here depending on your needs, so let's just do a order-of-magnitudes calculation.

An FM channel is roundabout 150 kHz wide. Invert that, and you get $6\frac23$ µs.

Electromagnetic waves travel with the speed of light in vacuum, and for all practical purposes here, in air, too. So those microseconds correspond to a distance difference to two transmitters of $c_0\cdot \frac{20}3\cdot10^{-6}\,\text{s}=3\cdot 10^{8}\frac{\text{m}}{\text{s}}\cdot \frac{2}3\cdot10^{-5}\,\text{s}=2\,\text{km}$.

Since we can be pretty sure that different FM radio transmitters aren't only 2 km apart, we can directly see that the resulting radio channel will be wildly non-flat and frequency selective. In fact, since the distance between these will be orders of magnitudes larger, we can assume that even given good channel information, we might have frequency-dependent deep fades that we simply can't properly "equalize away" (you get nothing through at those frequencies, so you can't even amplify anything to make it "flat" again).

In other words, FM and SFN don't play well together. There's been attempts to still implement SFNs for FM broadcast (I think Netherlands tried that), but as far as I know, they failed – they can only work by "aligning" the interference for limited locations, and if your network can't be too large and there's no strong variation of propagation environments for different transmitters (which would explain why small and flat Netherlands was an attempted adopter).

That's one of the most important reasons the world is leaving analog broadcasting behind. DAB and DAB+ are SFN-based, as their heavily multi-carrier oriented OFDM system with channel bandwidths of 1 kHz allows for much larger delay spread, and their frame structure and channel coding allow for heavy error correction, so that even deeply faded subcarriers don't matter that much.

So, yes, you've found a (if not the most important) spectrum management reason why FM should, if you ask me, die as soon as possible: it enforces the non-reuse of spectrum, and that's a waste of a precious resource.

Now, you mention multiple antennas:

Well, yes, if you have multiple antennas, you can use phasing to make them have some directivity, i.e. you can "simulate" having a single antenna that only receives from one direction. Hm. You buy multiple antennas, build a complex and large phased array and try to operate in a SFN, and the best-case scenario is you get the same performance as in a MFN with a single antenna. Doesn't sound too tempting, and, you'd need to have an idea in which direction you want to "look". It's notoriously hard to estimate that out of, again, an analogue signal with no known non-periodic structure. In a situation where you're at a fixed position (moving obviously breaks knowledge of directions) and adjust antenna phases just to get broadcast radio: Have you considered satellite radio?


Maximum Likelihood Estimation in Presence of Colored Noise

Samples of colored noise (taken at different times) generally
are correlated random variables because the
autocorrelation function of the noise process is not a delta function
as it is in the case of white noise. Thus, if we assume a zero-mean
process (noise is generally assumed to be regardless of its color),
then the covariance of two signals separated in time by $\tau$ seconds is
$R(\tau)$ where $R(t) = \mathcal F^{-1}(S(f)$ is the autocorrelation
function of the process (inverse
Fourier transform of power spectral density). Note that is is possible
for $R(t)$ to be zero for some values of $t$ (e.g. $R(t) = \operatorname{sinc}(t)$ is a valid autocorrelation function),
but it cannot be zero
for all nonzero $t$.

As far as the density function of any sample, if the process is
Gaussian, the sample is Gaussian even if the process has been filtered
with a linear filter before sampling. But if the process is not
Gaussian (it is,
let us say, LaPlacian), then while each sample will be LaPlacian, the
same cannot be said generally of samples of the process after filtering
of any kind. In other words, Gaussianity
survives linear filtering, LaPlacism generally does not.

So, how does maximum-likelihood estimation work when samples have
correlated noise? Consider the case when we wish to estimate the
unknown mean of a $\mathcal N(\mu, 1)$ random variable, and we
have two observations $x$ and $y$. In the standard case of
independent observations, the likelihood function is
$$L(\mu) = \frac{1}{2\pi}\exp\left(-\frac{1}{2}\left[(x-\mu)^2+(y-\mu)^2\right]\right).$$
The _maximum-likelihood estimator for $\mu$ is the number $\hat{\mu}$
that maximizes $L(\mu)$, which works out to be the number $\hat{\mu}$
that minimizes $(x-\mu)^2+(y-\mu)^2$. This is a quadratic in $\mu$
and the maximum-likelihood estimate turns out to be
$\hat{\mu}=\frac{x+y}{2}$. When the observations are correlated with
correlation coefficient $\rho$, then
$$L(\mu) = \frac{1}{2\pi\sqrt{1-\rho^2}}
\exp\left(-\frac{1}{2\sqrt{1-\rho^2}}
\left[(x-\mu)^2-2\rho(x-\mu)(y-\mu)+(y-\mu)^2\right]\right).$$
Once again we need to find the $\hat{\mu}$ where
$(x-\mu)^2-2\rho(x-\mu)(y-\mu)+(y-\mu)^2$ has a minimum.
We still have a quadratic in $\mu$ but now we get terms
like $xy$ in the coefficients. What $\hat{\mu}$ works out to
be is left for you to work out.

And what if we have $n$ observations where $n > 2$? All the above
still applies. For independent identically distributed Gaussian noise
in the samples, the sample mean
$n^{-1}\sum_i x_i$ is the maximum-likelihood estimate of $\mu$ but
in the case of correlated Gaussian random variables we get very
messy minimization problems because the quadratic that we are trying
to minimize depends on the inverse of the covariance matrix and
the result is a nonlinear function of the data instead of a simple
easy-to-remember result like the sample mean.

What if the noise is not Gaussian? The same principles apply -- set up the
likelihood function and find where it attains its maximum value --  but
the calculations are quite a bit different, all depending on
what you assume or know is the joint density of the observations.


How to see the amplitude of each frequency?

As I have not had enough reputation to comment, I will provide some of my interpretation of your question (Don't need to accept my answer).

It looks like that you want a spectrogram represented by a 500$\times$22000 matrix. This relates to the number of time and frequency slots. The number of time slots (columns) is determined by the sample number of the given signal. The number of frequency slots (rows) is determined by number of FFT output, which you can adjust by methods such as zero padding.

If you want to get better-looking spectrogram, you might need to adjust the window that is used in performing STFT (or maybe more overlapping in STFT). However, sometimes you just cannot get satisfying resolutions as the time-frequency resolutions are bound to the uncertainty principle. This is decided by your time-domain samples.


Reducing noise from the same frequency band as signal

Once you resample to the common Nyquist rate (or nearly), you can consider it a mixture model. Iirc, the optimal detector for a Laplacian signal in Gaussian noise is the soft threshold function.


AM-AM plot (input power vs output power) for power amplifier

I would assume that you are only monitoring the power in-band, and what is then occurring as the signal goes into saturation is that power is being translated to other frequencies due to the non-linearity. I would suspect that this would dominantly be power in harmonic frequency locations, notably the third and fifth harmonics.


Applying Superposition Property : $x^{2}[n]$ vs $x[n^{2}]$

Given the system I/O definition:
$$y[n] = \mathcal{H}\{x[n]\} = x[n^2] \tag{1} $$
you can easily show that it's a linear (but time-varying) system.
Following the standard procedure
let
$$y_1[n] = \mathcal{H}\{x_1[n]\} = x_1[n^2] \tag{2.1}$$ and
$$y_2[n] = \mathcal{H}\{x_2[n]\} = x_2[n^2] \tag{2.2}$$
then define
$$x_3[n] = a x_1[n] + b x_2[n] \tag{3}$$
and
$$
\begin{align}
y_3[n] &= \mathcal{H}\{x_3[n]\} \tag{4}\\\\
&= x_3[n^2] \tag{5}\\\\
&= a x_1[n^2] + b x_2[n^2]\tag{6}\\\\
&= a y_1[n] + b y_2[n]\tag{7}\\
\end{align}
$$
where Eqs(4)-(5) follow Eq.(1), Eq.(6) follows Eq.(3), and Eq.(7) follows Eq.(2).
Equation (7) indicates that the system is linear.


Interpretation of system response to unit impulse

In the figure that you have posted, the right plot shows the impulse response $h[n]$ according to the convention that only its nonzero samples are displayed. And therefore, it's understood that all the remaning samples of $h[n]$ are zero unless otherwise stated explicitly.

This means that $h[0]=0$, $h[1]=0$,...,$h[12]=0$. Which can be stated in the following:

$$ h[n] = \begin{cases}
0 &,& \text{ for } n \le 4 \\
\alpha_n &,& \text{ for } 5 \le n \le 11 \\
0 &,& \text{ for } n \ge 12 \\
\end{cases}
$$

In this setting the domain of support for the finite length sequence $h[n]$ is stated as $n \in [5,11]$ , which means that $h[n]$ is certainly zero for $n$ out of this interval, while it can have nonzero as well as zero values inside this interval.


Can edge detection be done in the frequency domain?

Since convolution in spatial domain is multiplication in the Fourier (frequency) domain, you can perform edge detection in Fourier domain by multiplying the spectra of image and the edge detection kernel and then perform IFFT on the result.

I think the high-pass filter alone is not appropriate for edge detection since it keeps all high-frequency features (e.g. sharp peaks and corners) which are usually not classified as edges.

More advanced edge detection methods would be tricky in frequency domain since edges are best described in spatial domain (in my opinion).

The question is why to do edge detection using FFT in the first place? Is it because performance considerations? If so, maybe the high-pass filtered image (produced quickly by FFT) can be quickly filtered again to remove non-edge parts.


What is the general formula for radon back projection for a javascript implementation?

To implement projection the simplest way is to rotate your image then sum over a row or column. The simplest way to implement back projection is to take a line of your sinogram which is a projection from certain view (angle) then repeat that line to form an image then rotate that image regarding the projection view, then sum all of backprojected views.
There is a unique formula for continuous case but for discrete case there are many different implementation, this one is the simplest.

Good luck.


It appears you're getting a bit lost between the time and frequency domains. If "Num" is the numerator of an FIR filter (as it appears), these are the coefficients (or impulse response) of your filter. The impulse response as a signal exists in the time domain.

If you want to perform a filtering operation using this impulse response (or filter kernel as some call it), you should perform a convolution (conv) between it and the signal you wish to filter in the time domain.

--

The impulse response you've created has an equivalent frequency domain response (which you will have seen in fdatool).  The time domain convolution described above is equivalent to multiplying this frequency-domain response by the frequency-domain representation of your input signal.

Since your signal x is complex I assume you mean it to a frequency-domain representation of some 'noise' signal. This representation in itself is strange, and it doesn't represent a 'white' noise signal, but a purely random arrangement of frequency and phase that I find it difficult to conceive.

If you want to observe low pass filtering I would suggest the following:

Create a pure real (no 'i' component) noise signal in the time domain using rand (remember that rand produces values in the range 0->1, you will need to scale and shift it to create a signal in the range -1 -> 1).  This signal will in theory be purely white, and hence have a flat magnitude response in the frequency domain (i.e. it would be a straight line).
Convolve this signal with your filter impulse response.

You could consider this operation to be equivalent to multiplying the response of your filter (that you saw in FDATool..) by the frequency domain of the input signal (which is effectively constant). The result should be a scaled version of your familiar low-pass filter response (although some other effects such as the trade-offs in your filter design come into play).

The output of this procedure is still in the time domain, so if you perform an FFT or some other transform to convert to the frequency domain, you can then plot it, and observe the low-pass filtering effect.

Good luck for your learning :)


A continuous time signal can be evaluated arbitrarily at any time,at least in theory, and can take on values at arbitrary points between any evaluation point. The Fourier transform, while not always well defined, can take on values at frequencies over $-\infty$ to $\infty$ and that value can be zero, and the signal in continuous time also takes on values over the infinite interval and those values can also be zero.

A discrete time signal only takes on values at discrete time intervals. It is only meaningful to consider values intermediate in time, intermediate to those  samples,  if for example the discrete time signal was sampled from a continuous signal, or if you have a reason (justification). In general, the values intermediate to the samples are not defined. You can consider the price of a stock or pork bellies on days where there are no trades. No one is buying or selling. The price isn’t the last traded price and it isn’t zero. It does have a Fourier Transform, computed as an infinite summation, and the frequency is uniquely (continuously defined) limited over $-\pi$ to $\pi$.

So there is only a subset of discrete time signals that we can meaningfully compare their Fourier Transforms, and those are signals that are Nyqvist Sampled.  Band limiting and periodic sampling alter the original continuous time signal, which alters the frequency spectrum, which almost answers your question.

The question of whether something is meaningful at intermediate points between samples in both time and frequency is more subtle.

I can have a periodic discrete time sequence and the frequencies in the Fourier transform are likewise values that are discrete in frequency, but are the intermediate values in frequency zero?

The most used DSP text book is  Oppenheim and Schaefer and the most recent edition has 3 interpretations of the DFT which I interpret as yes , no, and maybe.

I can also interpolate pork belly prices, but is it meaningful?

I think the answer lays with how you want to represent something.

If you want to describe a discrete time signal in the continuous time domain, a sequence of Dirac Delta functions is suitable where intermediate values are zero is OK. In other contexts, it may be meaningless.

Another example worth considering is a realization of a Poisson process.  You have a set of "arrival times" that you can store on a computer, limited by the amount of memory and numerical precision available on the machine. It is a list of numbers and the intermediate values denote no "arrivals" at those intermediate times.  The Poisson process is continuous time.  One typically doesn't perform Fourier analysis on a Poisson Process. There is some discretization in the computer representation, so one can quibble that it isn't strictly continuous but it is a practical representation of a continuous signal.  This is also an example of a continuous signal where Nyquist sampling is irrelevant.


The FFT is just an implementation of the DFT, which is an invertible operation on any vector composed of complex numbers.
Since invertible operations need to be surjective:
No, an FFT can't have less outputs than inputs.
You are, probably, looking for something that "condenses" your data, e.g. looking for features, or projecting onto a smaller vector space that is sufficient to describe your signals of interest; keywords you might be looking for would be "feature extraction" or "compressive sensing".
However, I think it might be sensible for you to post a new question, describing exactly what you actually want to do (and not focused on employing the FFT itself).


Along the line as Fat32's suggestion to consult IEEE papers, you might try to obtain a copy of "Digital Signal Processing" (IEEE Press Selected Reprint Series), Edited by: Lawrence R. Rabiner and Charles M. Rader, published by: IEEE press (1972).

That book contains reprints of many early (starting around the middle 1960s), and influential, scholarly papers on digital filtering. And of course those early papers contain references to even earlier papers.


Comparing two Accelerometer signals in time domain

I suggest using a correlation metric as a comparison. Remove any mean in the signals and after subtracting the mean, compare the signals using the Pearson Correlation Coefficient to have a single data point that scales between +1 and -1 with regards to their direct correlation assuming any time offset factors have been removed, and even better as a first evaluation use the cross-correlation function which would provide the above (without normalization) for every offset in time, in case their are delays between the results.


If $x(t)$ is deterministic its variance is zero. Note that the mean of $x(t)$ (defined as an expectation) just simply equals $x(t)$:

$$\mu_x(t)=E\{x(t)\}=x(t)\tag{1}$$

where the last equality follows from the fact that $x(t)$ is deterministic. From $(1)$ we get the variance

$$E\{[x(t)-\mu_x(t)]^2\}=E\{[x(t)-x(t)]^2\}=0\tag{2}$$


Can anybody recommend a good source for material of Gaussian and Complex Gaussian wavelet?

The notion of "Gaussian wavelet", whether real or complex, is not clear from the literature. Indeed, the Gaussian has a positive average, while wavelets are (or should be) zero-mean (to satisfy the admissibility condition).

Two standard wavelets may suit your need:

the Morlet wavelet, or Gabor wavelet, a complex exponential  multiplied by a Gaussian window envelope (or its real part),
derivatives of Gaussians (including  differences of Gaussians or DoG), sometimes called Ricker wavelet in 1D in geosciences or mexican hat in 2D.

Their complex extensions, or real restrictions, are relatively natural.


Characteristic and moment generating function of a random variable interpretation

As you well said, these functions are useful just for their mathematical properties that allow to compute the m-th order moment of a random variable.

Additionally, they happen to match Fourier transform theorems. However, not in the sense of a time-domain signal, since the variable in the probability density function (pdf) is the random variable itself, and not time.
A frequency interpretation only makes sense in time-dependent functions. Actually the definition itself of frequency is how often something happens.

Anyway, one could consider the "bandwidth" for $\Phi_X$ in the $\xi$-domain as how spread the pdf is around the mean.
A more spread pdf would imply a narrower "bandwidth", and less spread of the pdf would imply a wider "bandwidth".
For example, consider a uniformly distributed random variable:
$$f_X(x) = 1/\delta~\Pi(x/\delta)$$
$$\Phi_X(\xi) = sinc(\delta~\xi)$$
As $\delta$ gets smaller, the random variable is constrained in a narrower interval, while the "bandwidth" in $\Phi_X(\xi)$ gets wider, and vice versa (measuring the "bandwidth" as the width of the main lobe, $2/\delta$, as we do for communication signals).

Just by observing the characteristic function one may have an idea of how spread is the pdf. By comparing two characteristic functions with different "bandwidths" one may determine which pdf is spread the most.


I wish to create a GNU Radio program in C++ that doesn't need GRC or Python, is there a reference or template for this?

you can create flowgraph in c++ like gqrx and i think GNU Radio Manual and C++ API Reference documents will help you. This is gnuradio c++ document for top_block if you could not use it let me know to write a sample for you. These documents are for gnuradio version 3.7 and in gnuradio 3.8 there is c++ code generation and gnuradio blocks are porting to c++ and you can track them in this link.
i hope it helps.


The Fourier transform of $x(t)$ can be easily obtained if you realize that $x(t)$ can be written as the convolution of the pulse $R_T(t)$ and a Dirac comb:

$$x(t)=R_T(t)*\sum_{k=-\infty}^{\infty}\delta(t-2kT)\tag{1}$$

Since convolution corresponds to multiplication in the frequency domain, you get for the Fourier transform of $x(t)$

$$X(f)=\mathcal{F}\{R_T(t)\}\cdot\frac{1}{2T}\sum_{k=-\infty}^{\infty}\delta(f-k/2T)\tag{2}$$

because the Fourier transform of a Dirac comb is a Dirac comb (see e.g. this table). And since

$$f(t)\delta(t)=f(0)\delta(t)$$

Eq. (2) can be written as

$$X(f)=\frac{1}{2T}\sum_{k=-\infty}^{\infty}\hat{R_T}\left(\frac{k}{2T}\right)
\delta(f-k/2T)\tag{3}$$

where $\hat{R_T}(f)$ is the Fourier transform of $R_T(t)$. Eq. (3) shows that the Fourier transform of $x(t)$ is a weighted sum of Dirac impulses, where the weights are equidistant samples of the Fourier transform of the pulse $R_T(t)$.

Since the low pass filter $H(f)$ only passes frequencies lower than $3/4T$, only the contributions for $k\in\{-1,0,1\}$ of the sum in (3) are passed. This gives the final result for $Y(f)$. Note that there are two typos in your formulation of $Y(f)$. It should be

$$Y(f)=\frac12\delta(f)+\frac{1}{\pi}\delta(f-1/2T)+\frac{1}{\pi}\delta(f+1/2T)\tag{4}$$


How to calculate STFT of a function for a rectangular window

I don't understand what you mean by "positioned at $n=5$" with length $20$. Indeed:

is "positioned" the beginning  or the center?
is length an half-length around a central position, or a full length?
finally, at a given location, this is not an STFT anymore, but a mere windowed Fourier transform.

If I interpret your question in its most obvious sense, the window starts from $n=5$ to $n=20+5-1=24$, on an interval where $u[n]=1$, so you'd just have to compute a simple DFT of a cosine. The result can be obtained via Euler/De Moivre formulae, with two finite sums of geometric series.

If this is not the case, please provide more information, close forms are at hand