Combining spectrograms with different windows to getting arbitrary time and frequency resolution

Improving the time/frequency resolution balance remains a standard research issue, trying to adapt the STFT, inventing novel/better time-frequency distributions (Wigner, Cohen), using together or combining the above. For instance in music or speech source separation, a couple of STFTs is often used with both short and long windows, to be able to separate   transient from larger signals.

The ways you can combine such SFTF depend a lot on the exact processing or features you want to extract, as detailed by @hotpaw2. Here are a few pointers to references suggesting methods. The first one suggests a point-wise combination of two spectrograms as $$S_c(f,t) = (|S_1(f,t)||S_2(f,t)|)^{1/2},$$ and the second generalizes it to the geometric mean, the minimax or the reciprocal average.  The third encompasses it in a more generic merging approach.

Combined multi-resolution (wideband/narrowband) spectrogram, ICASSP 1991
On Combining Evolutionary Spectral Estimates, 1994
Adaptive Time-Frequency Resolution for Analysis and Processing of Audio, AES 2006


monochrome (1-bit black-and-white) image conversion

I would suggest for this case that you upsample and apply a slight blur and then a sharpen, then apply the threshold operation.  You won't get any more information from the pixel data, it's simply not there. But you'll get a smoother result out of the thresholding operation, and you won't need to dither. The end result is like a photocopier degradation.

Example:

Also, this is what it looks like when you just use a better diffusion dither algorithm ;)


How to filter very low amplitude staircase signal from a background noise

Linear filtering won't work here. If anything works at all then it might be something like total variation denoising, which is a non-linear technique for removing noise while preserving edges. This article is a very good read, and it also shows a Matlab implementation of the algorithm.


fft result order compared to direct implementation of DTFS using text book

Octave returns your same first result:

octave-3.2.3:1> fft([1 1 0 0])
ans =

   2 + 0i   1 - 1i   0 + 0i   1 + 1i

While 4*ifft([1 1 0 0]) results in what you claim is Matlab's result.

Please check you are not confusing direct and inverse FFT.


Zero padding is a useful interpolation tool using sinc function (or kernel). I will below explain it in 3 parts.
First: what is DFT; Second: Zero-padding; Third: Usage of zero-padding.

(a) $\textbf{What is DFT}$: By taking a DFT of a data set we are mapping the data values from the current discrete domain (many case happens to be time example in audio signals; it can also be spatial co-ordinates as for images) to a discrete frequency domain. With the transform doing this mapping being a matrix (hence linear) which is always full-rank the mapping is always invertible without any loss in information content. Example we can take idft of the obatined dft and get back the same data we started with again or no loss!

Why go to this other domain (so called transform (here frequency) domain) of representing the same data?
Well, (a real life note) when we are trying to make decisions we want to see the facts clear in front of us. We don't want unnecessary facts and likes and dislikes to wander in our mind. It similar in signal processing. When we want to analyze data we want to see the data in such a domain of representation where we have,

concentration of information only in localized points in that domain
the domain gives some sort of a physical interpretation to the data we have (this is a super plus point)

Thus it is clear by using DFT we want see if we can gain from the above two. Of course, if the data has one of the columns of the DFT matrix itself we have gained the maximum out of the DFT domain based representation i.e only one point in the mapped domain values will have a non-zero value ideally. Now suppose the signal is not representable using any one of the columns of the DFT matrix. Then DFT will try to represent the data by a minimal number of mixture of the columns of the DFT matrix i.e only few of the mapped domain values will have a non-zero values. So, which type of data will do this: those which share some characteristics with the DFT matrix columns such as are periodic.
$A\ moral$: hence DFT is good to represent periodic data as it gives a minimal representation based on number of non-zeros values (magnitude captures the contribution of each column in making up the original domain data).

(b)$\textbf{ Zero-padding}$: In zero-padding we have first padded zeros to the data in the original domain and then take the new zero-padded signal's DFT.
This means by zero-padding we have increased the number of columns in the DFT matrix (with the matrix now also being orthogonal) with no new data in the original domain being added. With this we can hope to see more in deep how much the new added columns contribute in representing the time-domain signal. Note what we earlier say with N-point DFT will also see the same at those N points in a subsequent 2N-point DFT of the same data. So the peak magnitude will not change. (Considering we are not normalizing the DFT matrix). Now when we take the iDFT normalization factor comes into picture and we have to normalize by 1/2N this time to get the same data back instead of by N in the previous case and not pick only N points in the 2N point data. The zero-padding in time-domain can also be interpreted as sinc interpolation in the DFT domain to get the contribution of the newly introduced columns in the 2N-pont DFT. This needs a separate post of detail description. But you can see it easily by drawing the signal and seeing the zero-padded signal as a multiplication of the original assumed to be 2N length signal with a rectangular window which results in convolution in the frequency domain with a sinc function.
It may be clear now that we are not increasing frequency resolution by zero-padding i.e we are not resolving two frequencies close by but we are filling new frequencies between the existing N-point DFT frequencies and finding their contributions for making the signal. Frequency resolution can be increased only by sampling the data more finely or taking new data points which is in our hand only if we are an experimenter.

(c) $\textbf{Applications of zero-padding}$: It is used as a frequency domain interpolation tool for getting the side lobe structure for filters. Also it is used to interpolate (or re-sampling) in time domain by zero-padding in frequency domain. Zero-padding in frequency domain needs care so as to preserve the original phase of the signal.

Hope this post is clear. May have a lot of redundancy if you already know much of it.


Compare two signals, find times when they are affected by the same processes

You can model your system as a hypothesis testing problem. For example, let $M_1$ and $M_2$ denote the two systems which are correlated. Let the observed signal in $M_1$ and $M_2$ be given as
\begin{align}
y_1(t) = s_1(t) + n_1(t) \\
y_2(t) = s_2(t) + n_2(t),
\end{align}
respectively, where $n_i(t) \sim \mathcal{N}(0, \sigma_i^2)$, for $i = \{1,2\}$, and $s_i(t)$, for $i = \{1,2\}$, is the original signal. Since the two systems are correlated, you can model them as
\begin{align}
\begin{bmatrix} y_1(t) \\ y_2(t) \end{bmatrix} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \begin{bmatrix} s_1(t) \\ s_2(t) \end{bmatrix} + \begin{bmatrix} n_1(t) \\ n_2(t) \end{bmatrix},
\end{align}
where $-1 < \rho < 1$ denotes the correlation constant between the two systems. Now, we can model the system as a hypothesis testing problem as follows:
\begin{align}
\mathcal{H}_0: \rho = 0, &\quad \mbox{uncorrelated model} \\
\mathcal{H}_1: \rho \ne 0, &\quad \mbox{correlated model}.
\end{align}
If you assume that $s_1(t)$ and $s_2(t)$ are known, then find the maximum likelihood estimate (MLE) of $\rho$ under $\mathcal{H}_1$. If you assume $s_1(t)$ and $s_2(t)$ are unknown, then use generalized likelihood ratio test, i.e., find the MLE of $s_i(t)$ under each hypothesis and replace $s_i(t)$ with its MLE.


compute a difference( dB) between two BER curves

For the segment you show, there's no single number that describes what you want – your lines need to be parallel for that to become a single number.
So, your question can't be answered in general from the diagram, as you did, you'll have to take a fixed BER; and then, I don't really understand the question – that's a horizontal line in your plot. You can read the difference directly from where the curves hits that line.
You can, however, talk about an asymptotic gain if you are able to figure out how these curves continue for large SNRs (you'll have to run with sufficiently many points at say around 23 dB Eb/N0, in my experience, for 16-PSK); at some point, these curves become sufficiently parrallel. Then you can measure the difference in Eb/N0 necessary to achieve the same BER. (that's a horizontal line through the graph; you measure the distance of intersection of curves with the horizontal line.)
Note that this does not work the same way for BER achieved with a given Eb/N0. (i.e. you can't just measure the vertical distance, that doesn't converge as well.)


Why the Hilbert transform is not correct using the hilbert function in MATLAB?

I don't think the Hilbert function in Matlab is incorrect - you can establish this by carrying out the same process for a stationary sinusoid (observe the discontinuities at the start/end). Further, the unwrapped phase of the Hilbert transform you've performed does have a sinusoidal characteristic, although with a significant offset.

I think it's more likely you're looking at an interaction of your rapidly changing phase with the windowing implicit to a filter implementation of the Hilbert Transform. Perhaps you could try slowing down the rate in change of phase to see whether you can isolate changes in phase from the window length? You could do this by expanding your time vector (using an increment smaller than 0.001).


Where should I set my anti-aliasing filter corner frequency for this signal?

Suppose if your sampling rate is fs, then the corner frequency of the Anti Aliasing filter would be fs/2


How does the $\mathcal Z$-transform's "region of convergence" work?

To be completely honest, I thought the theory behind the Z-transform was kinda opaque in college too. In hindsight, taking a course in complex analysis would have made it clearer. And I too dislike the notational conventions that seem to be used for this stuff. Strictly speaking, the usual convention here is that

$x[n]$ denotes a discrete-time sequence

$n\in \mathbb{Z}$
the brackets denote a discrete argument

$X(z)$ denotes a continuous-valued transformed function

$z\in\mathbb{C}$ (is a complex number)
the parentheses denote a function accepting a continuous-valued parameter
the capital $X$ denotes a transformed version of some other function/sequence $x$ (a similar notation is used for Fourier transforms: $F(j\omega)\leftrightarrow f(t)$

  What do they mean by z=0? Are they substituting z as 0, if so in which equation?

They mean, just plug $z=0$ into your usual definition of the Z-transform.

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

Generally (more precisely, when $x[n] \ne 0$ for some $n \ne 0$), this sum will diverge (to infinity) for some complex $z$. For instance, let $x[0] = 1, x[1] = 1$, and $x[n]=0$ for $n<0$ and $n>1$. Then $X(z) = 1 + z^{-1}$. The ROC does not includes $z=0$, for $\lim_{z\rightarrow 0} X(z)=\infty$

When your text says "When $n_2>0$ there will be a $z^{−1}$ term and thus the ROC will not include $z=0$", what they mean by that is, when $x[n]$ is nonzero for some $n>0$, it is unavoidable for the z-transform to include a $z^{-n}$ term, which diverges to infinity at $z=0$. That's all.

  How do we calculate the region of convergence for an infinite sequence?

Lots of math. Ha!

srsly, the way this is done is to obtain an algebraic formulation for the sequence in question, plug it into the Z-transform definition, and use the tools available from the analysis of geometric series (and complex power series) to determine where this Z-transform converges/diverges. In practice, determining whether $|z|=1$ converges is the most important question to answer, because that determines stability, and whether or not you can obtain a frequency response from the system, etc. But causality might matter too, depending on what you're doing.


You're running a DAC at 5 MS/s and have a reconstruction filter of 50 MHz bandwidth.

There's your error: You're not doing proper reconstruction filtering! Reduce the filter bandwidth to reflect your sample rate.

You can of course also interpolate your signal 10× and set your DAC clock to 50 MS/s. Interpolating an OFDM signal is potentially very easy: Just add enough unused carriers to the side (note: you might need to adjust your CP length and synchronization algorithms).


Controlling audio high/mid/low in frequency domain

while you could do this in the frequency domain, there are many problems with this approach, some of which are outlined here.

A better approach is to use a time domain filter. There is a tutorial here.


Nyquist Frquency isn't what it is supposed to be

The closer the sampling frequency gets to half the (highest spectral) frequency of the signal, the longer one has to sample to get a sufficient bounded fit (for some bound) of a reconstruction using those finite number of samples.  This required sampling time extends to infinity as one approaches Fs/2.

So, you simply didn't sample for long enough in your experiment.  Try 5 thousand and 5 million samples for sampling a 12.9 Hz sinusoid at 26 sps, and you will see a closer fitting looking plot over most of an X axis of 5 thousand or better million reconstructed sample points.  Then try a sinusoid even closer to Fs/2 and you will find you need to sample even longer to get the same visual quality of fit.


Obtaining the size of objects from an image

I would use a split-and-merge algorithm for segmenting and then simply count the pixel per segment. Details are provided here

Two assumptions are neccessary:

1) you objects do not overlap.

