Matlab: Help in understanding if the example for Maximum Likelihood Equalizer is properly functioning

I had a look at your code, and I know why you get 0 BER.

The thing is that you do not add any noise, you should use Eb/N0 to calculate the noise variance, and add a complex noise on your channel output.

The distortion added by "chanOutput = filter(chCoeffs,1,modSignal);" merely adds linear ISI, by convolving the channel impulse response with the input signal. This linear distortion is exactly known to the equalizer since you pass it the exact channel coefficients. This is never the case in real applications.

Even assuming you do know the channel coefficients, you still need to add noise immediately after the channel distortion and before the equalizer, to get a meaningful result of BER over Eb/N0.

Noise can be generated with randn having unit variance, then you need to multiply this noise by a constant.

I did not spend too much time understanding what you want to verify, but i did correct the script to include noise.

clear
N = 10^4; % number of symbols
N2= 10;   %repeat with same Eb/N0

% Eb_N0_dB = [0:25]; % multiple Eb/N0 values
Eb_N0_dB = (0:10); % multiple Eb/N0 values
nTx = 1;  %SISO
nRx = 1;

M = 4; %QPSK

for ii = 1:length(Eb_N0_dB)

    %% ----- ADDED
    EbN0 = Eb_N0_dB(ii);
    SNR_linear  = 10^(EbN0/10)*log2(M);
    SigPWR      = 1; %unit circle symbols
    NoisePWR    = SigPWR/SNR_linear;
    NoiseFactor = sqrt(NoisePWR/2); %half power complex half imaginary

    hMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
    hDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
 % Channel coefficients
    chCoeffs = [.986; .845; .237; .12345+.31i];
    hMLSEE = comm.MLSEEqualizer('TracebackDepth',10,...
                      'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
 % Create an error rate calculator
    hError = comm.ErrorRate;
    for n = 1:N2
      data= randi([0 3],N,1);
      modSignal = step(hMod, data);

 % Introduce channel distortion.
      chanOutput = filter(chCoeffs,1,modSignal);

      %% ----- ADD NOISE !!
      noiseVec   = NoiseFactor*randn(N,1)+1j*NoiseFactor*randn(N,1);
      chanOutput = chanOutput+noiseVec;

 % Equalize the channel output and demodulate
      eqSignal = step(hMLSEE,chanOutput);
      demodData = step(hDemod,eqSignal);
 % Compute BER
      errorStats = step(hError, data, demodData);
%       nErr=errorStats(1) ;
    end
    fprintf('Error rate = %f\nNumber of errors = %d\n', ...
      errorStats(1), errorStats(2))

    %assume 1 symbol error = 1 bit error (needs Gray coding in real system for this to be true)
    SER(ii) = errorStats(1); % SER ~= BER with undelrying Gray coding
end

%% changed, this is already divided by N*N2
simBer = SER; % simulated ber

EbN0Lin = 10.^(Eb_N0_dB/10);
theoryBer_nRx1 = 0.5.*(1-1*(1+1./EbN0Lin).^(-0.5));
p = 1/2 - 1/2*(1+1./EbN0Lin).^(-1/2);
theoryBerMRC_nRx2 = p.^2.*(1+2*(1-p));

close all
figure
semilogy(Eb_N0_dB,theoryBer_nRx1,'bp-','LineWidth',2);
hold on
semilogy(Eb_N0_dB,theoryBerMRC_nRx2,'kd-','LineWidth',2);
semilogy(Eb_N0_dB,simBer,'mo-','LineWidth',2);
axis([0 25 10^-5 0.5])
grid on
legend('theory (nTx=1,nRx=1)', 'theory (nTx=1,nRx=2, MRC)', 'sim (nTx=2, nRx=2, ML)');
xlabel('Average Eb/No,dB');
ylabel('Bit Error Rate');
title('BER for BPSK modulation with ML equalizer');


Separating/recovering base signal from two mixed signals, given phase information

Here's the best method I've found so far:
First, I shifted $B_2(x)$ by $\phi$; this shift aligns the underlying (but unknown) signals $w(x)$ and $w(x + \phi)$. Note that once $B_2(x)$ is shifted, the signals will not begin and end at the same $x$-values. I chopped off the values of each signal that didn't overlap. Parts (a) and (b) of the figure below show the original measured signals $B_1(x)$ and $B_2(x)$, and then the aligned signals, which I'll call $s_1(x)$ and $s_2(x)$.
The aligned signals $s_1(x)$ and $s_2(x)$ are not identical, because the base signal $b(x)$ has been added to them at a known phase shift. Therefore, finding the difference between $s_1(x)$ and $s_2(x)$ at each point, and then dividing by the phase shift $\phi$, will give us the numerical derivative of the base signal, $\frac{\triangle b}{\triangle x}$. In summary, the numerical derivative of the base signal is given by
\begin{equation}
 \frac{\triangle b}{\triangle x} = \frac{s_2(x)-s_1(x)}{\phi}.
\end{equation}
Then, to find the estimated base signal $b_{est}(x)$, we numerically integrate $\frac{s_2(x)-s_1(x)}{\phi}$. I corrected for integration drift by numerically integrating forwards and backwards, then averaging the two results.
Lastly, I shifted $b_{est}(x)$ forward by $\frac{\phi}{2}$. I think the reason for a shift of this value is because I computed the numerical derivative between points separated by a full phase, so the derivative approximate should be centered at half of that.
The results (base signal $b(x)$, estimated base signal $b_{est}(x)$, and error between the two) are shown in part (c) of the figure below.


Simplest way of detecting where audio envelopes start and stop

What you really want to do is essentially called as Voice Activity Detection or speech detection.

Basically any pure speech signal (which contains no music) has three parts.

The voiced sound - which is basically caused by Vowels
The unvoiced sound - which contains consonants.

The characteristic of human sound is such that while a lot of energy is used in voiced sound the real information is contained in consonants. Also, voiced sound is usually lower frequency where as unvoiced sounds are higher frequencies. [To be precise all voiced sound are resonated more or less a constant frequency for a given person which is his/her pitch].

Now, as any system there is noise. The voiced sound is usually quite powerful enough that it can be distinguished visible. When you apply a lower frequency filtering it is possible to collect good magnitude of voiced sounds however, the unvoiced sound (with all the rich information) will get lost.

Coming to the question how to solve it:

The trick lies in the fact that unvoiced sound still come from a resonating source; and inherently restricted over a certain frequency. Where as, the noise is rather uniform. So a simple measure that distinguish all three is "local power" or alternatively but equivalent is to take the windowed auto-correlation.

If you take at a time say 100 samples - and auto correlate itself, if it contains only noise the results will be pretty much zero (this is the property of white noise) where as for speech signal, this magnitude will be observable because the signal still has better structure.  This has worked for me in the past.

VAD has been an active research areas- because almost all Mobile phone communications wants to detect non speech part and remove them from encoding. But if they would remove non-voiced speech this would make telephony useless.

The G.729 standard computes VAD based on features like: line spectral frequencies, full-band energy, low-band energy (<1 kHz), and zero-crossing rate.

The GSM standard works as follows: Option 1 computes the SNR in nine bands and applies a threshold to these values. Option 2 calculates different parameters: channel power, voice metrics, and noise power. It then thresholds the voice metrics using a threshold that varies according to the estimated SNR. (from wikipedia)

For more advanced techniques i am listing some references on this subject.

Most sited reference: Jongseo Sohn;   Nam Soo Kim;   Wonyong Sung;  "A statistical model-based voice activity detection" Signal Processing Letters, IEEE, Jan 1999, Volume: 6 Issue:1 pp:1-3
Most relevant for you: Mark Marzinzik and Birger Kollmeier "Speech Pause Detection for Noise Spectrum Estimation by Tracking Power Envelope Dynamics" IEEE TRANSACTIONS ON SPEECH AND AUDIO PROCESSING, VOL. 10, NO. 2, FEBRUARY 2002 pp.109
Ramírez, J.; J. M. Górriz, J. C. Segura (2007). "Voice Activity Detection. Fundamentals and Speech Recognition System Robustness". In M. Grimm and K. Kroschel. Robust Speech Recognition and Understanding. pp. 1–22. ISBN 978-3-902613-08-0.
Introductory : Jonathan Kola, Carol Espy-Wilson and Tarun Pruthi "Voice Activity Detection"


Very interesting problem
I will focus on the second aspect of your question, how to improve it.
What I would do next, based on the experience you shared is to combine the two methods.
One problem of your first approach is that it accumulates the errors in the alignment between multiple images. A possible approach would be to use your focus measure to align the images iteratively.
$$ F_{i} = F_{i-1} + translate(image_{i}, o_{i}) $$
$$ o_i = \arg\min_{o}\quad  focus\left[ F_{i} + translate\left(image_{i}, o\right) \right]$$
This approach is maybe not the best at the beginning, but after the first step it have the advantage of not accumulating errors in the same way as if you compute $o_i$ as using only $image_i$ and $image_{i-1}$. Notice that $F_{i-1}$ aggregates all the previously aligned images, so if the one image is not aligned very well it will make so much difference since it sees how it was aligned before.
This is maybe not going to be as good as your global method. What I would try, to improve your global method is to make the displacement continuous with some type of interpolation, because a continuous optimization problem is much easier to solve than a discrete optimization problem, because can use the gradient information[1]. I am sure you are not doing the search for all the offsets independently that would be $s^N$, where $s$ is the number of possibilities in your search list.
Also instead of limiting the the displacements to a given range I think you could simply add some cost to the differences of the displacements $k \sum ||o_i - o_{i-1}||^2$, this makes sure that if have a high relative displacement between two frames the cost function will indicate that this is unreasonable.


No, if you do zero-padding of $H_1$ in frequency domain.

But yes if you do zero-padding of $h_1$ in time domain.
Basically you should have two FFT versions of $H_1$ ($H_{1,512}$ and $H_{1,1024}$). Because $h_1$ is static, you can precompute both, or you can compute $H_{1,1024}$ only when required.

Let $a = IFFT(FFT(x) \times H_{1,512})$ and $b = IFFT(FFT(x) \times H_{1,1024} \times H_2
 \times H_3)$. Size of a is 512 and size of b is 1024. For example, using overlap add method, we can freely change fft-size between 512 and 1024, but it still process 256-input per fft.

===a====
     +
    ===a====
         +
        ===b============
             +
            ===b============
                 +
                ===a====


Isn't the rect pulse filter implicit at the transmitter?

In most waveforms for data communication we do NOT want to transmit a rect for each pulse given the excessive bandwidth required (since a rect is a Sinc function in frequency which is relatively very wide-band). This is the motivation for pulse shaping (such as the commonly implemented root-raised cosine pulse shape): to limit the transmitted spectrum and improve the overall spectral efficiency of the waveform.


How to average coherences estimated using Welch's method in MATLAB

If you look at the documentation of mscohere which contains the definition of magnitude squared coherence, you will see that it is the absolute squared of the cross spectral density, divided by the product of the two separate spectral densities. Since power densities enter numerator and denominator of the expression, the MSC estimate from a full time series is not simply the average of the estimates from the parts.

To get a good "average" MSC estimate, do the following: estimate in each segment separately cross and single spectral densities, average these three over the segments, and then combine them in the end into an average MSC. Continuing your code:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);
sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);
sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;
sxya = (sxy1 + sxy2) / 2;
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