2) the picture does not have any background or further irrelevant objects


How to Subtract One Signal from Another?

Under the assumption that your background noise is statistically independent of the signal you want to recover and that your second microphone picks up virtually nothing from the body noise, you have the following model for your measurements:
$$ \begin{eqnarray} M_1(t) &=& L\{N(t)\} + B(t) \\ M_2(t) &=& N(t) \end{eqnarray}$$

Here, $M_1,M_2$ are the signals picked up by the microphone, $B(t)$ is the signal you want to recover, $N(t)$ is the environmental noise signal and $L$ is a linear time invariant system that maps the recorded environmental noise to what is picked up by the second microphone.

With $\langle\langle \cdot, \cdot \rangle\rangle$ as cross correlation, we also have the non-correlation constraint $\langle\langle N,B \rangle\rangle=0$ and $\langle\langle L\{N\},B \rangle\rangle=0$ as a consequence of statistical independence of $N$ and $B$.

These prerequisites guarantee a cross correlation matching pursuit to converge against $B(t)$:

Normalise $M_2$ with the standard inner product: $M_2(t)\leftarrow \frac{M_2}{\sqrt{\langle M_2(t),M_2(t) \rangle}}$
Find $c(t):=\langle\langle M_2,M_1\rangle\rangle$
Find the time associated with the maximum of the squared cross correlation: $\tau=\mathrm{argmax}\left( c^2(t) \right)$
Update $M_1$ to remove the contribution of the environmental noise: $M_1(t)\leftarrow M_1(t)-\langle M_2(t-\tau),M_2 \rangle\cdot M_2(t-\tau)$
If result not good enough, go back to 2.

In step 4 you may also use $c(\tau)$ instead of the inner product if the cross correlation uses the same inner product. Stopping the matching pursuit is a question of experience. Some of the assumptions you need for this algorithm are only approximately true, like the perfect uncorrelatedness of the environmental noise and the body signal. A good way to check if your matching pursuit runs into a non-productive state of signal modifications is to observe the natural square-norm of $M_1$. If it stops decreasing in step 4, terminate the loop.

After the algorithm ran, $M_1(t)$ will contain only your body signal $B(t)$.


Can you downsample unevenly with even distribution?

Absolutely positively -- uh -- maybe.
It depends entirely on the underlying process that generated the original series, what the series "means", how much unique information is actually present in those ten samples, how much you're willing to throw away, and what you know about how that information is structured.
At the extreme "no" end would be if those ten data points were essentially completely random -- i.e., if you recruited ten people on the street and asked them their favorite number, then there's no way of relating the 1st sample to the second.
The trivial extreme "yes" would be if those ten data points were super highly correlated -- the altitudes of ten points a meter apart on a nice smooth sidewalk, for instance.  Then your proposed simple interpolation scheme would work quite nicely.
What you're trying to do is a generalization of downsampling.  If you had an infinite sequence of points that it made sense to describe in the frequency domain, then as long as the sequence was band-limited to less than 7/20th of the original sampling rate, you could upsample by a factor of 7, filter the resultant signal, then downsample it by a factor of 10.
You can't do that, because filters suffer from end effects, and from the point of view of the necessary filters, your data sequence is all end and no middle.
The absolute best formal way to do the job would be to understand whatever the underlying process that generated those ten points.  From that, figure out how to work backward to the most likely location or the seven points, and voila, you either have an answer or your loved ones retrieve your body from under a mountain of scribbled on paper and give you a decent burial.
A pretty good ad-hoc way to do this, if you have reason to believe that those ten points represent some sort of a bandlimited process, is to just use a set of splines to make a smoothed line between the ten points, and then sample that line at the seven points.  This really only works well if the underlying process happens to match up well with a spline -- but if it's bandlimited it does -- and if it's something like the positions of a beam at those ten points then it definitely does, because the mathematical sort of spline is a pretty good simulation of the draftsman's spline, which is a beam that's pinned down at certain points to make a sparse set of points into a nice smooth curve, for boatbuilding.


So far I have not experimented with audio steganography, but anyway, maybe these ideas are of some use for you.

Adding echoes to a signal is pretty easy as this process just adds a delayed version of the signal back to itself. You probably want to use a gain factor $f_\text{gain}$ for the additional signal that is smaller than one, $f_\text{gain} \leq 1$, to avoid having an echo that is louder than the original sound. In DSP terms this process corresponds to FIR filtering with an impulse response that consists of a delta impulse at time zero ($\delta[n]$) and a time-shifted and weighted delta impulse:
$$
\mathbf{b} = \delta[n] + f_\text{gain}\cdot \delta[n-L_\text{delay}],
$$
with $L_\text{delay}$ the delay time in samples
and
$$
\mathbf{a} = 1
$$
because it's an FIR filter. To get multiple echoes simply use more than one delayed delta impulses.

For detecting hidden echoes in the signal, maybe it's worth a try whether you can achieve satisfactory performance with an adaptive filtering algorithm. For example, the (normalized) least mean-squares ([N]LMS) algorithm comes to my mind. I'm sure you can find implementations in a number of programming languages to download from the internet. Using the signals both as input signal and desired signal to the adaptive filtering algorithm would lead to a self-cancelling behaviour. Maybe, from analyzing the filter coefficients, you can find out about whether there are echoes in the signals and, if so, what the delay times are.

If you experiment with the Cepstrum method described in the web page that the link in your question points, take care to have your signal blocks long enough to capture the echoes.


In the context that you described, $t$ is a variable indicating time. That is, the wave will take on a value of $y(t)$ at the time instant $t$. If this is measured in seconds, then $\omega$ specifies the number of radians that the wave passes through in one second. And, as you noted, if the wave has a frequency of $\frac{100}{15}\ \text{Hz}$, then it will pass through 100 periods in 15 seconds.


Wavelet analysis results in information about both frequency and time localization.  Since the two filters are slightly offset in time, if your signal varies over time, one might expect the two wavelet results to be different.


I am not sure I understand your question, but for the particular question

  How are features of the EEG Signal extracted then?

CSP gives you a matrix $W$ that represents a spatial transformation of your EEG $X$. When using CSP, $W^TX$ will have the property of maximal variance for one class while minimal variance on the other class. Therefore, when using CSP for classification, the features are the variance of $W^TX$. In practice, the log-variance has a better distribution.

An interesting and brief explanation of CSP can be found on doi:10.1109/MMSP.2010.5662067 (until section II, then it goes on Riemannian geometry which might be out of your scope).


Power Spectrum Density (PSD) represents the spectral distribution of total power of a given WSS random signal. Conventionally, it's defined as the Fourier transform of the Auto-Correlation function (sequence) of the WSS random process as:

$$ P_x(\omega) = \mathcal{ F } \{ r_{xx}[m] \} $$

where the auto-correlation sequence (ACS) $r_{xx}[m]$ of the WSS $x[n]$ is defined by:

$$ r_{xx}[m] = E\{ x[n]x^*[n-m] \} $$

These are theoretical entities, you may not compute them in practice where only a single or a few realizations of the random process are available. In such a case, the PSD is estimated and there are different estimation techniques, one of the most common and simplest is known as the periodogram estimate of PSD and given by

$$ \hat{P}_{x}(\omega) = \frac{1}{N} |X(\omega)|^2 $$

where $X(\omega)$ is the discrete-time Fourier transform of the observed signal $x[n]$ of length $N$.


DFT of sum of sinusoids with random zeroed samples

The math is well known; it is the convolution theorem for the DFT.

In this specific case:

$$DFT\left\{f[n]\cdot z[n]\right\} = DFT\left\{f[n]\right\} * DFT\left\{z[n]\right\}$$

Where:

$f[n]$ denotes your original sequence of 128 samples

$z[n]$ denotes a sequence of 128 samples, each randomly selected from $\{0,1\}$

'$*$' denotes cyclic convolution, since we're talking about the DFT.

The results of this convolution are visible in your second plot, particularly at very positive and very negative frequency bins.


Fast/Slow Fading: How are doppler shift and coherence time related?

Here's an intuitive explanation.
Imagine a communication system (wireless) in a static environment (maybe picture the inside of a closed warehouse). There will be an amount of fading, but it is static: it never changes, because the channel doesn't change.
For the channel to change, the multipath scenario must change, and that means that objects must move around (it can be the reflectors/scatterers, the antennas, or both). If they move relatively slowly, you have "slow" fading; if they move relatively quickly, the fading becomes "fast".
Now, whenever you mix waves and moving objects, you get Doppler shifts. Furthermore, the faster the movement, the higher the Doppler. In consequence, you can measure how fast the channel is changing by either measuring how fast the objects are moving, or by how large the Doppler is.
Since the Doppler shift is easy to quantify, it is the preferred way to measure how fast a fading channel is changing.


A basic question in FIR filter design using MATLAB