Now max(abs(c-ca)) is almost zero.


To stand on @hotpaw2's answers, think of Mel as one kind of pyscho-acoustic scale, derived from a set of experiments on human subjects, others are Bark & ERB

Why have such a scale?

Imagine a 100 Hz sine wave playing in your head, .... wait for it to stick ...okay... Good.
Imagine now, a 200 Hz sine wave playing ... marinate in it. Good work.
Now compare the two, they have some "perceptual distance". Perceptual only because we're using your soft-grey matter as a measuring instrument.
Finally, if you were you repeat this experiment with 1100 Hz & 1200 Hz, and again with 10,100Hz and 10,200Hz, your opinion of this percpetual distance might diminish, i.e. 100 & 200 Hz can sound "farther apart" than do 10,100 & 10,200Hz, even though in all cases the differences are equivalent to 100Hz.

These pyschoacoustic scales try to capture these distances from low (20Hz) to high frequency (20kHz).


ICA - Exercise Problem from book - Matlab Code

In the line

C=(mnp*mnp')/3

you're dividing by 3 when you should really be dividing by the number of observations, which is 1000 in your case.


No can do. Your scene is not not planar, nor the camera motion you want to simulate is a pure rotation about the camera's center (a "pan-tilt"). These are the only cases where an image transformation with a pure homography will yield a realistic result.

One thing you can do is select one of the planes apparent in the image, and then
warp the whole image applying a homography to it - for example, you could make the left-side wall appear frontal. But your example has real parallax between the views, so it cannot be achieved with a homography.


Analysis and algorithmic complexity of the Constant Q transform

I didn't verify each exact step of your analysis, but it seems that you're overestimating the memory part by a huge margin. A lot of the kernels are identical except for the exact tiem moemnt at which they're applied. In fact, you only have one kernel per level. So for a typical music application you have only 12 levels per octave, and thus 12 kernels. Over a 10 octave range (20 Hz-20 kHz) that's still only 240 levels. The paper you reference uses quarter-tone resolution (24/octave) and 5 octaves. That's just a design choice.

Next, you state that "With their implementation, they store the kernels in full size". I doubt it. The paper explicitly mentions the effect of trimming the kernels (see Figure 3) and states that the kernels are trimmed to 6 values max (and sometimes just 1 !). This saves both memory and time. The exact amount does depend on the trimming applied.


Measuring Phase Difference at DC bin aliased frequency using FFT

The "DC bin" does not represent DC exclusively, but has also some low frequency energy. How much depends on the fft size and sampling rate. Simple example: if your fft size is 1000 and your sampling rate is 1000Hz, then the dc bin represents frequencies from 0 to 0.5Hz, as the distance between fft bins is 1Hz. So, a phase can be calculated for the DC bin.


Although the question could belong to SE.math, mastering inequalities for $\ell_p$ norms (for $p\ge 1$) or quasinorms (for $0<p< 1$), and their norm ratios and powers, is quite important in  modern signal/image processing: starting from Cauchy-Bunyakovsky-Schwarz, to power means, Rogers-Hölder, Young, Minkowski, Dresher-Gini etc.)  There are scalar, additive and convolutive versions, and reverse forms, that are used in data restoration and deconvolution, especially when using signal sparsity.

If $a_n$ and $b_n$ are positive finite sequences of length $N$, $u,v,w>0$, and $1/u+1/v\le 1/w$, the generalized Rogers-Hölder's inequality gives:

\begin{equation}
\left( \sum_{n=1}^{N} (a_n b_n)^w \right)^{1/w} \le \left( \sum_{n=1}^{N} (a_n )^u \right)^{1/u}
\left( \sum_{n=1}^{N} (b_n )^v \right)^{1/v}\,.
\label{eq_rogers-holder}
\end{equation}

You can find a proof in P. S. Bullen, Handbook of Means and Their Inequalities, Kluwer Academic Publishers, 2nd edition, 2003, p. 187, Corollary 8.
22,

If you now set $a_n=|z_n|^p$, $b_n=1$, $w=p$, $v=q$ and  $1/u=1/p-1/q$, you get
\begin{equation}
\ell_{q}{(|z|)} \leq \ell_{p}{(|z|)}  \leq \ell_{0}^{\frac{1}{p}-\frac{1}{q}}\ell_{q}{(|z|)}   \,.
\end{equation}

$\ell_{0}$ here, not a norm, nor a quasinorm, is the count index, ie the number of nonzero terms in $z$, or the sparsity. In other words, there is  an inclusion in $l_p$ spaces: if $p\le q$, $l_p\subset l_q$.
Your answer is incorrect, because you should have $a\le b$ instead.

One of the basic block is alluded to by  @robert bristow-johnson. If $0<c\le 1$ ($a_n$ positive as above):

$$ \left(\sum a_n\right)^c \le  \sum a_n^c \,.$$

Using additivity, prove it first with two terms:

$$(a_1+a_2)^c \le a_1^c+a_2^c\,.$$

By factoring by the largest of $a_1$ or $a_2$, you end up to prove that for $0<t<1$:

$$(1+t)^c\le 1+t^c\,.$$

Just differentiate the function $t\mapsto 1+t^c- (1+t)^c$ as $c\left(t^{c-1}- (1+t)^{c-1}\right)$, and verify the latter is positive.


Audio from multiple microphones mixed to mono and ASR

Record the audio from the three microphones separately and simultaneously using three audio cards. Use low sensitivity microphone to avoid audio pilferage. You can use ffmpeg command line to record the audio.


A few questions about Slepian and generalized gaussian windows

The Slepian sequences are a family of functions.  Most algorithms calculate 2*NW - 1  sequences at once for a given NW.  N is the number of points in the sequence and W decides half the width of the mainlobe in the frequency domain for the Fourier transform of a given Slepian sequence.  Typically you would use an NW of 3 or 4 for your signal processing.

In scipy, they are asking for the $N$ (m in the Python code) and the $W$ (width in the Python code) parameter separately while in matlab you input the time-bandwidth product $NW$ as a single parameter.  This makes sense because you're typically calculating the Slepian sequences for a fixed window size, $N$.

If you are estimating the power spectrum of a stationary time series the DPSS are the set of windows you should use.

The generalized Gaussian function returns a Gaussian-like function raised to successively higher powers depending on the p parameter.  As p is raised to successively higher powers the generalized Gaussian becomes narrower in the time domain.  The nice property of a Gaussian is that it is it's own Fourier transform and it is the function that achieves the limit concerning the uncertainty principle.  A Gaussian function may be useful if you would like to compute a short-time Fourier transform or spectrogram as an estimator for a time-varying power spectrum of a non-stationary time series.


Yes it is true. Your FFT will compute those DFT values for k=0 to k=99.999 all of which are exactly zero but $X[0] = 7000 \cdot 100000 = 7 \times 10^8 $ for k=0. (the DC)

In order to see why this is the case, we can simply look at the definition of the DFT (as FFT is nothing but a fast algorithm to compute the result of DFT)

The DFT $X[k]$ of a signal $x[n]$ of length N is:

$$ X[k] = \sum_{n=0}^{N-1}{ x[n] e^{-j\frac{2\pi}{N}nk }} ,~~~~~~~~~~~~~~0 \leq n,k \leq N-1 $$

Now assuming the signal $x[n]$ is constant K:
$$x[n] = K ~~~~~~~~~~,~~~~~~ \text{for} ~~~~~~ 0\leq n \leq N-1 $$

The DFT becomes:
$$ X[k] = \sum_{n=0}^{N-1}{ K e^{-j\frac{2\pi}{N}nk }} ,~~~~~~~~~~~~~~0 \leq n,k \leq N-1 $$

There are a number of ways to evaluate or interpret this finite sum of a periodic complex exponential sum all of which will yield the following:
$$X[k]=
\begin{cases}
N \cdot K ~~~~, k = 0 \\
0 ~~~~, 1 \leq k \leq N-1 \\
\end{cases}$$


Generate projection image by OpenGL or Matlab

You would not need OpenGL for such perspective transforms.

First create a 3D cylinder with Coca-Cola texture on it. (This is found even in 3D warehouse).

Then, start by creating a pseudo-camera matrix $\mathbf{K}$ and then generate arbitrary rotation angles with desired range for each axis. Then compose the final rotation by $\mathbf{R}=\mathbf{R}_x\mathbf{R}_y\mathbf{R}_z$. Also synthesize a translation $\mathbf{t}$. Your projection matrix then reads $\mathbf{P}=\mathbf{K}[\mathbf{R} | \mathbf{t}]$. Project your 3D points with this matrix and you will get the warping effect.


If you think of the spectrograms as filter banks it might be easier to conceptualize. With a STFT, each bin in the bank has the same time/frequency resolution. With the Wavelet transform, the bins at lower frequencies have higher frequency resolution, which means they must also have lower time resolution.  There is an inherent trade off between time and frequency resolution such that one cannot be increased without decreasing the other.  This can be thought of intuitively as an impulse in the time domain has an infinitely wide frequency response, or infinitely high resolution in the time domain and zero resolution in the frequency domain.


Sparse Channel Estimation in OFDM System

Are you still there, I've seen this question right now.
First let me correct your question,

Y(k)=H(k)X(k) === > Y(k) received signal in frequency domain

That's true, but it's not necessary that signal in frequency domain in order to be used in compressed sensing, it might be in time domain also.

Now, regarding your questions,
the representation of h channel can be sparse too but it's not always, it can be also in other forms .. taking an example of multi-path environment, the received signal can be sparse, and h channel can be either sparse or not, that's not important.

how to design the dictionary,  that's the goal of the compresses sensing, it's how can you design it,  you can check online, you will get many dictionary.

Estimation is different about using compressed sensing. you may use pilot or not, that's different topic.

reference for compresses sensing, is An Introduction to Compressive Sensing
Collection Editors:
Richard Baraniuk
Mark A. Davenport
Marco F. Duarte
Chinmay Hegde

Thanks


How to convert a spectrogram back to a signal

I am not a good Python coder, but did similar processing in Matlab in the past. The subject has been discussed in SE.DSP in several instances, for instance Librosa stft + istft - Understanding my output.

Tools seem to exist:

scipy.signal.istft: Perform the inverse Short Time Fourier transform (iSTFT).
librosa.core.istft: Inverse short-time Fourier transform (ISTFT).

How you should process complex coefficients is a complicated manner. The most simple consists in modifying the modulus, and keeping the phase. However, the result depends on the type of enhancement you want to perform. You could find inspiration on papers talking about shrinkage (esp. with complex, dual-tree or quaternionic wavelets) and homomorphic filtering. I could update later if you can clarify objectives.


Differences between DSP programming vs MPU programming

first of all, you certainly can use an ARM or Intel or Arduino or whatever.  DSP chips are MPUs that are optimized for the most common DSP operations.

consider an FIR filter.

$$ y[n] = \sum\limits_{i=0}^{N-1} h[i] x[n-i] $$

to perform that in $N$ instructions the DSP has to do, in one instruction, the following:

multiply coefficient $h[i]$ and sample $x[n-i]$.
add the product $h[i] x[n-i]$ into an accumulator.
pre-fetch, for the next instruction both $h[i+1]$ and $x[n-(i+1)]$ and
post-increment the pointers to point to $h[i+2]$ and $x[n-(i+2)]$

if the accumulator is fixed-point, the width of it should be wide enough that there is no rounding on the right and with enough guard bits on the left for overflow.  then there should be a saturate instruction (or the saturation is done automatically when the word is moved to the output).

finally, the pointer that points to $x[n-i]$ must have a circular or modulo-addressing mode so that it "wraps around" since the data in the $x[n]$ is stored in a circular FIFO buffer.

now, a DSP chip will probably do all that with its specialized hardware, but you need not use a DSP as long as you, the programmer, are willing to deal with these details.  given that some MPUs run at a much faster clock rate than most inexpensive DSPs, even if you must do all of these operations with multiple instructions, it still might run faster on, say, an ARM than on a SHArC.  maybe not.

in doing circular addressing, the simplest way to do that is to use a regular array (like in C) but define the size of that array to be a power of 2 (say $2^p$).  then use the AND operation (like "&" in C) to mask off all but the lower $p$ bits.  that will accomplish the modulo arithmetic on the index or address of the data.

other than that, it's mostly a matter of knowing the math for DSP, understanding a little about the numerical issues regarding finite-width words, whether they be fixed or floating point and being careful.  DSP programming is really no different than any other mathematical programming.


Example of Entropy and Channel Capacity Computation

Since all transitions are equally likely, there is no information in the channel output about the source. Hence, the mutual information, and, consequently, the channel capacity, are zero in this case. You might as well look at the output of a random number generator and from that guess the source symbols.

Expressed in formulas you have

$$H(Y)=H(Y|X)\tag{1}$$

because $X$ and $Y$ are independent, and, consequently,

$$I(X,Y)=H(Y)-H(Y|X)=0\tag{2}$$

i.e., there is no information in $Y$ about $X$ (which makes sense because the two random variables are independent).

Since mutual information is symmetrical, i.e., $I(X,Y)=I(Y,X)$ we also have $H(X)-H(X|Y)=0$, i.e.,

$$H(X)=H(X|Y)\tag{3}$$

which means that the uncertainty about the outcome of $X$ is the same whether or not the channel output $Y$ is given.


Why is convolution or FFT breaking in this case?

Yes, zero-padding is the issue, however not as you would except...

@Soren is only partially right: you can add zeros at the beginning, middle or end of the signal just as you did, but absolutely not as you wish.

What's the general rule?

  Look at how the signal repeats with and without zeros and choose the
  place that distorts the signal the least.

What does it mean?

In the case of a real, discrete and aperiodic time-space signal, there is a causal relationship between the value at n and the value at n+1. So all these values should be considered as a bloc, adding zeros in the middle would distort it. Here is an obvious reason from your code, you had a proper Gaussian, not something like this:

In the case of your convolution kernel, it's a bit more complicated. What do the last taps (second half) mean?

In the aperiodic DTFT (i.e. linear convolution, as opposed to circular), these mean that the convolved values depends of both [n, n+1, n+2, ...] (as expected) but also [len - 1, len - 2, len - 3, ...] (where len corresponds to the end of the array without padding). This "classical" linear convolution can be simulated by padding the end/beginning of the signal.
In the periodic DFT (circular convolution), the taps of the second half actually correspond to the values [n-1, n-2, n-3, ...] (not the values by the end of the signal).

np.convolve and bunch try to do a "classical" linear convolution while you seem to want a non-classical one (zeros added in the middle). Here is a comparison:

As I understood, the actual kernel you want is the one described in "your kernel".

You can easily get the right kernel by fftshifting:

kern = np.fft.fftshift(kern)

This has the effect of centering your DC (using linear phase shifting) and then you can just pad zeros at the beginning/end as much as you wish and use the standard convolution/filtering tools.

In the end, this was an interesting question: people usually just pad blindly using convolution functions or options from the FFT transform without ever wondering what they are doing. Some even put fftshifts out of habit every time they deal with anything resembling a kernel.


I think the confusion comes from the fact that the command hilbert in Scipy (and also in Matlab/Octave) does not just compute the Hilbert transform, but its output is the analytic signal. So if $x(t)$ is the (real-valued) input to such a function, its (complex-valued) output is
$$y(t)=x(t)+j\mathcal{H}\{x(t)\}\tag{1}$$
Clearly, if you want to obtain $x(t)$ from $y(t)$, you just need to take its real part.
Concerning the inverse Hilbert transform, as you've already found out, the following holds:
$$\mathcal{H}^{-1}\{x(t)\}=-\mathcal{H}\{x(t)\}\tag{2}$$
The relation $(2)$ holds no matter if $x(t)$ is real-valued or complex-valued. However, if you use the function hilbert on a complex-valued signal, you throw away information that cannot be retrieved, because you throw away the negative frequencies, which are not redundant for complex-valued signals (unlike for real-valued signals). So the function hilbert  cannot be inverted if the signal is complex-valued.


correct for shift from butterworth filter?

The delay or "shift" of an IIR butterworth filter varies with frequency.  If you want a constant delay that you can correct for in a simple manner, you need to instead use a linear phase FIR filter, not an IIR filter.


Volterra Series modeling and its adaptive implementations can be used to perform some simple nonlinear system identification tasks... See if it fits to your problem.


Speech comparison algorithm for rating on similarities

First, you will have to correct for differences in timing. For example, if one utterance is "--heeelloooo---" and the other "hellooooooo----" (- representing silence), a direct pairwise comparison of MFCC frames will show differences simply because the two samples are not aligned. You can use Dynamic Time Warping to find the best alignment between the two sequence of feature vectors - and compute the corresponding distance.

The second problem is that if the two recordings are not from the same speaker, you will have to compensate for the differences in timbre. The MFCC of a female saying "aaa" are not the same as the MFCC of a male speaker saying the same phoneme! A relatively simple model to account for variations in voice timbre is to assume that there exists a linear transform $\Gamma$ which "maps" the MFCCs of one speaker onto the MFCCs of another speaker (to be fair, only one small subset of these transform accurately model how changing parameters such as age, gender, etc "shifts" the MFCC). Once two recordings have been aligned, even roughly, you can use a least square procedure to estimate $\Gamma$. This procedure is known as speaker normalization or speaker adaptation.

Your comparison procedure will thus consists of the following steps. $A$ and $B$ are your original MFCC sequences.

Align the two utterances using DTW. This yields $A'$, a matrix with the observations from $A$ warped/shifted to be aligned with the observations in $B$. You can stop here if $A'$ and $B$ are known to be from the same speaker.
Estimate the transform $\Gamma$ which minimizes the difference between $\Gamma A'$ and $B$.
Use the distance between $\Gamma A'$ and $B$ as your metric.

A last thing that comes to my mind is that you have to discard the first MFCC coefficient (which roughly expresses the signal loudness) to improve the capacity of your system to match utterances pronounced at a different volume/recording level.


How to do base-band simulation of coherent detection of M-ary orthogonal signals?

Since the answer was given by vaz in the comments I figured I would just post it as an actual answer...

To perform coherent detection for M-ary orthogonal signals in a complex base-band simulation simply take the real part of each matched filter (correlator) output and select the symbol with the largest value.

Also, don't make my mistake and use the absolute value of the real part.  I made those changes and verified against the theoretical coherent equations MATLAB uses.


What is used to initialise difference equations that require data from the output?

If there is no good reason to choose otherwise, you would initialize with zeros, so in your case $y[-1]=0$. Other initializations may be useful in certain situations, e.g., when processing blocks of data to avoid transients between blocks.

Note that if the filter is stable (which it has to be in almost all useful applications), the influence of the initial condition on the output becomes negligible after a while.


NCO loses lock after a while for certain values of alpha

I will assume that your PLL is purely digital. Is that correct?

I assume your PLL consists of 3 blocks

The first block calculates the phase of your incoming signal (the signal you wanna lock to).