All digital filter frequency parameters (passband begin frequency, passband end frequency, and stopband begin frequency) are stated in terms of an input signal sequence's Fs sampling frequency. For a lowpass filter example, if I said a lowpass filter's passband width (it's "cutoff" frequency) is 0.2, I'm saying that the cutoff frequency is 0.2 times Fs. So if an input signal's Fs was 1.0 kHz that lowpass filter's cutoff frequency is 200 Hz. However, if an input signal's Fs was 20 kHz then in this situation that same lowpass filter's cutoff frequency is 4 kHz.

For a different scenario, if you said, “Rick, please use Matlab and design a digital lowpass filter whose cutoff frequency is 300 Hz.” It would be impossible for me to do that unless you told me what was the filter’s input sequence’s Fs sample rate. But if you then said “Fs = 1000 Hz”, then I’d know the filter’s ‘normalized’ cutoff freq was 300/1000 = 0.3 enabling me to design the lowpass digital filter.


sampling a sinusoidal signal at the Nyquist rate

Noha, when the classical sampling theorem (Nyquist-Shannon) for baseband signals
is applied to an ideal sinusoidal of infinite extend of frequency $f_m$ Hertz
$$x(t) = \cos(2 \pi f_m t)$$ it's found that the signal must be sampled at a sampling frequency greater than $2 f_m$; i.e.
$$f_s > 2 f_m$$

Note that even though the selection of $f_s = 2 f_m$ would sufficiently work for most of the typical engineering signals, it will not work for an ideal sinusoidal signal. The reason is that the ideal sinusoidal has an exact impulse at the $f_s / 2$ frequency whereas most typical signals will not.

The practical solution is to use a sampling frequency slightly greater than the theoretical value given by the sampling theorem unless you are sure that there wont be any components at $f_s/2$.


Find the impulse response from the frequency response

HINT: (because this is a homework type problem)

Using the sign function the frequency response can be expressed as

$$H(\Omega)=-j\,\textrm{sgn}(\Omega)\tag{1}$$

Direct computation of the inverse Fourier transform is tricky because of convergence issues of the corresponding integral.

If you can use Fourier transform tables, find the inverse Fourier transform of the step function $u(\Omega)$ in the frequency domain. If you can't find it, use the Fourier transform of the step function in the time domain $u(t)$, and figure out a simple transformation to obtain from it the inverse Fourier transform of $u(\Omega)$. Then express $\textrm{sgn}(\Omega)$ in terms of the step function $u(\Omega)$.

As a final hint, that system is called a Hilbert transformer.


Confusion between samples per second and samples per symbol

A digital communications signal such as QAM can be generally written as $$s(t) = \sum_{n=0}^{N}a_n p(t-nT),$$ where $a_n$ are (possibly complex) amplitudes taken from a finite alphabet, and the symbol rate is $1/T$. The pulses $p(t-nT)$ are often called "symbols". Matlab commands that ask you to specify "samples per symbol" are asking how many samples you want to use to define $p(t)$.

Of course, since you know the pulse duration, specifying the samples per symbol is an indirect way of specifying the sampling rate.


Why FFT of a filtered step increases the lower values of the transform with respect to the unfiltered signal?

So taking your filter and looking at its magnitude response we see that the magnitude close to zero will increase but all other magnitudes (on the other side of the dip) should be approximately unity gain (no magnification / increase in amplitude).

The fact that the DC gain of the filter is 0.94 means that the steady-state output for a step input will be 0.94, not 1.  That may be the effect that you're seeing.

Also, the notch is not very deep. You may want to redesign the filter to make it attenuate more.


Stability of the Low pass filter after Spectral Transformation

The filter is stable if the all the poles of the transformed filter, i.e., filter obtained after applying spectral transformation, lie inside the unit circle. For example, if $H(z)$ is the original filter and $F(z)$ is the filter response of the transformation filter, then the spectrally transformed filter response is given as
\begin{align}
G(z) = H(F(z)) = H(z)|_{z^{-1}\leftarrow 1/F(z)}
\end{align}
where
\begin{align}
\frac{1}{F(z)} = \frac{1-\beta z}{z - \beta}.
\end{align}

Check if the poles of $G(z)$ lie inside the unit circle. If yes, then $G(z)$ is stable. Note that $F(z)$ is an all-pass filter. Therefore, the regions of stability and instability of $H(z)$ are preserved in $G(z)$.


Image Processing method for spiky selection

There are more corners on the borders of your "spiky object", so one approach would be to tune a corner detector for this.

For example, I calculated the determinant of the structure tensor (Mathematica code below) of a distance-transformed image:

Binarizing with hysteresis yields this image, which should be a good starting point for the segmentation algorithm of your choice:

Mathematica code (src is the source image you posted)

At first, i calculate a distance transform of the input image. This creates contrasts over the whole object area (instead of just the border), so the whole object can be detected.

dist = ImageData[DistanceTransform[src]];

Next I prepare the components of the structure tensor. The filter size for the gaussian derivatives if 5, the window size is 20.

gx = GaussianFilter[dist, 5, {1, 0}];
gy = GaussianFilter[dist, 5, {0, 1}];
gx2 = GaussianFilter[gx^2, 20];
gxy = GaussianFilter[gx*gy, 20];
gy2 = GaussianFilter[gy^2, 20];

To calculate the corner filter at each pixel, I simply plug these into the symbolic determinant of the structure tensor:

corners = Det[{{dx2, dxy}, {dxy, dy2}}] /. {dx2 -> gx2, dxy -> gxy, dy2 -> gy2};

Which is basically the same as:

corners = gx2 * gy2 - gxy * gxy;

Converting this to an image and scaling it to 0..1 range yields the corner detector image above.

Finally, binarizing it with the right thresholds gives the final, binary image:

MorphologicalBinarize[Image[corners], {0.025, 0.1}]


How to pad a array representing music to pad the stft of this array?

After some reading, I found the solution to my problem. In the default mode, windows overlap on half their size. Scipy starts by adding zeroes at the boundaries on the left and on the right so that the first and last window are centered on the first and last point of data. Then it pads the end with zeroes so that the last windows on the right is full of point. So, in order to add a window on the left/right, you have to add nperseg/2(127 in my case) points of data. So to add 127 windows on the right and on the left, you just need to add 127^2 = 16129 points on the right and of the left of the original data.

Edit : Based on how the stft is computed, if the original signal is zero-padded on the right and left, any new window created will be on data points with value 0, and as a result will have a value of 0. This means that it is exactly equivalent to 0-pad the original signal of 127*n points or to 0-pad the stft of n points.


pitch detection: which is the most robust way to distinguish voiced from unvoiced audio

SWIPE is a pitch estimator that worked well for me in the past. You can find dissertation & code in MATLAB here. It might give you some ideas.


Is Image Noise Removal Image Enhancement or Image Restoration Operation?

It depends on the final goal and the model used.
If we assume we have an image - $ I $ which is a result of a sum on of an image $ W $ with noise $ n $: $I = W + n$ and our goal by reducing the noise is estimating $ W $ than this process is restoration.
Yet I could try achieving a different goal by denoising, such as making it look better (Subjective target) or use prior such as Piecewise Constant image which means the denoising process if more of an enhancement to the image.


What are some python packages I can use to cut audio files

The librosa library is usually used. You can find the files here:

https://github.com/librosa/librosa


The Fourier Series Of This Triangle Wave

Check out Wikipedia's Fourier Series page.  Their "Example 1" shows how to derive the Fourier series of a sawtooth wave.  All you have to do is normalize the results for your particular time and amplitude values.


Selecting a Digital Signal Processor for a Specific Image Sensor

Good approach to first do a rough calculation of the bandwidths you need.
Couple of remarks on that:

You forgot a factor of 3, that camera has three "color pixels" per image pixel
Your use case screams "I should be using a commercial off-the-shelf USB camera"; don't engineer something very complex if it doesn't have a value proposition over what you can buy. If the value proposition is "I can learn how to work with cameras, and DSP", I'd frankly start smaller, and learn how to work with camera images on a PC in software before moving out part after part to an embedded processing platform.
You try to downplay the camera interfacing, but it's really a critical part here: in need to interface digital buses for which their processors have no dedicated interface, engineers turn to FPGAs. You'll need one, just as the "glue" between your camera
Your approach of using an MCU betrays that you're not overly familiar with what these can and can't do. Wrong order of magnitude: for this kind of problem you need what is usually sold as "application processor", with plenty RAM, and things like interfaces for network hardware and the opportunity to run an OS reasonably fast.
You picked an oldish low-power DSP core. That 400 MHz thing can do way less than modern application processor CPUs can do. For example, every SSE-enabled x86 CPU can do two multiplications at once, just like that DSP core, but they tend to run much faster.
You're planning to stream something between 1080p and 4k video .. to the internet. You really need a fast internet uplink to do that. You might want to download a 4k video snippet and look at the amounts of bits per second that needs.

So, all in all:

You need an application processor, which take your MCU completely out of business (you will probably still need an MCU, to do things like turning on power to the different power rails of your other components in the right order, but not at all involved in handling the camera data)
You need a device with a MIPI interface. If you're able to buy in big numbers, you'd typically get a camera or smart phone System-on-Chip that has a CPU (see above) which has such an interface. In your case, it's going to be an FPGA, or a MIPI-to-USB3 bridge (cypress sells these). In any way you do it, getting the interface right isn't trivial.
I don't really see how a DSP core fits in here in practice: for the data / camera handling, your application CPU does the most. For the video encoding, you'd typically use a dedicated encoding accelerator that does the computationally hardest part of it.

All in all, this sounds like you'd start with a small x86 mainboard (there's way smaller form factors than µATX, but honestly, µATX isn't so large for prototyping), and get a camera that you can attach to that - and that would usually simply be a high-quality webcam or DSLR.
For highly integrated systems, going the FPGA route is very attractive, and you can do a lot of computation by writing the right hardware description - problem is, that takes a lot of experience. For even higher integration, you'll find ASICs that are both pretty camera-specific and use-case-specific that take on the task of handling everything from camera to radio interface. In fact, DJI has bought a semiconductor manufacturer just so that they have that kind knowledge in-house, although they are by no means a semiconductor company – they sell video drones.
Honestly, I'd call the FPGA approach "unrealistic within the next year", given your current understanding of the matter, sorry. The ASIC route is impossible, since production of a single high-speed wafer costs millions, and you don't have a silicon design team.


A propagation loss expressed in decibels is equivalent to an attenuation in linear scale. If you have a loss of l_dB (in dB), that's equivalent to an attenuation of l = 10^(l_dB/20). So, if your transmitted signal is S, the received signal is R = S ./ l.

Note that the factor of 1/20 in the conversion from dB to linear scale is needed because we're scaling the amplitude (voltage) of the signal.


In the context of signal processing, a continuous-time system

$$y(t) = \mathcal{T} \{ x(t) \}$$

is said to have memory if its output at time $t=t_0$ has a dependency on the input values for any time other than $t=t_0$. Otherwise; if the output at time $t = t_0$ depends only on the input value at time $t=t_0$, then it's a memoryless system.

So the system defined by:
$$y(t) = x(2t)$$
has memory, which can be proven by a single example, as $y(t)$ at $t=5$ is given by $x(t)$ at $t=10$...

Furthermore, a system is invertible if unique (distinct) inputs produce unique (distinct) outputs, otherwise it not invertible.

So the continuous time system defined by:
$$y(t) = x(2t)$$
is invertible, as its inverse is given by $x(t) = y(t/2)$.


These problems are most easily solved by using the $\mathcal{Z}$-transform. But it is also possible to do it in the sample domain. You have two equations for expressing $w[n]$ and $w[n-1]$ in terms of $x[n]$ and $y[n]$:

$$x[n]=w[n]-aw[n-1]\tag{1}$$
$$y[n]=w[n]+bw[n-1]\tag{2}$$

From $(1)$ and $(2)$, $w[n]$ is obtained as

$$w[n]=\frac{1}{a+b}\big(bx[n]+ay[n]\big)\tag{3}$$

Plugging $(3)$ into $(2)$ gives

$$y[n]=\frac{1}{a+b}\big(bx[n]+ay[n]+b^2x[n-1]+aby[n-1]\big)\tag{4}$$

Bringing all terms with $y[n]$ to the left side gives

$$y[n]\frac{b}{a+b}=\frac{1}{a+b}\big(bx[n]+b^2x[n-1]+aby[n-1]\big)\tag{5}$$

which finally results in

$$y[n]=x[n]+bx[n-1]+ay[n-1]\tag{6}$$


You can describe a frequency response in terms of its real-valued amplitude and its phase:

$$H(f)=A(f)e^{j\phi(f)}\tag{1}$$

Note that $A(f)$ is not the magnitude, but a bipolar amplitude function. Equivalently, you can express $H(f)$ in terms of its magnitude and its phase:

$$H(f)=M(f)e^{j\tilde{\phi}(f)}\tag{2}$$

Now we have $M(f)=|A(f)|\ge 0$, and, consequently, the phase $\tilde{\phi}(f)$ in $(2)$ is generally different from the phase $\phi(f)$ in $(1)$. For frequencies for which $A(f)<0$ is satisfied, we have to add $\pi$ (or $-\pi$) to the phase to compensate for the sign change:

$$\tilde{\phi}(f)=\begin{cases}\phi(f),& A(f)\ge 0\\\phi(f)\pm\pi,&A(f)<0\end{cases}\tag{3}$$

For a linear phase filter of type I or type II (even symmetry), the phase $\phi(f)$ is a linear function of $f$, where $\phi(f)$ is defined as in $(1)$, i.e., with a bipolar amplitude function. The phase $\tilde{\phi}(f)$ as defined in $(2)$ has jumps wherever $A(f)$ changes its sign.

You can extend the definition to type III and type IV linear phase filters (odd symmetry) by expressing $H(f)$ as

$$H(f)=jA(f)e^{j\phi(f)}\tag{4}$$

where $A(f)$ is again a real-valued bipolar amplitude function, and $\phi(f)$ is a linear function of $f$.

Note that by linear function I mean a strictly linear function without an additive term, i.e., $\phi(f)=a
\cdot f$, with $a=-2\pi\tau_g$, where $\tau_g$ is the filter's group delay.


Denote the array manifold (what you call array pattern) by a $N \times 1$ vector ${\bf{a}}\left ( \theta \right )$ where $\theta$ is the DOA and $N$ is the number of array elements. Further denote by the $N \times 1$ vector $\bf{x}_k$ the $k$-th snapshot of the array, then you are measuring the following signal:

\begin{equation}
{\bf{x}_k}= {\bf{a}}\left ( \theta \right ) s_k + {\bf{n}}_k
\end{equation}

where $s_k$ is the signal reflecting from the target and $n$ is a random noise.

Denote the beamformer by a $N \times 1 $ vector of weights ${\bf{w}}$, then the array response after beamforming would be

\begin{equation}
r(\theta) = {\bf{w}}^H {\bf{x_k}},
\end{equation}

Note that $r$ is a scalar for a given DOA $\theta$ and $\bf{w}$ is a function of $\theta$. Now, to get the array response for a variety of DOAs just change $\theta$.


Is the DCT prone to spectral leakage like the DFT?

A window function other than rectangular can be applied to suppress sidelobes also with the discrete cosine transform (DCT). Window functions are also sometimes used together with some flavors of DCT as a lapped transform, with pre and post windowing providing some protection against reflection artifacts arising from processing (such as quantization) of the DCT data.
Here is tested to multiply a real sinusoidal input by 1) a rectangular window or 2) a Hamming window, before DCT analysis. The tested inputs are cosines with different phase shifts. Phase shifts that are multiples of 90 degrees give identical results in the magnitude of DCT. Each plot shows the magnitude of DFT in decibel scale for multiple different fractional offsets of the sinusoid frequency from a DCT bin frequency, in units of bin width. Octave's dct was used for the calculations.
DCT of cosine with 0 degree phase shift:
Rectangular window:

Hamming window:

DCT of cosine with +45 degree phase shift:
Rectangular window:

Hamming window:

DCT of cosine with -45 degree phase shift:
Rectangular window:

Hamming window:

Note that depending on the phase of the real sinusoid, its DCT, and DCT of the sinusoid multiplied by a window function, will be more or less spiky, which is another phenomenon than the window's much lower sidelobes.
Octave source:
pkg load signal;
phase = 2*pi/8; # Phase shift of the cosine
N = 256; # Sequence length
M = 8; # Number of fractional shifts that are in range 0..1 of DCT bin width
n = [0:N-1]'; # Zero-based index
omega = 2*pi*(0.25 + 0.5*[0:M]/M + 0.5*round(N/M)*[0:M]); # Cosine frequency
x = 2*real(exp(i*(n.*omega/N + phase))); # Cosine with phase shift
plot(n/(M*round(N/M)), 20*log10(abs(dct(x)/sqrt(N))));
xlim([0, 1]);
ylim([-60, 0]);
title(["rectangular window, cosine phase shift = " num2str(360*phase/(2*pi)) " deg"]);
xlabel("bin (not labeled), for each peak: fractional offset of cosine from bin center");
ylabel("|DCT(x)[k]| (dB)");

y = x .* hamming(N);
plot(n/(M*round(N/M)), 20*log10(abs(dct(y)/sqrt(N))));
xlim([0, 1]);
ylim([-60, 0]);
title(["Hamming window, cosine phase shift = " num2str(360*phase/(2*pi)) " deg"]);
xlabel("k (not labeled); for each peak: fractional offset of cosine from bin center");
ylabel("|DCT(x)[k]| (dB)");


Finding audio "tics" in stereo wav files using python

I don’t know what DAW you use, but I’m pretty sure Waves and probably some other companies have plug ins for this kind of application.
After reviewing your images it seems that the tics are high frequency and high amplitude so there are some simple things you could try.
You could check each sample and compare the absolute value to some threshold, such that if it is higher, you have a tic.
Second, you could high pass filter the signal to isolate the tic, then do the previous step.
Third, you could calculate the RMS of the signal, then compare the absolute value (peak) of each sample to the RMS, such that a tic occurs whenever the ratio or difference exceeds a given threshold.
Lastly, you could do frequency magnitude analysis.  For this you would set up a window size N and step size M.  Every M samples you take N samples and apply a window function, then calculate the frequency transform magnitude of the N samples.  One idea would be to look at the rate of change of individual frequency components such that if it exceeds a threshold you have a tic.  This is probably the most CPU intensive and complex implementation.


Is it possible that adding noise to a signal can increase detection possibility?

Our ears have a hearing threshold- i.e. a minimum amount of sound power that is required for them to detect.  Thus, a pure tone that is less than the threshold will not be detected.  If white noise is added to the tone such that the total power of the tone and noise exceeds the threshold, then the ear can hear the tone.

Perhaps you have some sort of similar non-linear effect going on.


How to smooth the echo duration change to get rid of audio clicks?

Maybe, when $D$ changes, you can crossfade from the audio at the old $D$ to audio at the new $D$.
If you ran a pitch detector alongside of the delay (and stored the detected period length in a corresponding array) you could crossfade one period at a time until you got to within a period of your new $D$.  Then slide very slowly to the exact $D$.


Using MATLAB's Spectrogram Function for Analysis

See the MATLAB documentation: s = spectrogram(x) returns the short-time Fourier transform of the input signal, x. Each column of s contains an estimate of the short-term, time-localized frequency content of x.. Namely each column of the matrix s is the result of an fft() on some samples of the input. So the plot you see is the magnitude of the columns of s.

Spectrogram is about analysis of Non Stationary signals. So something is changing over time which means it makes no sense to look on the DFT of all samples. The window length is the time you think the signal has the same properties over time. The overlap time should be similar to the time the signal is changing. Something like the time of Fade Out / Fade In, the transient length.

Example
The following code will recreate the figure form the function (Up to the Colorbar and the units of the Axis):
t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

figure();
spectrogram(x,128,120,128,1e3);

s = spectrogram(x,128,120,128,1e3);