The second one will compare the phase of the signal to the phase of your "internal oscillator" yielding the error e[n]

The third one is your "controller" which I assume is some kind of PI controller with alpha being the proportionnal gain and beta being the integral gain.

Your phase detection algorithm that detects the phase of your input signal has some kind of low-pass effect. This creates some kind of delay between the output your PI controller and the measurement. This delay puts a limitation on the values of your alpha and beta parameters. I.e. if you increase the values too much, your PLL can be unstable and you will lose the lock.

If you change the phase of your incoming signal, how many samples does it take for the change to be effective? This will tell you approximately what your delay is. From this delay, you should be able to find alpha and beta gains that will give you a good gain and phase margin. You should try at least 6 dB and 45 degrees.

edit : How you obtain the value y*? Is is some kind of delay chain? Hilbert Filter? I assume there is some kind of delay involved...

edit 2 : Your controller seems ok, your gains are low. Your bandwidth is low and thus the closed-loop should be stable. I suspect that your error happens when the phase measured by algorithm transitions between pi and - pi. Picture this case

Phase NCO = (pi - 0.001)

Phase ref = (pi + 0.001)

The phase difference should be 0.002 radian right?

However, phase ref will not actually be (pi + 0.001) but rather (- pi + 0.001)

Yielding error = (-2*pi + 0.002). Which is not the real error, this has the effect of suddenly increasing the bandwidth of your controller, which is something you do not want to do. This is probably why it is stable at lower values of alpha. The glitch would be there but not as big.

I'm assuming your algorithm yields a phase between -pi and pi, if the algorithm yields a phase between 0 and 2*pi, the thinking remains the same.


Rate offset and jitter concealment for streaming audio

at +-2ppm, we'd win or lose a single sample any 45 seconds – my pragmatic approach here would be ignoring the problem, which, of course, has extremely low latency :)

But yes, especially if the producer is a bit slower than the consumer, you'll have a problem in the long run.

Now, there's no easy solution to this – first of all, your receiver needs to be able to see that problem at all.

But let's assume it does – simply by counting (recovered?) clock phases of the receive signal and comparing to its own clock cycles (this gets easier for serial audio data buses, e.g. I2S, because there's many bits per sample, so your bit clock is higher, so the counting offset gets visible much faster, but only if the transmitter continously clocks the bus).
If that can't work for technical reasons (non-continuous bursting of audio data, bus clock is decoupled from ADC clock, no access to recovered clock, etc), the only way I could think of would be the one of keeping track of buffer levels – I assume your receiver does have one, albeit small, buffer between its digital input and the ADC, because otherwise, the ADC would have to run directly off the bus clock, right?

You could observe the buffer "fillage" or clock offset, low pass filter and or caclulate a derivative of that info, define upper and lower boundaries, and simply drop a whole sample if you've got a next-to-critical source-faster-than-sink mismatch, or add a zero sample in the opposite case.

Of course, the "right" way to do this would be resampling – but from your question it seemed to me you're more interested in latency than in accuracy; you use the same frequency offset info as gotten from either bus/clock shift or buffer fillage to drive your own arbitrary resampler:

For example, what's pretty popular among these parts (PC-style computer based SDR), is a resampler based on MMSE-combined filters. They're not overly long.
The next more elegant thing is the Polyphase filterbank arbitrary resampler: again, filterbank with different filters, this time actually used just for a rational interpolator/decimator, but followed by logic that "combines" the output of the next lower rational resampler "branch" with the next higher linearly to get the wanted real resampling ratio. Advantage is that with a big filterbank, your error gets pretty small at reasonable CPU cost, and latency is bound by the length of the rational resampler bank filter - which typically won't be hundreds of samples long :)


Highest frequency of collected force signal

In order to determine the spectral content of your signal you can compute and display its Fourier transform. Assuming you have $N$ samples of the signal $x[n]$ sampled at $F_s = 1500$ Hz, then you can simply determine its frequency content by the following Matlab command:

Xk = fft(x ,N);

Then you can plot its spectrum to see if there is significand energy in the frequency region above $150/2 = 75$ Hz by the following code:

figure, plot ( linspace(0,Fs,N) , abs(Xk) ) ;
xlabel('Frequency [Hz]')
title('DFT Spectrum Magnitude |X[k]| of x[n]');

If there is any signifcand energy above $75$ Hz then there will be significand aliasing. However note that the sensor might have applied a $75$ Hz lowpass filter before sampling the signal, so that you can still avoid alising, but your signal will be distorted, nevertheless, as it would have beeen sampled from a lowpass filtered version of the true signal.


Applying discrete Wigner-Ville distribution to Lena picture?

The developers of Lucid, which is based on time-frequency analysis, report interesting results. Technical details are provided in this PowerPoint file.


How to detect rectangular images inserted into another image?

In the two positive examples, the rectangular box is "likely" to be noticed visually, notably by the frame, but also the difference of the content inside and outside the box, especially in the color diversity.
However, the compression (if lossy, JPEG_like) may induce local fake edges (quantization of the DC and higher frequency), local noise and fake colors along $8\times 8$ squares.
Hence, I would suggest a local/global algorithms in three steps.

First, estimation: a crude segmentation with closed contours. For instance with geodesic active regions or level sets. The principle is as follows: from a starting contour, the border evolve until it meet a region that either has a gradient (edge) and/or the statistics of the regions on both sides are different enough. Generally they are performered on colored gradient-tensors. Although the contour has some regularity, it is allowed to split and separate into several region. I think the is doable by maximizing the color contrast and allowing sharp edge
Second: model matching. You get a set of connected contours, which can be fitted by Cartesian rectangular boxes. The latter can be probed with your size constraints. One question remains: how to select only one box?
Third: decision/machine learning. Given a set of images, and a training ensemble of negatives and positives examples, you can train a logistic classifier on features related to the inner/outer boxes.
If you cannot learn examples, a distance between the feature-vector of the inner box to the feature-vector of the outer box could provide a quantitative metric related to some probability of having a rectangular image  inserted.


Multiplication of two signals is formally called as modulation; one of them is modulating the other.

In image processing most typically this will be called masking, when one of the images is a special selective mask for pixels.

However, when neither of the images is a mask, then their multiplication will blend (mix) them into each other according to their bright and dark regions...


Characterising a system from its input and output

Use the matlab cpsd (cross power spectral density) function on your two signals.
That will give you the frequency response.
Do an xcorr (cross correlation) to find the time delay between the two signals.
Do this for several pairs of input and output files and see if the results match between different pairs.
If they are (reasonably) similar then you could create an FIR filter that reproduces the frequency response and delay.

You could then apply the filter to a set of input data to approximate the expected output.

If there is a PID controller in the middle, however,  I don't know if this procedure will be enough.


If we FFT a signal with number of points far less than the length of the signal, the result is unreliable?

As per the Matlab help, y = fft(x, 1024) where length(x) > 1024 will result in just taking the FFT of x(1:1024).

So if talking the FFT of just the first 1024 samples of your 1,000,000 point signal is OK, then sure, what you say is correct.


Unsupervised Way to Choose K in K-Means Algorithm

The idea underlying the k-means algorithm is to find clusters that minimize the intra-cluster variance and maximize the inter-cluster sum-of-square. Note that the total variance is fixed so minimization of intra-cluster variance is equivalent to maximization of inter-cluster sum of square. Basically various methods aim to achieve such minimization.

However, there is no "correct" K number or "best" method in k means. One clustering may be better than another by one metric, but worse by another metric. And depending on your own project, sometimes one clustering works but sometimes it does not.

You had better try those method to see whether they are suit in your case. If your data is not high-dimensional, you can just try kernel density estimation to observe the number of peaks. Otherwise, the easiest and heuristic one is rule of thumbs, then I would suggest Elbow method (cease increasing the number of clusters when the new generated cluster is closed to some of the existing).

There are many discussions on stackoverflow and cross validation on this topic that you may be interested in as well.


Confusion regarding model order and lags

Hi: As far as what you wrote in confusion 2):

The AR expression should not have an $x[t]$ on both sides. For example, below is an AR(2):

$x[t] = a1*x[t-1] + a2*x[t-2]+ e[t]$

Similarly, the MA expression should always have a coefficient of 1 on the  $e[t]$ term. So, for example, an MA(2) would be:

$x[t] = e[t] + a1*e[t-1] + a2*e[t-2]$

As Stanley said, it is only a matter of convention but the order of an AR(p) is generally thought of as $p$ so there are $p$ AR coefficients if the order is p. Similarly, for an MA(q) there are q MA coefficients so the order is considered to be q ( the $e[t]$ is not included in the order for an MA because there is no coefficient associated with $e[t]$).

Note that, for both MA and AR models,  the number of parameters to be estimated is not always equal to the order of the model because there may be non zero mean terms or trend terms. For example, below is AR(2) with unknown mean $\mu$:

$x[t] = \mu + a1*x[t-1] + a2*x[t-2]+ e[t]$.

Note that there are also ARMA models where the number of coefficients is then p+q. For example, an ARMA(2,1) is:

$x[t] = a1*x[t-1] + a2*x[t-2]+ b1*e_[t-1] + e[t]$.

But usually the term "order" is not used for ARMA models.

I hope this makes it clear.


Detect lines transitioning to blobs in images

If you are not interested in picking out the frame that bubbles start forming and since the transition is (clearly) from order (clearly formed shapes) to disorder (mostly noise), you can use a feature based approach. That is, characterise the whole image with a few numbers.
An obvious feature here is image entropy. The typical (but not the only) way to calculate entropy of an image / signal is through its normalised histogram (you need to normalise it, to get to a probability).
Are there any differences between the histograms of the images?
Yes there are (at least, the ones that come with the question):

From left to right, we have the image with the well defined shapes (ordered), the middle state and on the right the final state where there is mostly background plus some noise.
The striking feature here is that when the image is ordered, the ridges that are well defined because they create contrast with their background create a bimodal histogram. One mode corresponds to the background (the lower one) and one mode corresponds to the "foreground" (the one that forms a "shelf" towards the right) or the worm-like ridges.
You can catch these differences in many different ways including (of course) distance from Gaussian distribution, skewness, etc. Where the image is going to end up is a set of pixels distributed normally, any deviation from that "profile" must be due to formation of those ridges (in this particular context).
If you try to measure (Shannon) entropy over these images, you get the following three values:
    Ordered: 7.2331079063161425
    Middle:  6.734680695379271
    Noise:   6.591988661960253

(These measurements come from skimage.measure.shannon_entropy() )
This number looks at the frequencies of occurence of the pixels and returns the number of bits required to describe the value of a pixel.
You might feel here that these numbers are contradicting: Low entropy means more order, high entropy means less order. "Less order" here means closer to an absolutely uniform distribution. That is, you cannot really predict what comes next, the image is entirely random.
But at the same time, this description might also help understanding what is going on here: The "Ordered" image alternates between two means (more preference on the lower ofcourse) and spans a wider range of highly probable grey tones and the "Noisy" image is basically $128 \pm \approx 70$ with diminishing probability.
Shannon entropy evaluated in this way looks at the pixel values independently of each other. But they are not. The meandering structure forms straight and curved points to the extent that if one pixel is "on", it is highly likely that the pixel next to it is on as well. So, we are not talking about single pixel probablities any more.
To "catch" this kind of differences, you need a measure of entropy that is sensitive to them. One of them is Entropy as calculated through the co-occurence matrix.
The name is kind of self-descriptive but very briefly: The co-occurence matrix describes how the values of two different pixels that are some distance appart co-vary. Additionally, since we now talk about the distance between two pixels, we can also define a bearing. This offset, implies that the co-occurence matrix describes co-occurence along a direction.
The corresponding values now are:
    Ordered: 14.197350713914734
    Middle : 13.128723041575485
    Noisy  : 12.824963198002461

(These measurements come from mahotas.features.haralick() with return_mean=True. Otherwise, you get the value of each feature for each direction considered (typically, the 4 directions around a pixel))
We observe a similar trend.
The difference between the two ways of computing the entropy here is robustness and sensitivity. We do not have a long enough sample here, but it is likely that the haralick features would be more sensitive to structural changes than a single histogram.
Therefore, you can get a long enough set of images from the experiments, "map" how the value of this feature varies and then set up "alarms" by setting thresholds around specific ranges.
Hope this helps.


How to report how a filter was designed (scientific papers) + how to calculate filter attenuation?

Five of those (low/high-pass, cutoff freq, and IIR/FIR, how the filter was designed, attenuation) can be easily described using words or numbers.  The others (phase response, transition bands) are often better described with a picture.

I would recommend describing those five "describable" characteristics and, if appropriate, put pictures of the phase and magnitude response along with the actual filter taps in an appendix or something.


Is it possible to open up this particular convolution equation?

Here is a re-write of a previously mistaken answer:

Let us write $y[n] = z[n] * \Big(x[n] \cdot w[n] \Big) = \displaystyle\sum_{m=-M}^{M} {z[n-m] \cdot (x[m] \cdot w[m])}$ and rearrange as $y[n] = \displaystyle\sum_{m=-M}^{M} {(z[n-m] \cdot x[m]) \cdot (w[m])}$. After this point, we cannot take the terms out of the brackets.


A counter-example would be an ideal LC-circuit, which is passive but not BIBO stable: if excited at its resonance frequency, the output grows without bounds. Such systems are sometimes called marginally stable, because they have poles on the stability boundary (the imaginary axis).
Note that if a system is not lossless, i.e., if there is a positive resistance such as in an RLC-circuit, the system is strictly passive, and as a consequence it is also BIBO-stable. I.e., strict passivity implies BIBO-stability (all transfer function poles are in the left half-plane).


Because DMRS is only sent if there is data to transmit and, therefore, narrowband while SRS, as the term sounding suggests, is not associated to any data transmission and also wideband.
For example, SRS can be used by base stations to choose the best portion of the UE channel bandwidth to schedule next PUSCH.


How to determine the parameter of a Complementary Filter?

This is the difference equation for a low pass filter.  Five time constants (5 * 1 time constant) is the time it will take to for the output, to reach 99.33% of the value of the input, from when the input changes from 0 to its final value, and stays there (a step response).  This comes from 1-e^(-5).  As you are specifying time constants, it implies you want a time domain definition for use as a moving average filter, so specify the period of averaging you would like, assume that as 5 time constants and divide by 5 five to get one time constant for use in the equation you have quoted above.


Raised-Cosine Filter - calculating phase response

Just like the sinc pulse, a raised cosine pulse with zero delay (centered around $t=0$) has a real Fourier transform. You can find its expression in section 11.4 of "Telecommunications Breakdown" by Johnson and Sethares (a free draft version is available here).

For a time-shifted pulse, you can apply the time-shift property of the Fourier Transform.


How to check if the poles of transfer function are stable?

It's not $|z_{\inf,i}<1|$ it's $|z_{\inf,i}|<1 $ which means that, as @Mohammad Mohammadi said in the comments, your pole on a complex plane must be inside the unit circle.

Analytically, you get the magnitude with
$$ r = |z| = \sqrt{a^2+b^2} $$
If $r<1$, your pole is stable.

Note that $z_1 = b + ja$ and $z_2 = b-ja$ have the same magnitude, which means that if one is stable, the other is stable as well (in the complex plane, they are symmetrical to the horizontal axis).

Illustration

Matlab code for the illustration

%%
clear
close all
clc;

% STABLE

a = 0.73;
b = 0.24;

P1 = a + b*1i; % Pole 1
P2 = a - b*1i; % Pole 2

r = sqrt(a^2+b^2); % Magnitude

% UNSTABLE

a2 = 1.23;
b2 = 0.84;

r2 = sqrt(a2^2+b2^2);

P12 = a2 + b2*1i; % Pole 1
P22 = a2 - b2*1i; % Pole 2

zplane([0 0; 0 0], [P1 P12; P2 P22]); % Plot the poles on a complex plane
legend('', '', ...
    ['z = ' num2str(a) ' \pm ' num2str(b) 'i (r = ' num2str(r) ')'],...
    ['z = ' num2str(a2) ' \pm ' num2str(b2) 'i (r = ' num2str(r2) ')'], ...
    'Unit Circle (r=1)')


What does this notation $g(x) = S\{\delta(t) \}(x)$ mean?

My understanding is that a system $S$ is a notation for   a generic object, which takes instances depending on its arguments, given between $\{\cdot\}$. Here, $S\{\delta(t)\}$ instantiates it as its effect on the zero-centered impulse response, with free variable $t$, taken on variable $x$.

Now, if this system $S$ is linear and shift-invariant, the operator given by $S\{\delta(t)\}$ defines $S$ it completely.

If the system applies on discretized inputs (that we can call $x$), it can be rewritten as a linear operator, represented by an infinite Toeplitz matrix. Each row implements the (linear) convolution, and the Toeplitz structure emulates the shift-invariance. I suggest to read Toeplitz and circulant matrices by R. Gray. The $0$ in $A_{x,0}$ may refer to the way an infinite matrix is matched with a signal centered at $0$, but I still wonder if this notation is clear.


Transform matrix creation in compressed sensing

I recommend not forming the transform matrix directly. That will cost more memory and take longer to compute the transform. Instead, use the idwt/dwt or waverec/wavedec functions in place of multiplying by A/A.T.
Many reconstruction solvers and convex optimisation toolboxes support supplying function handles for this instead of a matrix. If you are implementing your own solver it should also be easy enough to integrate these functions instead of the matrix-vector products.


What's the difference between SISO and SIMO in using matched Filters

I don't know why we can't use RRC filter with OFDM ?!! Yes, of course, we do can. The only disadvantage for that is the complexity. however the gain will be in spectral efficiency improvement and bandwidth efficiency because of avoiding the needs of guard interval which is used in the regular OFDM. In Engineering, we should always ask, why and why not? why should we use it and why can't we use it?

So, Firstly, it's important to know that ISI in context of OFDM is different about ISI in the other types of modulations. As known,  the ISI happens when the symbol, in time domain, has a tail to the time of next symbol, but  that's not the same case for OFDM. In term of OFDM, that issue can be solves in two steps, because of the orthogonality, (and remember the use of FFT and iFFT), and second we insert the guard interval like CP which create circular convolution to eliminate the ISI. Be careful, in multi-path channels and large delay spread, that ISI will be exist in OFDM !! so we should present new solutions to deal with. "That exactly what the example reference,you provided above, tries to deal with that issue too".

But the question you will ask now, Why regular OFDM doesn't use RRC filter?  In reality, As I mentioned above, Regular OFDM is using Rectangular filter instead of RRC filter for minimal complexity, but the cost of that low complexity is a.) inserting CP which reduce the bandwidth efficiency b.) spectral efficiency reduction by sidelobes.

These two issues in regular OFDM, CP insertion and low spectral efficiency, can be solved by using RRC-OFDM, this RRC-OFDM doesn't use CP and the spectral efficiency is improved very well there. That technique was first presented in 2001, and then patented as US7065150 in 2006.

Coming back to the question, Should we use one matched filter in received SIMO signal or Multi filters. That depends on the simulation you are doing. you can use only one, but be careful how implement it. Or try to filter each copy of signal separately and then collect them back again and continue your work.

thnx


First of all - don't use 80kHz (I don't know why you came up with this number as a limit for MATLAB, whereas it is 1e6 in fact) but something that is more reasonable and standardized, i.e. 96kHz (or 88.2kHz if you wish...). Another thing is that your sound card might not be able to serve such high sampling frequency. You should check that in playback options --> Advanced tab and set it correctly. Last thing - you should drop using obsolete sound and use audioplayer - unless you use some ooold MATLAB.


How many bps would be required to reproduce the words in speech recognizably?

This all depends how you interpret the requirements and what is "good enough".

For example, a really bandwidth efficient  method would to simply encode with a speech to text system, send the text and re-encode with a text to speech system. That would be about 0.075 kb/s (150 words per minute, 6 letters per word, 5 bits per letter).

That would make a pretty terrible speech coding system but whether that meets the "rules" or not is more a semantics than a technical question.


Deriviation of the "Twiddle Sum" property

The most straightforward way to see this is to note that for $k=mN$