figure();
hA = axes();
imagesc(20 * log10(abs(s).'));
set(hA, 'YDir', 'normal');
```


Calculate moving RMS with exponential averaging

n_d =  pd.DataFrame((data)**2)

np.sqrt(n_d.ewm(alpha = 1).mean())


Detect presence of a specific frequency (or lack thereof)

Regarding the implementation question, maybe you could tell us more about your hardware: ARM family (cortex M3 or M4?), clock rate...

A fixed-point implementation of the Goertzel algorithm is a dozen instructions, so this would indeed be doable with a sample rate of 1 MHz on a low-end ARM micro (say STM32F1 or LPC1343) clocked at 72 MHz (assuming you don't have much overhead for reading your ADC, eg. by letting the DMA subsystem handle this).

However:

You might hit the upper range of your ADC sample rate, and/or will have little margin to design your anti-aliasing filter (say 105kHz is the upper frequency of interest... You can sample at 210kHz but will need the steepest anti-aliasing filter... Or at 800kHz and have more margin for the transition band of your anti-aliasing filter).
This still looks like a waste of CPU cycles, especially if you are only interested in a narrow band of the signal. Keep in mind that if the only thing that interests you is the 99kHz to 101kHz band, you actually need a sample rate of just 4kHz to capture all the useful information - provided you have a mean of modulating the input signal by a 100kHz sine wave.

In such situations, and if you control the target hardware, a bit of pre-processing in the analog domain could be handy, such as:

Frequency shifting your signal by 100kHz so the band of interest becomes the base-band and so that you can sample it at a low sample-rate (requires a sine oscillator, a multiplier, and a low-pass filter).
Applying a narrow band-pass filter in the area of interest + rectifier + low-pass filter - this would measure the energy of the frequency band of interest.


Evaluate expressions without computing DTFT

For the first two questions you just need to use the definition of the DTFT
$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$
and use $\omega=0$ and $\omega=\pi$, respectively.
For 3. you just need to use the definition of the inverse DTFT:
$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$
Finally, for 4. and 5. you need Parseval's theorem. In addition, for 5. you need to express the sequence corresponding to the derivative of $X(\omega)$ in terms of $x[n]$. You should be able to either derive this yourself or look it up in a table.


This is not an answer.  I am just trying to restate the question so that the equation is dimensionally homogeneous:

I don't have an idea how to start it. Help pls :)

$$m(t) = z(t) + \tfrac{1}{4}\operatorname{sinc}\left(\tfrac{t-T}{T}\right)$$

We have some signal $z(t)=\left(\operatorname{sinc}(\tfrac{t}{T})\right)^2$. How often should the $m(t)$ signal be sampled to insure correct signal reconstruction?

this question can be meaningfully answered.

let's define:

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\
1 \qquad & u = 0 \\
\end{cases} $$

$$ \Pi(u) \triangleq \begin{cases}
1 \qquad & |u| < \tfrac12 \\
\tfrac12 \qquad & |u| = \tfrac12 \\
0 \qquad & |u| > \tfrac12 \\
\end{cases} $$

$$ \Lambda(u) \triangleq \begin{cases}
1 - |u| \qquad & |u| \le 1 \\
0 \qquad & |u| > 1 \\
\end{cases} $$


How to get the X and Y Coordinates of the circles from an image?

It is indeed possible to detect the center of the circles by various image processing techniques.

For example, convert to grayscale, apply dilation until the edges disappear, then connected components and compute the centroids.

But this only gives you points, in no specific order, and this might be of little help to you.

It is possible to match the new map against the old one, match the circles, to transfer existing info, but you'll have to check the correct transfer "by hand", anyway.

I wonder why this should be done in JavaScript, as I think that this is an offline operation.


Frequency domain does not really work here. A flanger is NOT a linear time invariant system. You could try to implement this is "piece-wise" LTI and change the delay only at a frame boundary but this will mostly likely result in unacceptable artifacts.

A flanger will require a time variant fractional delay. Rounding the delay the nearest integer delay results in unacceptable "zipper" noise.


Division of complex signal by frequency response

There are several things that could be going wrong here.

Check your FFT lengths

If you do multiplication or division in the frequency domain, your frequency domain data needs to be at least as long as $N+M-$ where $N$ is your signal length, $M$ is your weighting length.

In matlab pseudo code this means you need:

N = length(x);
M = length(w);
freq_length = N + M - 1;

X = fft(x,freq_length);
W = fft(x,freq_length);

y = ifft(X./W);

So your resulting y will be longer than x.

If you don't do this, you will get time-domain aliasing because of the circular nature of convolution using the FFT.

Is your weighting real-valued in the time domain?

You may need to use conj(W) in the above, but you also need to check that w (the time domain version of the frequency weighting) is real-valued. If it's not, you will get a complex-valued signal.

Also note that you may get a small imaginary component.  Often times all that is required to remove this is:

yr = real(y);


Derivative of equation containing Gaussian filtered image

You can represent the convolution as a matrix multiplication with a $N\times N$ Toeplitz matrix, where $N$ is the length of $\mathbf{w}$. Check out this explanation.
So you can write
$$
E(\mathbf{w})=\frac{1}{2}\|\mathbf{AHw}+\mathbf{b}\|^2_2
$$
where $\mathbf{H}$ is your Teoplitz matrix filter. Now if you define
$$
\tilde{\mathbf{A}} = \mathbf{AH}
$$
then
$$
E(\mathbf{w})=\frac{1}{2}\|\tilde{\mathbf{A}}\mathbf{w}+\mathbf{b}\|^2_2
$$

And you can see that the solution is exactly as you wrote but with $\tilde{\mathbf{A}}$ instead of $\mathbf{A}$.


FFT radix-2 in audio - resolution issues - logarithmic scale - sample rate VS buffer size

If you overlap, you can refresh an N long FFT every sample, or every M samples, no need to wait an entire 1 second, or for N samples.

No need to find a sample rate that is 2^n, as most modern FFT libraries support any reasonable product of several small primes.

For analysis, no need to use the same size FFT for both low and high frequencies (or mid, or upper mid, or an entire octave, etc.)

Why FFT?  O(NLogN) provided a speedup that in ancient pre-history allowed DFTs to be run on chips with far far less than a few million transistors, and today still allows real-time performance on power sources that fit in a watch or hearing aid or implant.  Could even allow some Fourier analysis of data when computers were women who used paper and pencil.  The square matrix derivation of a DFT requires less chalk on the classroom chalkboard (today: less iPad scrolling) than more modern audio tricks.

But if you have a larger power source with modern hardware/chips/iPhones, you can use piles of wavelets, large filter banks, psycho-acoustic models, Goertzels, and several other algorithms that are far less efficient than an O(NLogN) sped-up DFT, but together can still easily run faster that real-time requires.


How can I get the polyphase components of a signal

No, that's not right:

Your polyphase components need to be ... well, phases of your signal, and hence, only take every $N$th component. What you do is something different.

To say "these are polyphase components of $h$", you need to define how many components you'll have. You're mixing different $N$s in your code, but overall, your code is just confused – it never selects the $N$th input samples.

Hence, you'd need something like

\begin{align}
p^4_0[n] &= h[4n]\\
p^4_1[n] &= h[4n+1]\\
p^4_2[n] &= h[4n+2]\\
p^4_3[n] &= h[4n+3]
\end{align}

for a $N=4$-phase polyphase decomposition of $h$ (of a specific type; "type" specifies in which order you take the $+0$,$+1$, $+2…$ components and declare them as $p_0$, $p_1$ and so on).

Feels like you need a better text book to start with! In English, the classic would be fred harris' Multirate Signal Processing for Communication Systems. It does an introduction that I'd explain differently these days, but I've worked with such systems a bit, and my initial understanding was guided by that book, so I consider it relatively good.


How to display the output as the input with a fractional delay?

The problem is that you are using "upsample", which merely inserts zeros in between the original samples, instead of "interp", which upsamples and then low-pass filters.


Frequency spectrum of complex equivalent Baseband Signal - Frequency Content around carrier frequency

You have a bug in the formulation of your complex heterodyne oscillator signal inside your complex multiplication.  This seems to be creating both sum and difference frequency results.  Thus the extra high-frequency spectrum.


First of all, it's important to understand that there is no single best way to transform a continuous-time system to a discrete-time system. The method you're using is called backward Euler method, and it is defined by the mapping

$$s\leftarrow\frac{1-z^{-1}}{T}\tag{1}$$

Note that in $(1)$ you scale by $1/T$, where $T$ is the sampling interval (i.e., $1/T$ is the sampling frequency). Eq. $(1)$ just means that continuous-time differentiation (multiplication by $s$) is approximated by a first-order backward difference ($z^{-1}$ corresponds to a delay by $1$ sample). This and other methods for converting continuous-time systems to discrete-time systems are discussed in this answer.

Here is an example where I chose an arbitrary frequency range of interest $f\in[0,10]$, and an appropriate sampling frequency $1/T=20$. The figure below shows the frequency response of the continuous-time system (in blue), and the frequency response of the discrete-time system (in green), which was transformed from the continuous-time system according to $(1)$.


How to modify an LTI differential equation to be acausal or anti-causal?

It's important to realize that generally the differential equation (DE) alone doesn't tell us anything about causality. You claim that the system given in your question is causal. However, an anti-causal system with impulse response $h(t)=-u(-t)$ (where $u(t)$ is the step function) also satisfies the very same DE! So the given DE actually describes two systems, a causal one and an anti-causal one.

What you need to enforce causality (and linearity and time-invariance) is the auxiliary condition that the system is initially at rest. "Initially" means the moment in time when the input signal becomes non-zero. If we denote this moment by $t_0$, these conditions are

$$y(t_0)=y'(t_0)=y''(t_0)=\ldots =0\tag{1}$$

where $y'(t)$ denotes the first derivative of the output signal, etc. Solving the DE with the initial rest conditions $(1)$ for $t>t_0$ gives you the output signal of the causal system described by the DE.

A similar thing is true when using the Laplace transform to compute the transfer function of an LTI system. The expression for the Laplace transform of the transfer function does not uniquely identify the corresponding system. You also have to specify the region of convergence (ROC). In your example, the system's transfer function is

$$H(s)=\frac{1}{s}\tag{2}$$

Depending on the ROC, this function corresponds to two different time domain functions:

$$h(t)=u(t),\quad \text{Re}(s)>0$$
and
$$h(t)=-u(-t),\quad \text{Re}(s)<0$$

You need to choose which solution makes sense in your case.


Audio Processing Used in Voip Mobile Apps

How would a VoIP app developer decide what audio processing should be
  used in an app ?

experience, testing, looking at what others do (aka. typical software development).

So, you, for example, know that you need to work in situations where you'll have echos of your speaker in your microphone, so you add echo cancellation.

Do that with all negative effects that you know beforehand.

Then build a prototype, identify problems, find method to solve problem, repeat.

All in all, engineering.

  Is there a list of all processing options somewhere
  online ?

Sorry, that list would be infinite, so no.

But: You could read a (good) book on speech processing. That would make (somewhat) of an expert on the topic, and experts find solutions to problems by applying their knowledge!

  Is there a methodology involved in selecting the needed
  audio processing for an app ?

See first point. Engineering.

  If I were to identify the minimum set
  of audio processing needed for a VoIP app, how would I go about doing
  it ?

I'd argue that you start with a decoder for the audio codec and call "none" a minimum amount of post-processing that. Then you work up until the application works as intended.


Standard Savitzky-Golay filters are linear phase (type I) FIR filters. So they have an odd number of filter coefficients $2N+1$, and the delay equals $N$.

For a good overview of Savitzky-Golay filters see this article by Ronald Schafer.

For the definition of the four types of linear phase FIR filters see this answer.


In the digital domain a higher sampling will have a better SNR

No, that's not the case.

In fact, if you oversample a bandlimited signal, you don't get overproportionally more signal power than you get more noise power, but as long as you don't use perfect filters which don't exist, you do see more noise, introduced at every step of the analog-to-digital conversion.

You can get an improved SNR from an oversampled signal. But it's not due to higher sampling, but due to processing gain of what you you do with the sampled signal.

  Similarly does a narrow filter give a better SNR in analog domain?

If you assume your noise has a wider bandwidth than your signal, then quite directly follows that the less bandwidth that contains only noise you let through, the better your SNR gets.

Don't forget that SNR is usually just a "helper figure", not inherently something you care about without looking at what has happened to your signal. For example, inevitably, very narrow analog filters have very nonlinear phase, which leads to very strong distortion of the signal's phase. At the end you might have a very good SNR, which only compares powers of signal to noise, but your signal might still be very, very broken.

  Once analog signal is converted into digital signal can SNR be increased? By using a digital filter?

That is one way of achieving processing gain as outlined above, yes. So, no, your higher-sampled signal doesn't have better SNR, your digitally postprocessed signal has.


How many coefficient sections should I select when designing filters with MATLAB FDATOOL?

If you had infinite-precision math, then no it would not change anything to split your order-24 IIR in 3 cascaded order-8 sections. You could even implement it directly as a single difference equation, i.e. no cascade.

However, in the real world, we don't have infinite-precision math, we typically use 64-bit doubles, 32-bit singles or even integers (fixed-point). In those cases, implementing an order-24 IIR as a single difference equation would lead to quantization issues. Your filter would likely be unstable. Even if it were stable, it would probably not have the expected frequency response.

Therefore what we usually do is split the order-2M IIR filter into M cascaded order-2  IIR filter called byquads. The quantization issue is much easier to deal with when using byquads, at least when using doubles or singles. With Fixed-point arithmetic, even when using byquads, you need to double-check everything.

As a bonus, in DSP libraries, you usually have optimized functions for byquads IIR.

Care to explain why you need an order-24 IIR? I fail to see why you would need something more than an order-8.


Use of equalizers in communication systems for Intersymbol Interference

A (linear) equalizer tries to compensate for the (linear) filtering effect of the channel. This filtering introduced by the channel creates intersymbol interference (ISI), so the goal of the equalizer is to reduce (or, ideally, eliminate) ISI, such that the symbol error rate is reduced to an acceptable level. Linear filtering creates ISI by smearing out the symbols over more than one symbol period, making neighboring symbols interfere with each other.

Time-domain and frequency-domain equalizers are just specific implementations of equalizers, their basic function is the same. Furthermore, compensating for the impulse response of the channel is the same as compensating for its transfer function because impulse response and transfer function are just two different descriptions of the same phenomenon: linear distortion.

Unknown channels can be estimated, but usually adaptive equalizers are used, which can adapt to unknown and/or time-varying channel characteristics. Usually, a training sequence of known symbols is used for an initial adaptation of the equalizer. There are also methods that do not require a training sequence; these methods are referred to as blind equalization (e.g. constant modulus algorithm).

The output of the equalizer is indeed an estimate of the transmitted signal; the final decision is done by the slicer following the equalizer.


Direct-Form II is the canonical structure to implement an LCCDE, using minimum number of delay elements and multiplications. The input $x[n]$ is placed on the left, and the output $y[n]$ is on the right, complying with the signal processing conventions.
DF-II structure, implements the poles (IIR part) before the zeros (FIR part), however this may be undesirable for numerical issues with high gain poles.
The transposed DF-II structure, which maintains the canonical advantage of DF-II while avoiding the potential numerical instability of it, is obtained by reversing all branch directions, and replacing sum nodes with signal nodes, also replacing (and renaming) the input output nodes too; hence input becomes at the right side and output emerges at the left node. Transposed DF-II implements zeros first, and poles later.
However, for convenience, to comply with the convention that the input is applied at the left node, the transposed DF-II structure is horizontally flipped. This way, input $x[n]$ appears at the left again.
Both forms of transposed DF-II are equivalent.


HINT:

Use the property

$$z_0^nx[n]\Longleftrightarrow X(z/z_0)$$

with $z_0=1/2$. You can find this property in Table 3.2 on page 126. It's quite straightforward to show this property by using the definition of the $\mathcal{Z}$-transform.


Visualizing the measure of each FFT frame as an image in an audio signal

Very simple:

The Fourier transform of every real-valued signal is (hermitian) symmetric.

You can thus, for your real-valued signals, simply not visualize half of the FFT bins.


Draw the three discrete time signals (shown on the left side) on top of one another. You will see the summation leading to $\delta[n]$ and $\delta[n-1]$ term as follows. For $n=0$, only the first term survives and hence $3\delta[n]$. For the second term, put $n=1$ that makes $3/2+2=7/2$ at index $1$.

Then, obviously third term extends in $u[n-2]$ form as a summation from all three terms at each subsequent index. This is because $n=2$ brings $3/4$ from the first term, $1$ from the second term and $1$ from the third term, totaling to $11/4$ onwards.


What is the difference between object detection and target detection? Are they the same thing?

"Target" is a military term; historically, that's why it's the word used for "identifiable thing from a radar system".
So, I wouldn't use "target detection" outside of radar systems, whereas the word "object detection" is generally used in image processing (for anything that isn't a person). I might be a bit alone with this, but if I look at an image and someone asks me to perform "target detection" on that, I'd assume they want to shoot at that.
I think "target" in radar has a more "physical meaning", as in, it's actually a single contiguous physical object which emits/reflects radio-frequency emissions that should be recognized / tracked as such.
"Object" can be anything, from a round ball, to scattered remnants of a white blood cell; it's slightly less "physical", more "appearance-based".  (Of course, reflecting light is still reflecting EM waves; the difference lies in the number of sensor elements in a camera vs radar array.)


s-Domin or z-Domain - What to Use for Mixed systems

First choice :

Convert the Laplace transform of your process to the Z-domain using the ZOH method as it models your DAC. In your case, your DAC is a PWM.

Second choice :

Work in the Laplace domain but take in account the delay due to your PWM, which should be T/2 where T is your PWM period. To model the delay, you can use the Padé approximation.

I personnaly used the second technique. It was easier, more intuitive. But both techniques will work.

For the numerical issues :
The order of your transfer function might get large when including the controllers, I used minreal too, but you can convert your tranfer functions to ZPK, this solves many of the numerical issues you might have.

Edit : As Tim Wescott mentioned in the comments. The ZOH modelization is accurate when using a real DAC. Not so for a PWM. It's better off to work in the s-domain and try to model the delays (sampling, PWM, computation, etc.) as accurately as possible.


Relationship between sampling and number of pixels?

Yes, that's correct. It's a natural extension of one dimensional up and down Sampling. However, one can define lattices (arbitary shapes) for sampling images, so one could downsample more in a particular axis than the other. Essentially instead of a sampling period, we have a sampling grid. There are benefits of using different freqeuency grids. Of course all of this holds without aliasing only when the sampling theorem holds. Otherwise you get aliasing.


Normalization purpose in signal processing

Normalization is basically bringing the two signals to the same range or a predefined range. A typical example of a predefined range is the statistical perception of the normalization, which is transforming the signal so that its mean is $0$ and standard deviation is $1$. After such transform the canonical form is obtained. In layman's terms, such the operation consists of subtracting the minimum value and dividing by the range (difference between the maximum and minimum values).

Transforming all signals to such canonical form eases and robustifies the process of comparisons as well as serving too different needs such as visualization, and analysis. If we are to speak in terms of image normalization the effect would be the following:

Un-normalized:

Normalized:

Note that, the distinction between the colors is clearer and more apparent.


First pay attention that the Nyquist Shannon Sampling Theorem requires
sampling at rate which is strictly larger than twice (Though it is for edge cases, yet this is the accurate theorem).

Sampling is 1 step to get signals to be Discrete, which means we can store them in a digital form.
Sampling is the discretization of the Time Axis and we need to apply Quantization which is discretization of the Value Axis.

Sampling Theorem states that given a sampled signal (Which is a set of the signal values at discretized grid of time samples) what are the cases it can be fully reconstructed given a single property of the signal - Band Limited Signal.

So what's changed?
You don't have the whole signal but only set of samples from it on a time grid.


Design and a notch filter using two Butterworth filters in parallel (MATLAB)

Cost-effectiveness and ease of implementation of the filter is emphasised.

It that's a your goal then you really shouldn't be using two Butterworth filters. Use a notch filter instead. If there is a "requirement" to use other filters, please explain why this is required. I think your best bet here is to change the requirement.


Windows Audio WAVEFORMAT buffer BYTE to 16-bit int

This doesn't seem to be very good code, and it certainly isn't well commented to explain why it is doing what it is doing.

To answer your question at the "*****  WHAT'S GOING ON HERE  *******" mark:

1) The "* 2" in the chunk.size calculation is clearly because the data is assumed to be stereo.

2) Each sample value is a 2 byte int (16 bit), the chunk needs to have chunk.size samples.

3) The multiplier is to compensate for the memcpy working differently on different platforms.  A char may be one or two bytes,  memcpy moves char sized pieces.  It used to be strictly bytes, but that isn't the case with this code.  This problem would be more efficiently solved with a global DEFINE.

4) The memcpy copies the data from the input buffer to the storage location

5) The loop looks for a silent chunk, allowing a little tolerance -1 to 1 is silent.

Now, here is why I think this is poor code:

For ordinary audio, there are going to be a heck of a lot of Log messages generated.
If it is non-zero (not silent) the chunk will be kept (pushed onto the queue).  If logging every value isn't needed, the loop should break as soon as the chunk is determined not to be silent
The logic of zeroing 1s and -1s in a kept chunk makes no sense to me
If the chunk is not kept, the chunk.chunk array is not deleted (assuming the kept chunks are deleted downstream).  This is a memory leak.

That's enough for now, hope it was helpful.

Ced


I believe that it will be worth for many people to do it from scratch. Let's call the signal after input and feedback summation $E(s)$:

Then question arises what $E(s)$ is equal to? Well that's simple:

$$E(s)=5R(s)-(5+K_vs)C(s) $$

So we have our backward path, now the forward one. What $C(s)$ is equal to?

$$C(s)=E(s)K_p \dfrac{2}{(Js+B)s} $$

Let's substitute previously calculated $E(s)$ to the above equation:

$$C(s)=\dfrac{2}{(Js+B)s}K_p[5R(s)-(5+K_vs)C(s)] $$

Because we have our output $C(s)$ on both sides of the equation, we must eliminate that. Expand the right side of the equation:

$$C(s)=\dfrac{2}{(Js+B)s}5K_p\cdot R(s)-\dfrac{2}{(Js+B)s}K_p(5+K_vs)C(s) $$

Now move the part with $C(s)$ to the left side and move it before parenthesis, this yields:

$$C(s) \left[1+ \dfrac{2}{(Js+B)s}K_p(5+K_vs)\right]=\dfrac{2}{(Js+B)s}5K_p\cdot R(s) $$

Now you can divide right side of the equation by the factor from the left side:

$$C(s) =  \dfrac{\dfrac{2}{(Js+B)s}5K_p}{1+ \dfrac{2}{(Js+B)s}K_p(5+K_vs)}\cdot R(s) $$

Great! Now we have direct relation between the input and the output. Let's make this fraction more pleasant:

$$C(s) =  \dfrac{\dfrac{10K_p}{(Js+B)s}}{1+ \dfrac{10K_p+2K_pK_vs}{(Js+B)s}}\cdot R(s) $$

Now expand the one in denominator to: $\dfrac{(Js+B)s}{(Js+B)s} $, and add it:

$$C(s) =  \dfrac{\dfrac{10K_p}{(Js+B)s}}{\dfrac{(Js+B)s + 10K_p+2K_pK_vs}{(Js+B)s}}\cdot R(s) $$

We can now remove common $(Js+B)s$:

$$C(s) =  \dfrac{10K_p}{(Js+B)s + 10K_p+2K_pK_vs}\cdot R(s) $$

Finally let's do one last multiplication:

$$C(s) =  \dfrac{10K_p}{Js^2+Bs + 10K_p+2K_pK_vs}\cdot R(s) $$

And voila!

$$C(s) =  \dfrac{10K_p}{Js^2+ s(B +2K_pK_v)+ 10K_p}\cdot R(s) $$


I am sorry that I cannot see where your problem is. What do you mean by "how can I save the result?" If you want to use $oim$ in matlab, then it is already there. If you are interested in extracting one specific class, you can simply get the corresponding $n$th mask by calling $mask=(oim==n)$, which returns a binary mask with only the pixel belonging to the $n$th class marked as 1s.

If you are going to save the classification result outside of matlab, for example a tif image, simply use imwrite function in matlab. Late you can read it via imread function.


Recommendation for courses / studies on digital signal processing

I've kind of grouped your subjects into larger overall subjects.
Note that there's a lot of overlap here, with the possible exception of actually making it work in a microprocessor (except -- in my opinion the best person to implement something is someone who understands it.  So -- overlap).
Specifically, you could claim that it's all applied math.  Or all signal processing.  You could even stretch this list a bit and wedge it into control theory, but it wouldn't be an easy fit.

Applied Math

Singular Value Decomposition
Principal Component Analysis
Machine Learning for Audio Processing

Signal Processing

Digital Filter Design
Kalman Filters (especially in the context of inertia measurement units and nonlinear problems)
Wavelet Transforms

Control Theory

Kalman Filters (especially in the context of inertia measurement units and nonlinear problems)
Control theory (implementing a control algorithm and tuning it for a problem)

Embedded Software (or systems) Engineering

Running a control algorithm on a microprocessor

A weakness of the US tertiary educational system is that all of these subjects will be used in engineering departments, but they're all really applied math.  So different universities will group them differently, and it can be hard to find a program that will cover them all.  I'm not sure if that applies to where you live.  I happened to end up -- by luck more than design -- at an institution that let me get enough of a grounding in the related subjects that I could learn the rest on my own (Worcester Polytechnic Institute -- it's a great place).
I would suggest that you approach this by finding a university program that includes the classes you need to get a grounding in this, and get your second master's.  Don't expect to learn all of this within a Master's program -- I can lay claim to about 90% of this list, but about half of it are things that I learned on the job, over the last three decades.
Classes to look for are below.  Note that you're not going to have time for them all, so you'll need to pick and choose.  If you want to grow into your full list, take as much math as you can and expect to do self-study:

Stochastic signals and systems (basically, this is the study of finding the optimal system to process a signal with random components).

If you haven't taken it, a 4th-year course on statistics from the math department.

Real analysis (AKA "Advanced Calculus", but I think that's an obsolete name).  It's an invaluable help when you're sitting at your desk at work wondering if the nifty new thing you're applying is actually mathematically sound.
If it's offered, estimation and detection theory.  This is basically a follow-on to stochastic signals and systems.  It's a deep dive into the math underlying the detection of signals, signals in noise, etc.  The one that I took had you derive the basic Kalman filter from first principles in a homework problem -- which gives you an idea of the depth of the subject.
Any class on optimal state estimation (i.e., Kalman filtering and all its variations).
Signal processing, assuming you haven't taken a class in it.
Digital signal processing, ditto.
Control theory classes will just be called "control theory", but you may not be able to fit much in.  If there's a class offered in state-space systems and you've already got classical (transfer function) control under your belt -- take it.
If you want to specialize in control theory, go for nonlinear control or whatever is offered.  But if you do that, you probably won't have time for the Detection & Estimation class, so think hard.

I'm not sure what to suggest on the implementation front -- that's basically yet another specialty.  I just picked that up on my own, on the way.  But, I started out as a kit with a hobby in electronics and a microprocessor board, so I was already doing basic embedded programming when I was 13.  Then when I got an EE degree, any class that involved digital circuits was either "oh, yes, I've done this", or "oooh!  what a neat formal way to do what I've struggled through intuitively".


How to change volume of a PCM 16 bit signed audio?

Changing the volume of an audio signal must be done by applying a gain (multiplication) - and optionally clipping if your system has a limited dynamic range. This is as simple as that. Applying a non-linear function to an audio signal will cause distortion and add harmonics, and you don't want this to happen - you want to modify the loudness of the signal, not its timbre. [To be fair, there are non-linear processings designed to change the perceived loudness of the signal without affecting the timbre, within a given dynamic range constraint (eg. multiband compression), but it doesn't look like this is what you need.]

Where non-linearity and fancy response curves come to play is when designing a user interface - when deciding on the relationship between the position of the control (knob or slider, whether on a GUI or as physical hardware) and the gain applied to the signal. This is where perception matters, because the users will expect a mapping between the position of the slider and their perception of loudness. Please note that even if the relationship between the position of the volume control and the gain applied to the signal is non-linear, the process of applying the gain to the signal is linear, and non-linearity would be unwanted there!

When it comes to physical volume controls, eg in hifi systems or personal audio players, the relationship between the knob position and the attenuation is closer to an exponential curve, though its shape has been tweaked and is constrained by the manufacturing process - sometimes it's just two or three linear segments. You can find those curves in the datasheets from manufacturers ("A" taper). Mixing console faders usually have their response compressed so that the upper half of their travel covers the useful range of -20 dB..+6dB.

In the software world - at least for music production - it is most common to have volume/gain knobs calibrated in dB. For example, if you have a 100 pixels long volume slider graduated from -48dB to +6dB, the gain applied to the signal would be $10^{\frac{-48 + 54 \frac{x}{100}}{20}}$.


Applying DFT twice does not actually reverse an array. Instead, the first element stays in place while the rest of the array is reversed. Why?

The sequence is exactly what you should expect:
$$x[-n]=x[N-n]\tag{1}$$
Clearly, for $n=0$ $x[n]$ and $x[-n]$ have the same value.
It seems like you were expecting to see the sequence $x[N-1-n]$ instead of $x[N-n]$. If you wanted that sequence you would need to modulate the result of the first DFT:

N = 16;
x = 1:N;
c = exp(-1i*2*pi/N*(0:N-1));
x2 = real( fft( ( fft(x) .* c ) ) ) / N;    % real() just to remove rounding errors
x2 =

  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1


Low Level Audio Engine Audio Data Reformatting - Sample Rate and Bit Depth Conversion

Sinc interpolation introduces no aliasing artifacts to data that was originally properly bandlimited before sampling.  (But needs to be slightly modified, as a Sinc or Dirichlet kernel has an infinite duration.)
One canonical reference on why reconstruction works and how to interpolate is here: https://ccrma.stanford.edu/~jos/resample/
I have a simple implementation in Basic here: http://www.nicholson.com/rhn/dsp.html#3
The other common method is low-pass filtered integer ratio upsampling followed by low-pass filtered integer ratio downsampling by the appropriate two ratios.  Another common method is to generate a polyphase FIR filter bank (of multiple phases of windowed Sinc (or otherwise optimized FIR filters) interpolation kernels.  Another common method is to use polynomial interpolaters of each of the Sinc kernel's humps, resulting in a Farrow filter bank.
You need to carry enough bits of arithmetic intermediate values in order to keep your quantization noise floor below whatever your spec might be.


Breaking down the conversion into an upsampling and a downsampling step as explained in mbaz's answer is not necessary - band-limited interpolation achieves the same result in one single process, for any conversion ratio, and its quality is easy to fine-tune (through various approximations of the band-limited pulse).

libsamplerate is an example of open-source implementation that does not rely on upsampling/downsampling. It should be possible to bind it to NodeJS through an add-on, or even to rewrite it in javascript - once stripped of the code handling corner cases (variable rate, compensation of processing delay), the code fits in a 50-line loop.


Deconvolving a 1d Signal Using a Lookup Table of Kernels

The way I understand the problem is each sample of the output is a linear combination of the samples of the input.
Hence it is modeled by:
$$ \boldsymbol{y} = H \boldsymbol{x} $$
Where the $ i $ -th row of $ H $ is basically the instantaneous kernel of the $ i $ -th sample of $ \boldsymbol{y} $.
The problem above is highly ill poised.
In the classic convolution case we know the operator matrix, $ H $, has a special form (Excluding the borders) - Circulant Matrix. With some other assumptions (Priors) one could solve this ill poised problem to some degree.
Even in the case of Spatially Variant Kernels in Image Processing, usually, some form is assumed (Usually being block circulant matrix, and the number of samples of each kernel is larger than the number of samples in the support of the kernel).
Unless you add some assumptions and knowledge into your model the solution will be Garbage In & Garbage Out:
numInputSamples = 12;
numOutputSamples = 10;

mH = rand(numOutputSamples, numInputSamples);
mH = mH ./ sum(mH, 2); %<! Assuming LPF with no DC change

vX = randn(numInputSamples, 1);

vY = mH * vX;

mHEst = vY / vX;

See the above code. You will always have a perfect solution yet it will have nothing to do with mH.
Now, if I get it right, you say I don't know $ H $ perfectly, but what I have is a pre defined options.
So let's say we have a matrix $ P \in \mathbb{R}^{k \times n} $ which in each row has a pre defined combination:
$$ H = R P $$
Where $ R $ is basically a row selector matrix, namely it has a single element with value $ 1 $ in each row and the rest is zero.
Something like:
mP = [1, 2, 3; 4, 5, 6; 7, 8, 9];
mH = [1, 2, 3; 7, 8, 9; 7, 8, 9; 4, 5, 6; 4, 5, 6];

% mH = mR * mP;

mR = mH / mP;

So our model is:
$$\begin{aligned}
\arg \min_{R, \boldsymbol{x}} \quad & \frac{1}{2} {\left\| R P \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} \\
\text{subject to} \quad & R \boldsymbol{1} = \boldsymbol{1} \\
& {R}_{i, j} \geq 0 \quad \forall i, j \\
\end{aligned}$$
It is still exceptionally hard (Non convex) problem but with some more knowledge it can be solved by utilizing alternating methods where we break the optimization problem as:

Set $ \hat{\boldsymbol{x}}^{0} $.
Solve $ \hat{R}^{k + 1} = \arg \min_{R} \frac{1}{2} {\left\| R P \hat{\boldsymbol{x}}^{k} - \boldsymbol{y} \right\|}_{2}^{2} $ subject to $ R \boldsymbol{1} = \boldsymbol{1}, \; {R}_{i, j} \geq 0 \; \forall i, j $.
Solve $ \hat{\boldsymbol{x}}^{k + 1} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \hat{R}^{k + 1} P x - \boldsymbol{y} \right\|}_{2}^{2} $.
Check for convergence, if not go to (2).

Now each sub problem is convex and easy to solve.
Yet I still recommend you add better assumptions / priors.
Such as the minimal number of contiguous samples that must have the same PSF (Similar to 2D in images where we say each smooth area is smoothed by a single PSF).
Remark
We didn't employ the fact each element in $ R $ is either 0 or 1 as the straight forward use of that will create a Non Convex sub problem.
In case the number of PSF's is small we can use MIP solvers. But the model above assumed each row is a PSF so for large number of samples even in case we have small number of PSF'w the matrix is actually built by shifting each PSF as well. So we'll have large number in any case.
Another trick might be something like Solving Unconstrained 0-1 Polynomial Programs Through Quadratic Convex Reformulation.
Yet the simplest method would be "projecting" $ R $ into the space (Which is not convex, hence projection is not well defined). One method might be setting the largest value to 1 and zero the rest.
Update
In comments you made it clear you know the kernel per output sample.
Hence the model is simpler:
$$ \boldsymbol{y} = A \boldsymbol{x} + \boldsymbol{n} $$
The least squares solution is simply $ \boldsymbol{x} = {H}^{-1} \boldsymbol{y} $.
For better conditioning and noise regularization (Actually prior about the data, but that's for another day) you can solve:
$$ \hat{\boldsymbol{x}} = {\left( {A}^{T} A + \lambda I \right)}^{-1} {A}^{T} \boldsymbol{y} $$
This is a MATLAB code for proof of concept:
load('psfs.mat');

mA = psfs;
vY = y;
vX = x;

vParamLambda    = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1];
numParams       = length(vParamLambda);

vValMse = zeros(numParams, 1);

mAA = mA.' * mA;
vAy = mA.' * vY;
mI  = eye(size(mA));

for ii = 1:numParams
    paramLambda = vParamLambda(ii);
    vEstX = (mAA + paramLambda * mI) \ vAy;

    vValMse(ii) = mean((vEstX(:) - vX(:)) .^ 2);
end

figure();
hL = plot(vParamLambda, 10 * log10(vValMse));
xlabel('Value of \lambda');
ylabel('MSE [dB]');

This is the result:


What's the Intuitive Description of Circular Symmetric Complex Zero Mean White Gaussian Noise?

To expand on Dilip's comment:

The noise being white implies only a certain property of the noise itself, not its correlation with any other signal.
In digital communications, the independence and/or lack of correlation between the signal and the noise is assumed, if not explicitly stated. The reason is that this is what is observed in practice. Even if you can mathematically conclude that a given message signal and a given noise are independent and uncorrelated, there is not much point in doing so.
What is meant by circularly symmetric Gaussian noise is that the noise looks the same in all directions. Imagine a certain Gaussian complex noise that (for some physical reason) produces larger noise in specific directions; for example, along a line at an angle of $\pi$/4 on the complex plane. You could then design a QAM constellation (and detector) that takes advantage of this behavior to get better than expected performance given the noise power in the channel. When a problem states that the noise is circularly symmetric, it's saying that you can't play this kind of game.


Is it possible to calculate/estimate the storage space a certain duration of .mp4 video files will take up?

No.

The MPEG4 video codec allows for a large ranges (easily a dynamic range of 50 for "useful" quality) of compression ratios for the same raw video material. The same applies to audio compressors.

And of course, there's simply video material that is easy to compress, so using the same type of quantization, motion prediction and frame structure, some videos will compress way better than other.

In fact, in many cases, you ask your encoder to aim for a target average bitrate even, so that you don't have to estimate anything, but set the amount of bits per second. If you have a technical necessity to know the bitrate beforehand, that would be the way to go!


Linear and Circular Convolution in Fourier Domain (DFT)

the following matlab/octave code gives the linear convolution result using frequency domain :

A = ((-1).^[0:79]').*hamming(80);    % input one
B = blackman(100);   % input two

C1 = conv(A,B);     % A * B (convolution) in time domain
C2 = real( ifft( fft(A,179).*fft(B,179) ) ); % convolution using freq domain

The output will be identical of length 179 samples:


Opencv Android Perspective Transformation not working as expected

@Override
public boolean onTouchEvent(MotionEvent event){
    boolean value = super.onTouchEvent(event);
    switch ( event.getAction()){
        case MotionEvent.ACTION_DOWN: {
            return true;
        }
        case MotionEvent.ACTION_MOVE:{
            //Point 1
            float x = event.getX();
            float y = event.getY();
            double dx = Math.pow(x - mCircleX,2);
            double dy = Math.pow(y-mCircleY,2);
            if (dx + dy < Math.pow(mCircleRadius,2)){
                //Toucher
                mCircleX = x;
                //Toast.makeText(getContext(),String.valueOf(mCircleX),Toast.LENGTH_SHORT).show();
                mCircleY = y;
                postInvalidate();
                return  true;

            }
            //Point 2
            float x2 = event.getX();
            float y2 = event.getY();
            double dx2 = Math.pow(x2 - mCircleX2,2);
            double dy2 = Math.pow(y2 - mCircleY2,2);
            if (dx2 + dy2 < Math.pow(mCircleRadius2,2)){
                //Toucher
                mCircleX2 = x2;
                mCircleY2 = y2;
                postInvalidate();
                return  true;

            }
            //Point 3
            float x3 = event.getX();
            float y3 = event.getY();
            double dx3 = Math.pow(x3 - mCircleX3,2);
            double dy3 = Math.pow(y3 - mCircleY3,2);
            if (dx3 + dy3 < Math.pow(mCircleRadius3,2)){
                //Toucher
                mCircleX3 = x3;
                mCircleY3 = y3;
                postInvalidate();
                return  true;

            }
            //Point 4
            float x4 = event.getX();
            float y4 = event.getY();
            double dx4 = Math.pow(x4 - mCircleX4,2);
            double dy4 = Math.pow(y4 - mCircleY4,2);
            if (dx4 + dy4 < Math.pow(mCircleRadius4,2)){
                //Toucher
                mCircleX4 = x4;
                mCircleY4 = y4;
                postInvalidate();
                return  true;

            }
            return true;

        }
    }
    return value;
}

}


Why is oversampling helpful to an anti-alias filter?

You might remember Nyquist's theorem. Given a signal which is band limited to $f_1$, we must sample it at least at $2f_1$:

$f_S>2f_1$

So if you check you favourite Signals and Systems book (e.g. Oppenheim's), you might recall that, once sampled, we can consider the signal's discrete Spectrum (which is periodic every $2\pi$ radians, by the way).

(In the image, Hz are considered instead of radians, being a 2*pi scale factor the only difference)

Then, we would go back to continuous time by low pass filtering with $f_{cutoff}=f_S/2$. So let's look into a few cases:

if $f_S<f_1$, the original spectrum and the subsequent replicas would overlap ($f_1>f_S-f_1$). This would result in distortion over the original signal after the low pass filter. This is the so called aliasing.
if $f_S>f_1$, as in the image, there is no overlapping, and therefore no aliasing. Furthermore, we can recover the original signal by applying a low pass filter. Fortunately, we don't need the sharpest transition band in our filter.

(please excuse my lack of accuracy here, the filter would also be periodic as we still are in discrete frequency domain. Also, please take into account that 3 dB cutoff should be on $f_S/2$)

Last, imagine how this would be if we adjusted $f_S=f_1$. Both spectra would converge at $f_S=f_1$, so we would need an ideal (infinite slope) filter to recover the undistorted original signal. Ideal is no good, you know.

Now, let's think about those LP filters.
 (wikipedia)

The wider the transition band is, the flatter the pass band is. On the other hand, the pass band is curlier when the transition band is tighter. Non-flat frequency response means distortion. Don't forget that the signal that we are "deleting" in the stop band would produce distortion, both in the anti-alias and the D/A filter.

So we want to keep the sampling frequency as low as possible, while not distorting the signal, either by the filter ripple or by aliasing. But, on the other hand, we want to keep the filter as cheap and accurate as possible. By allowing a wider transition band, we avoid the need of incrementing the order of the filters, i.e., the number of filter coefficients (in the D/A filter) or of electronic devices (in the anti-alias).

So we confront a cost/distortion trade-off, in which it seems preferable to increase sampling rate (e.g. CD audio is sampled at 44.1 kHz while we hear no more than 20 kHz), so we can keep reasonable distortion and avoid huge order filters.

Take into account that increasing filter's order also leads to an increment in number of operations, although it is less than the one produced by increasing sampling rate.

Filter terminology here, just in case: http://www.dspguru.com/dsp/reference/filter-terminology


How to do image reconstruction from array data?

You will not be able to reconstruct the image without the phase information.  The frequency magnitudes alone are not sufficient.


What you're describing is call deconvolution. It's a concept frequently used for distorted images and for equalization in communications. As such you should be able to find a number of resources for your specific application. In general, the original image may not be recoverable exactly (it often isn't).

However, you can do a pretty good job depending on how severe the image has been undesirably filtered. In addition, many methods will attempt deconvolution by estimating the filter which it appears the signal has been affected by. If the assumption that the signal was convolved with a linear filter is a poor one, your performance may be severely inhibited.


Study resource for Digital communication

I highly recommend you the book Schaum's outlines of Analog and Digital Communications.

It will be the fastest refresher for the signal processing and probability concepts which are essential prerequisites for digital communications. It's not using that heavy math also, nevertheless, not a trivial one either ;-)


Would a visible light beamforming "telescope" be feasible?

Not in the way you think, no.

Digital Beamforming works because the digital receiver is able to sense the phase of the incident radio wave; that phase is retained.

Photodiodes don't sense the phase – they only sense the intensity.

Also, this requires the whole incident radio or light wave to be phase-coherent across its front – that's the case for laser light, but not for two random photons emitted from a huge fusion fireball, I'm afraid. So, the second requirement also falls flat. Generally, even non-laser sources of light have a coherency period – but the more different wavelengths contribute the smaller that gets.

In TX, that principle works – if you feed a laser through an array of digital phase modulators (optical crystals with electrically adjustable refractive indices), you can build phase-coherent multiemitters.

In RX, that might also work, if the received coherent light can be made to fall through a microscope objective, and then onto an acousto-optic deflector, you might, via modulating that deflector's state, be able to build a selective/phase-shifting summer, leading to beamforming. But: that's particle physics sized electronics normally, i.e. suitable for situations where you have a very strong laser illuminating very small points (in fact, I learned about AODs just today, and they're usually used the other way around: take one beam of laser, split it into multiple parallel beams, focus these, and trap quantums in the focus point. Cool stuff.).

Anyway, there's also mathematical accuracy limits to digital beamforming, and I'd must admit that I don't know whether they'd even allow higher precision with the kind of timing accuracy we can do today. Then again, phyisicists invent the coolest devices, and maybe you can


Find beginning of signal for BPSK Demodulation

This problem is typically solved by preceding your actual payload data with a synchronization pattern, or sync word. This is a sequence of symbols that are known to the receiver ahead of time and are used to aid synchronization with the payload data frame. It could be detected in a couple ways:

Sync pattern matched filtering: When looking for the synchronization marker, instead of just integrating over one data bit, you could create a matched filter that represents the entire modulated sync word. This has the advantage that the SNR at the output of the sequence matched filter is higher than you would have otherwise, by a factor of the number of bits in the pattern. This can help if you're trying to detect the frame start in a very low-SNR environment.

Disadvantages of this approach are that it is more computationally complex (because you're correlating against an entire sequence of bits), and also more sensitive to frequency offset between the transmitter and receiver.
Hard decision pattern matching: Instead of the "soft-decoding" approach of correlating against the entire modulated sync sequence, you could instead just look for the pattern's sequence of symbols at the output of your existing receiver. As an example, you could decide that all of your frames will begin with the sequence 01001100101010 (I make no warranty that this is a good choice). You would then slide through the bits that come out of your BPSK decoder and declare the start of a frame when you see the sync pattern. It may also help to allow some small number of bit errors in the pattern when looking for the sequence.

Compared to coherent integration across the entire sync sequence, this technique will not perform as well in low-SNR environments. However, it is much simpler and may work just fine for your application.