$$W_N^{kn}=e^{j2\pi mnN/N}=e^{j2\pi mn}=1$$

So the sum for the case $k=mN$ is simply

$$\sum_{n=0}^{N-1}1=1+1+\ldots+1=N$$

Note that the solution using L'Hopital's rule is a bit dubious because for $k=mN$ the formula for the geometric series is not valid because the terms are all equal to one (note the $r\neq 1$ condition in the link above). So you arrive at an expression that is only valid for $k\neq mN$, and then you somehow massage it into a solution that also works for $k=mN$, but this is bad math, even for engineers.


Audio signal processing book suggestions

If you can get it:

Zölzer, Udo (1997). Digital Audio Signal Processing. John Wiley and Sons. ISBN 0-471-97226-6

or

Orfanidis

might have to buy them used.

is your interest music processing and synthesis?

there are some nice online books.  Look for anything Julius Smith.  But there are others, but i can't remember any names to search for.

Remembered one:

Davide Rocchesso

here's another link

another old one:

Ken Steiglitz


It is an easy concept but a sometimes it is hard to express it compactly.
The difficulty arises because of additional concepts to be known for any formal description.

As you may know, mathematical analysis of noise variables and signals are handled via the tools of probability theory, in particular by the stochastic (random) processes.

The usual progress is such that first properties of random variables are described and then properties of random signals are derived based on it.

The mean for a random variable, X, is defined to be its Expected Value E(X) and computed based on its Probability Density Function $f_X(x)$ such as $\mu_X = \int {xf_X(x) dx}$ with the limits from minus to plus infinity for a continuous random variable.

So a zero mean random variable is that one for which the above integral is zero. (for a discrete random variable the integral is replaced by a sum and pdf is replaced by a pmf). This is also called as the average value.

Random signals are considered to be made up of an infinite set of random variables each making a single value at a single instance of observation of the instance being created. Hence there is the concept of ensemble and time averages.

Now coming to random signals X(t) (noise) expected value of a random signal is also expressed as E(X(t)) which for a Stationary (or at east Weakly Stationary up to second order) Process is a fixed value.

Your answer is that a zero mean noise is that one for which E(X(t)) is zero for all t.

However above is a theoretical description of Mean. In practice to compute mean and other theoretical parameters, you will refer to statistical computation which deal with practically observed data and computations based on that.

The key which connects statistically made (practical) computations to theoretical computations is the Ergodicity Theorem which states that for an ergodic random process theoretical computations are equal to statistical ones.


Generate the Matrix Form of 2D Convolution Kernel

I created a function to create a Matrix for Image Filtering (Similar ideas to MATLAB's imfilter()):

function [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% Generates an Image Filtering Matrix for the 2D Kernel (The Matrix mH)
% with support for different operations modes (Convolution / Correlation)
% and boundary conditions (Zeros, Symmetric, Replicate, Circular). The
% function should match the use of MATLAB's 'imfilter()' with the same
% parameters.
% Input:
%   - mH                -   Input 2D Convolution Kernel.
%                           Assumed to have odd dimensions.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numRows           -   Number of Rows.
%                           Number of rows in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - numCols           -   Number of Columns.
%                           Number of columns in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - operationMode     -   Operation Mode.
%                           Sets whether to use Convolution or Correlation
%                           for the operation mode.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2}.
%   - boundaryMode      -   Boundary Condition Mode.
%                           Sets the boundary condition mode for the
%                           filtering. The options are - Zeros, Symmetric,
%                           Replicate and Circular.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, 4}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output filtering matrix. Multiplying in
%                           the column stack form on an image should be
%                           equivalent to applying 'imfilter()' on the
%                           image.
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx2()' - https://www.mathworks.com/help/images/ref/convmtx2.html.
% Remarks:
%   1.  The height and width of 'mH' are assumed to be odd number. In case
%       either or both are even the user should pad the kernel with zeros
%       (Either a row, column or both). according to the anchor of the
%       kernel the user do the padding pre or post the kernel.
% TODO:
%   1.  Refactor the code to share the common operations of different
%       boundary modes.
% Release Notes:
%   -   1.0.001     30/12/2019  Royi Avital
%       *   Fixed some bugs related to using 'numCols' instead of 'numRows'
%           in the calculation of 'pixelShift' for the cases 'jj + ll >
%           numCols' and 'jj + ll < 1'.
%   -   1.0.000     16/01/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

OPERATION_MODE_CONVOLUTION = 1;
OPERATION_MODE_CORRELATION = 2;

BOUNDARY_MODE_ZEROS         = 1;
BOUNDARY_MODE_SYMMETRIC     = 2;
BOUNDARY_MODE_REPLICATE     = 3;
BOUNDARY_MODE_CIRCULAR      = 4;

switch(operationMode)
    case(OPERATION_MODE_CONVOLUTION)
        mH = mH(end:-1:1, end:-1:1);
    case(OPERATION_MODE_CORRELATION)
        % mH = mH; %<! Default Code is correlation
end

switch(boundaryMode)
    case(BOUNDARY_MODE_ZEROS)
        mK = CreateConvMtxZeros(mH, numRows, numCols);
    case(BOUNDARY_MODE_SYMMETRIC)
        mK = CreateConvMtxSymmetric(mH, numRows, numCols);
    case(BOUNDARY_MODE_REPLICATE)
        mK = CreateConvMtxReplicate(mH, numRows, numCols);
    case(BOUNDARY_MODE_CIRCULAR)
        mK = CreateConvMtxCircular(mH, numRows, numCols);
end

end

function [ mK ] = CreateConvMtxZeros( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if((ii + kk <= numRows) && (ii + kk >= 1) && (jj + ll <= numCols) && (jj + ll >= 1))
                    vCols(elmntIdx) = pxIdx + pxShift;
                    vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);
                else
                    vCols(elmntIdx) = pxIdx;
                    vVals(elmntIdx) = 0; % See the accumulation property of 'sparse()'.
                end
            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxSymmetric( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (2 * (ii + kk - numRows) - 1);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (2 * (1 -(ii + kk)) - 1);
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((2 * (jj + ll - numCols) - 1) * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((2 * (1 - (jj + ll)) - 1) * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxReplicate( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (ii + kk - numRows);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (1 - (ii + kk));
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((jj + ll - numCols) * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((1 - (jj + ll)) * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxCircular( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - numRows;
                end

                if(ii + kk < 1)
                    pxShift = pxShift + numRows;
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - (numCols * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + (numCols * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

The code was validated against MATLAB imfilter().

Full code is available in my StackOverflow Q2080835 GitHub Repository.


Is the definition of recursion and iteration in signal processing different from computer science?

In my understanding, the term iteration in DSP is usually used in the mathematical sense, i.e., for describing iterative methods solving certain numerical problems. The power method for computing the largest eigenvalue of a matrix and its corresponding eigenvector is a typical example of such an iterative method used in signal processing.

The term recursion is normally used in the sense of recurrence relation, where a value $y_n$ is generally a function of $n$ and of its values $y_k$, $k<n$:

$$y_n=f(n,y_{n-1},y_{n-2},\ldots)$$

In signal processing applications, such as filtering, the index $n$ is interpreted as a time index. E.g., when implementing a (causal) filter, the current filter output can be computed from the current input and previous input and output values. This is a specific form of a recurrence relation. On the other hand, when talking about iteration, the iteration index is usually not related to time. I think it is this latter point that your professor referred to.


this is an old resource, but good.   in fact, i have a very old contribution (another 3 decade old design) that i will repeat here at SE:

A method that Orfanidis (Introduction to Signal Processing mentions came from an old comp.dsp post of mine. (here's a pdf, check Problem B.9)  It's just a simple "pinking" filter to be applied to white noise.  since the rolloff is -3 dB/octave, -6 dB/octave (1st order pole) is too steep and 0 dB/octave is too shallow.

An equiripple approximation to the ideal pinking filter can be realized by
alternating real poles with real zeros.  A simple 3rd order solution that i
obtained is:

$$ H(z) = \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

where

$p_1$ =  0.99572754

$p_2$ =  0.94790649

$p_3$ =  0.53567505

$q_1$ =  0.98443604

$q_2$ =  0.83392334

$q_3$ =  0.07568359

The response follows the ideal -3 dB/octave curve to within $\pm$0.3 dB error over a 10 octave range from 0.0009$\times \pi$ to 0.9$\times \pi$.  ("Nyquist" = $\pi$.)


I don't have the privilege to comment, so I'm going write my comment as an answer.

I'm guessing your function, as well as imfill in Matlab, performs the "filling" iteratively.

Two basic improvements to speed these kind of iterative functions up would be to 1. increase the seed points, i.e. using other possible locations in addition to the border area, maybe base on image intensity before the binerization process, or whatever characteristics the images have, and 2. try to parallelize your function (see Dirk-Jan Kroon's Region Growing function for an example).  If this doesn't speed things up, or you just have too many seed points, then I would suggest you start implementing a function that uses your GPU for parallel processing.

Now there may be "non-iterative" functions that fill 3d holes... so I'll leave that to the experts here on dsp.se.  It would also be nice if you can show us a slice of your 3D image for a better understanding of what you are dealing with.


The function varies as you increase the time step ∆t , t will be zero from the start, So at ∆t =0, the initial condition L*(0) is the gray scale value obtained by converting the color image into L* a* b*  color space.


How to use low pass filter in Matlab for discrete data?

Check conv function. Use conv(f,x) where f is the filter and x is the signal.

You can also implement the filter process in frequency domain:

X=fft(x);
F=fft(f);
X_LP=real(ifft(X.*F));

You can use filter too. In y = filter(b,a,x), the filter is described by numerator coefficient vector b and denominator coefficient vector a,which is different in conv where the input is the time series of the filter. But conv(f,x) and filter(f,1,x) gives the same result in your case.

EDIT

For dfilt.dffir type filter, use the code below (example page):

x = randn(100,1);    % Original signal.
b = fir1(50,.4);     % 50th-order linear-phase FIR filter.
hd = dfilt.dffir(b);    % Direct-form FIR implementation.
y = filter(hd,x);

filter can handle FIR and IIR systems, while conv only works for FIR filter.

Besides,  filter can also return the filter states, so that it can be used in subsequent calls without incurring filter transients.


How is cross-correlation related with orthogonality?

I suppose you mean the cross-correlation at lag zero. Well take an Hilbert space $H$ (i.e. a metric space in which you can define a scalar product $\langle\cdot ,\cdot\rangle$). Then $x,y\in H$ are orthogonal if $\langle x,y\rangle=0$, by definition.
If your Hilbert Space is $L_2(\mathbb{R})$ (the space of real square integrable functions) then the scalar product of $f,g\in L_2(\mathbb{R})$ is defined as
$$
\langle f, g\rangle = \int_{-\infty}^\infty \overline{f(t)}g(t) dt
$$
If the signals are assumed to be ergodic, which is actually a costumary assumption, then the cross-correlation is simply a sliding scalar product, you have
$$
\text{corr}(f,g)(\tau) = \int_{-\infty}^\infty\overline{f(t)}g(t+\tau) dt
$$
hence
$$
\text{corr}(f,g)(0) = \langle f ,g \rangle
$$
and uncorrelation means orthogonality.
The same applies for discrete time signals in $l_2$


Comparing 2 sinusoids: Discrete case vs Continuous case

Because of the $2\pi$-periodicity of the cosine and the discrete nature of $n$ ($n,k\in \mathbb{N}$):
$$
x[n] = \cos[\Omega_\mathrm{x}] = \cos[(\Omega_\mathrm{y}+2\pi k)n] = \cos[\Omega_\mathrm{y}n+2\pi kn] = \cos[\Omega_\mathrm{y}n] = y[n]
$$
This important "phenomenon" is called aliasing and only occurs for discrete-time signals. The normalized frequency $\Omega$ is therefore usually restricted to the interval $[-\pi, \pi]$.


There are a number of ways to do this, but I like this way:

// very rough pseudocode
step_size = 9000000 / width_of_window;  // 9000000 or however many points
for each pixel in width of window
  figure out what portion of the audio stream this x value corresponds to
  take step_size shorts corresponding to position in audio stream
  find the minimum and maximum values in that portion of the stream
  draw a vertical line at the current x value with position and length corresponding
      to the min and max values you found

So basically you're displaying the min and max values of successive portions of the stream in each column of your draw box, or canvas, or whatever you're using.


Classifying a time series with a Lorentzian power spectral density

While it doesn't exactly give the power spectral density formulated by you, a bimodal (when considering also negative frequencies) power spectral density is given by the time evolution of the position or displacement of a harmonic oscillator driven by Brownian motion. The only modification compared to the Ornstein–Uhlenbeck process is inclusion of inertia, which enables oscillations when the harmonic oscillator is underdamped. An Ornstein–Uhlenbeck process is governed by the differential equation:

$$\dot x(t) = -\theta x(t) + F(t) \tag{1}$$

The dot above the $x(t)$ indicates the derivative of $x(t)$. In an infinitesimal time step $dt$, displacement $x(t)$ relaxes a step $\theta\,dt\,x(t)$ towards zero and is also affected by white noise external forcing $F(t)$. There is no inertia: neither the relaxation term or the external forcing term describe "memory of past velocity". I have poor differential equation handling skills, but my hunch is that if a bit of the past velocity $\dot x(t-dt)$ was mixed into the current velocity, it would give this modified form with inertia:

$$\ddot x(t) = -\gamma_m \dot x(t) - f_m^2 x(t)+F(t) \tag{2}$$

Here $\ddot x(t)$ denotes the second derivative of $x(t)$. Now there are terms both for damping of velocity and for relaxation of the displacement. For the modified process, the power spectral density is of form (Gröblacher 2012, Coffey & Kalmykov 2012):

$$S(f) = \frac{a}{\left(f_m^2-f^2\right)^2 + f^2\gamma_m^2}, \tag{3}$$

where $\gamma_m$ is the width of the resonant peak at half maximum, and $a$ is a normalization constant which has also absorbed the power spectral density of the driving force which is white (constant in frequency) for a Brownian heat bath. The true peak frequency $f_\text{res}$ is (Gröblacher 2012):

$$f_\text{res}=f_m\sqrt{1-\frac{\gamma_m^2}{2f_m^2}}.$$

Figure 1. Power spectral density of the displacement of a massive (as opposed to massless) harmonic oscillator in a heat bath, with $a = 1$, $f_m = 1$ and $\gamma_m=\frac{1}{2}$.

$S(f)$ of Eq. 3 is identical in form to the square of the magnitude frequency response of a two-pole resistor-inductor-capacitor (RLC) filter:

$$|H(if)|^2 = \frac{a}{(f^2 - f_0^2)^2 + f^2\,(f_0/Q)^2}$$

If the filter is fed white noise, the power spectrum of the output will be identical to $S(f)$. It is characterized by a -12 dB/octave (exactly -40 dB/decade) asymptotic slope:

Figure 2. Square of the magnitude frequency response of the two-pole filter, with logarithmic axes. The constants are equivalent to those in Fig. 1.

I'll return to the "just add inertia" argument once more. Consider a discretization of the Ornstein–Uhlenbeck process defined in Eq. 1:

$$\dot x[k] = c_0 x[k] + F[k]$$

We can define velocity and acceleration as:

$$\dot x[k] = x[k] - x[k-1]$$
$$\ddot x[k] = \dot x[k] - \dot x[k-1] = x[k] - 2 x[k-1] + x[k-2]$$

If we add damped inertia to the process as a term $c_1 \dot x[k-1]$:

$$\begin{array}{rrcl}
&\dot x[k] &=& c_0 x[k] + c_1 \dot x[k-1] + F[k]\\
\Rightarrow& x[k] - x[k-1] &=& c_0 x[k] + c_1 x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& 0 &=& (c_0 - 1) x[k] + (1 + c_1) x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] + (1 - c_1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] - (c_1 - 1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 - 1) \dot x[k] + c_0 x[k] + F[k]\\
\Rightarrow& \ddot x[k] &=& \frac{c_1 - 1}{c_1} \dot x[k] + \frac{c_0}{c_1} x[k] + \frac{1}{c_1}F[k]\\
\end{array}$$

We have arrived at the discrete equivalent of Eq. 2.

References

Gröblacer, Simon (2012) "Quantum Opto-Mechanics with Micromirrors: Combining Nano-Mechanics with Quantum Optics" DOI: 10.1007/978-3-642-34955-3, chapter 2, page 7.

Coffey, William T.; Kalmykov, Yuri P. (2012) "The Langevin Equation: With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering" 3rd ed. ISBN: 978-981-4355-66-8, page 241.


"Image processing method" here really boils down to ... "for" loops.
You'd walk the main diagonal (or the first non-main diagonal next to it), and as soon as you hit a value, you count the values right of that aren't 0. Soon as you know how many these are, you're done. Half that number, add it to the number of the row, and tadah, center of square.
To deal with overlap is trivial: once you've reached the end of a row go one element further and then you go up until you hit either the diagonal or another non-zero value.
Note that you say "sparse" matrix: depending on the matrix storage format, the "how many entries are set in this row" might already be how the whole thing is stored, making this even easier.


Finding the Best Gaussian Smoothing Kernel to Minimize the Discrepancy Between Two Images

This is closely related to Blind Deconvolution.
The only difference is we limit our self to a very specific type of blur kernels.

The nice thing about the Gaussian Kernel is being defined by single parameter - The Standard Deviation of the kernel.
The less nice thing is the connection isn't linear.

Optimization Problem

Let's define a classic non linear model for this problem:

$$ \arg \min_{\sigma} \frac{1}{2} {\left\| A \left( \sigma \right) x - b \right\|}_{2}^{2} $$

Where $ A \left( \sigma \right) $ is the convolution matrix generated by a Gaussian Kernel parameterized by $ \sigma $, $ x $ is the original image ($ A $ in your question) in a vector shape, and $ b $ is the blurred image ($ B $ in your question).
This is a classical Non Linear Least Squares problem which can be solved by MATLAB using lsqnonlin().

Code Sample

This is the main part of the code:

mA = rand([numRows, numCols]);

gaussianKernelStd = kernelStdLowerBound + ((kernelStdUpperBound - kernelStdLowerBound) * rand(1));
mB = ApplyGaussianBlur(mA, gaussianKernelStd, STD_TO_RADIUS_FACTOR);

% Objective Functions
hObjFun = @(kernelStd) reshape(ApplyGaussianBlur(mA, kernelStd, STD_TO_RADIUS_FACTOR) - mB, [numPx, 1]);
estKernelEst    = lsqnonlin(hObjFun, initKernelStd, kernelStdLowerBound, kernelStdUpperBound, sSolverOptions);

Results

As can be seen, the estimation is almost perfect.
On large images it might take time (Using some tricks of the Gaussian Filter in the Fourier Domain one could do that there with major speed up), but still, it is not free.

The full code is available on my StackExchange Signal Processing Q49121 GitHub Repository.


Why Isn't the ML Estimator (MLE) in MIMO Spatial Multiplexing Obtained by the Least Squares Solution?

You must attention to the written text.
It doesn't say the ML isn't optimal, what it says is that the problem isn't regular LS problem but Least squares problem with Constraints.

The constraints make analytic solution infeasible and hence in order to solve it one must go through any signal in the space of valid signals and mark the one which minimizes the LS problem.

What they mean that the regular LS solution result doesn't necessarily yield a solution in the space of valid space of signals and even taking the closes one to the solution isn't optimal.

So we have the problem the constrained problem:

$$\begin{align*}
\arg \min_{s} \quad & \frac{1}{2} {\left\| H s - y \right\|}_{2}^{2} \\
\text{subject to} \quad & s \in {QPSK}^{M}
\end{align*}$$

Which its solution is optimal.

What they say isn't optimal is solving:

$$\begin{align*}
\arg \min_{s} \quad & \frac{1}{2} {\left\| H s - y \right\|}_{2}^{2} \\
\end{align*}$$

And then find $ s \in {QPSK}^{M} $ as you'd do in MIMO.
Namely, the greedy method (Do LS while ignoring the constraint and then apply the constraint) isn't optimal.

Yet the problem is that there is no analytic or efficient way to solve the Constraint LS problem but doing brute force search which isn't feasible.

I agree the text isn't clear about making this observation.
Yet it well known in the optimization field that solving in this greedy 2 step method isn't guaranteed to be optimal.

Remark
I found a link to the PDF - Doron Ezri - MIMO OFDM Lecture Notes. What I talked about can be seen in part 5.4.


What image processing operation should I use to remove fuzzy black area

You can perform local averaging/blurring of the picture and then use a threshold to remove all areas that are too dark:
img = cv2.imread('k1aqp.png')
figure(figsize=(20,20))
subplot(2,2,1);
imshow(img)

subplot(2,2,2);
B = cv2.blur(img, (55,55), 0)
imshow(B)

subplot(2,2,3)
mask = B > 150
imshow(mask)

subplot(2,2,4)
img2 = img.copy()
img2[~mask] = 255
imshow(img2)

you need to play around with the blur size and the threshold. But in general, this is one way to go.


Break higher order filters down into cascaded second order sections (or biquads). This is covered in every text book on DSP such as https://ccrma.stanford.edu/~jos/filters/


Computing the gains of a piecewise affine system

Add a DC signal to your sinusoidal input, so that the input falls in the middle of each line segment of the input-output curve. Then choose the amplitude of the input sinusoid in order to avoid "spilling over" into the other segments. Measure the amplitude of the output sinusoid (neglecting its DC component, of course). Lastly, compute the gain.


What is the expression for intensity of Laser beam focussed at a point at an angle theta?

Methodology
We can apply simple ray optics to this problem, and it gets easy!
What you're describing is an Orthographic Projection (OP).
So, what you'll do is

realize circles subject to OP become ellipses
take the center of the circles of equal intensity (in your case, $x^2+y^2=\text{const.}$), and project that. It'll be the center of all the ellipses of equal intensity.
figure out what the major and minor axes of the ellipses are going to be
simply write down the formula for points on ellipses with given axes.

Orthogonal Projection
For the math behind OP itself, I'd like to refer you to the wikipedia link given above. It's really just a single matrix multiplication, if your original (circular) image plane and the projected plane both contain the origin $(x,y,z)=(0,0,0)$; it becomes especially easy if the intersection of these two planes is e.g. the $y$-coordinate axis, which we can assume here; you can calculate the normal vector of the slanted plane simply by multiplying the cosine (or sine, depending on your construction) of the incident angle to unit vector on the $z$ axis, then.
Major axis of the ellipses of equal intensity
For the question of what the major axis is going to be: it's the projection of the diameter of the original circle that's parallel to the plane intersection (source). Luckily for you, you've defined that intersection to be a diameter of the circle! Since all points on the intersection of original and image plane are invariant to the OP, you have to do nothing to find the major axis of any circle of given intensity: it's simply the distance on the $y$-axis of the original!
Minor axis
The minor axis is simply perpendicular. You'll find that the projection is also especially easy, and is simply the "complementary" trigonometric function to what you've used to find the normal of the projection plane.
Overall formula
Now that we know how to convert perpendicular incidence to slanted incidence circles, we just have to figure out how to do the inverse. Luckily, since the projection is such a well-formed matrix, that's just taking the inverse– and that's trivial, given the structure of the matrix!


The partial fraction form helps in calculating the z-transform inverse since we can get the inverse of each term in partial fraction by inspection and hence get the inverse of the whole transfer function. The factored form can directly give poles and zeros by equating the numerator and denominator to zero.


Why does the control effort increase if the bandwidth of the system is higher then the bandwidth of the plant?

Intuitively (and I leave it as an exercise for you to work it out with math), it's because you have to push harder on the plant to get it to move fast.

Consider a bowl of honey and a spoon -- moving the spoon around in the bowl slowly takes little effort, but moving it fast requires a lot of effort, because honey is viscous.

It may be most informative to take a simple plant model, like a 1st-order low-pass filter, and without ever designing a closed-loop controller, design a signal for it to follow that has significant energy content above it's rolloff frequency.  Then look at the relative amplitudes of the components of that drive signal compared to the resulting components of the motion of the plant.

Even in open loop your control effort is higher at higher frequencies, for reasons very much analogous to that spoon in honey.


Ok, I found out from a book that linear prediction always results in minimum phase filters, which implies stability.

This can be shown using the mean-squared error (MSE) criterion which is to be minimized in the autocorrelation method. The proof is done by contradiction; if an all-pole filter had a single pole outside of the unit circle, the derivative of the MSE in regards to the pole radius is always greater than zero. In terms of the MSE, this implies the pole position is not optimal, so it cannot be a result from linear prediction (under ideal conditions, I assume).

The actual proof can be found here.


Deconvolution to Remove Gaussian Blur in 1D Signal (Wiener Filtering?)

Approaches
There are many methods for Deconvolution (Namely the degradation operator is linear and Time / Space Invariant) out there.
All of them try to deal with the fact the problem is Ill Poised in many cases.
Better methods are those which add some regularization to the model of the data to be restored.
It can be statistical models (Priors) or any knowledge.
For images, a good model is piece wise smooth or sparsity of the gradients.
But for the sake of the answer a simple parametric approach will be takes - -Minimizing the Least Squares Error between the restored data in the model to the measurements.
Model
The least squares model is simple.
The objective function as a function of the data is given by:
$$ f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$
The optimization problem is given by:
$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$
Where $ x $ is the data to be restored, $ h $ is the Blurring Kernel (Gaussian in this case) and $ y $ is the set of given measurements.
The model assumes the measurements are given only for the valid part of the convolution. Namely if $ x \in \mathbb{R}^{n} $ and $ h \in \mathbb{R}^{k} $ then $ y \in \mathbb{R}^{m} $ where $ m = n - k + 1 $.
This is a linear operation in finite space hence can be written using a Matrix Form:
$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$
Where $ H \in \mathbb{R}^{m \times n} $ is the convolution matrix.
Solution
The Least Squares solution is given by:
$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$
As can be seen it requires a matrix inversion.
The ability to solve this adequately depends on the condition number of the operator $ {H}^{T} H $ which obeys $ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)} $.
Condition Number Analysis
What's behind this condition number?
One could answer it using Linear Algebra.
But a more intuitive, in my opinion, approach would be thinking of it in the Frequency Domain.
Basically the degradation operator attenuates energy of, generally, high frequency.
Now, since in frequency this is basically an element wise multiplication, one would say the easy way to invert it is element wise division by the inverse filter.
Well, it is what's done above.
The problem arises with cases the filter attenuates the energy practically into zero.  Then we have real problems...
This is basically what's the Condition Number tells us, how hard some frequencies were attenuated relative to others.

Above one could see the Condition Number (Using [dB] units) as a function of the Gaussian Filter STD parameter.
As expected, the higher the STD the worse the condition number as higher STD means stronger LPF (Values going down at the end are numerical issues).
Numerical Solution
Ensemble of Gaussian Blur Kernel was created.

The parameters are $ n = 300 $, $ k = 31 $ and $ m = 270 $.
The data is random and no noise were added.
In MATLAB the Linear System was solved using pinv() which uses SVD based Pseudo Inverse and the \ operator.

As one can see, using the SVD the solution is much less sensitive as expected.
Why is there an error?
Looking at a solution (For the highest STD):

As one could see the signal is restored very well except for the start and the end.
This is due to the use of Valid Convolution which tells us little on those samples.
Noise
If we added noise, things would look differently!
The reason results were good before is due to the fact MATLAB could handle the DR of the data and solve the equations even though they had large condition number.
But large condition number means the inverse filter amplify strongly (To reverse the strong attenuation) some frequencies.
When those contain noise it means the noise will be amplified and the restoration will be bad.

As one could see above, now the reconstruction won't work.
Summary
If one knows the Degradation Operator exactly and the SNR is very good, simple deconvolution methods will work.
The main issue of deconvolution is how hard the Degradation Operator attenuates frequencies.
The more it attenuates the more SNR is needed in order to restore (This is basically the idea behind Wiener Filter).
Frequencies which were set to zero can not be restored!
In practice, in order to have stable results on should add some priors.
The code is available at my StackExchange Signal Processing Q2969 GitHub Repository.
Resources

Deconvolution of Synthetic 1D Signals - How To?
Estimate the Filter Coefficients of 1D Filtration (Convolution).
1D Deconvolution with Gaussian Kernel (MATLAB).
How Is the Formula for the Wiener Deconvolution Derived?


What is the difference between general microprocessors, DSP and FFT chips?

In addition to the memory differences that David pointed out, DSP chips have multiple multiply/accumulate (MAC) units.  MAC is a very common signal processing operation that is the basis of filter and FFT/iFFT operations.  Some DSP chips that I used in the past had six MAC units, which allowed them to speed up the FFT operations in two ways: they could do six MAC operations at once, and the MAC operation itself would take one clock cycle while on a general purpose microprocessor each MAC would take at least two clock cycles since it would be a multiply followed by an add.


You are right that in scenarios where the channel is unknown and/or time-varying, we cannot implement a fixed matched filter. The first filter in the signal chain is usually an analog band pass filter eliminating out-of-band  noise. Then after A/D conversion one would use an adaptive equalizer that performs part of the function of the matched filter, but that will usually also be optimized for further noise suppression. This optimality criterion for equalizers that combines minimization of intersymbol-interference and noise power is commonly referred to as minimum mean-square error (MSE) criterion.


Here is (I think) what you are supposed to be doing

Draw a diagram of the poles and the unit circle and zoom in on the area around the resonance frequency
Estimate the relative gain as the inverse of the distance between the pole and the unit circle (assuming that conjugate pole is too far away to change the gain much over a narrow frequency range)
Pencil in the -3dB points by, again, using distance as a representative for gain
Ignore the curvature of the circle (since r is close to 1)
Do some trigonometry to calculate the distance between the two -3dB points
Convert this distance to a frequency difference.


Assuming that they are adjacent and have a "red" energy/pixel value distinct from not projected areas, you may use a simple 2D filter, that fits the size of expected dispersion, and that way get the local maximum that represents the center of the collection of points that represents a single projection.