It's not completely clear from your question, but I interpret your situation as:

You have a set of OFDM symbols that have been passed through a multipath channel.
You have generated an equalization matrix of phase shifts (in radians) that you want to apply to each subcarrier of each symbol.
You're not sure how to apply the phase shifts to the OFDM symbols.

Being able to implement an equalizer via simple multiplication is one of the benefits of OFDM, since one can assume that fading is flat across each subcarrier's channel. To apply a phase shift of $\theta$ radians to a complex number $x = a + jb$, do the following multiplication:

$$
\tilde x = e^{j\theta}x
$$

$$
\begin{align*}
\tilde x &= \left(\cos(\theta) + j \sin(\theta)\right) \left (a + jb\right) \\
&= \left(a \cos(\theta) - b \sin(\theta)\right) + j \left(a \sin(\theta) + b \cos(\theta) \right)
\end{align*}
$$


Approximate a System Frequency Response with a Filter in MATLAB

Convolve the input signal with (Real Part) harmonic signal $ \alpha {e}^{-2 \pi i {f}_{0} t + \phi} $.
Let's call this output as $ {y}_{1} $.
The parameters $ \alpha $ and $ \phi $ are from your reference system response.

This will give you the output of the reference system on your input signal only at your desired frequency.

Now convolve the input signal with (Real Part) harmonic signal $ 1 {e}^{-2 \pi i {f}_{0} t} $.
Let's call this output as $ {y}_{2} $.
This will extract this component of your input signal.

Calculate $ {y}_{3} = {y}_{1} - {y}_{2} $ and all you need is to add $ {y}_{3}  $ to your signal.

There are trickier and more efficient way to do so.
Yet I wanted to keep it simple and intuitive.


TL;DR:

CWT's $\partial_t$ is taken with respect to wavelet timeshifts, STFT's with respect to window time samples, with CWT's FT taken along rows, and STFT's along columns.
The $2\pi$ is per assuming a different 'oscillatory factor' of Fourier Transform's kernel, i.e. $e^{-j2\pi f}$ vs. $e^{-j\omega}$. Numeric implementations can work with either.
"Why modified STFT" is a separate question, addressed in second answer.

Motivating phase transform
To derive CWT's phase transform, original paper uses $f(t) = A \cos(\Omega t)$ as a motivating example, or for our purposes equivalently $A e^{j\Omega t}$, then works backwards from the result of $\partial_t W_f / W_f$ to define $w$ such that $\Omega$ is extracted. The idea is, this $f(t)$ is the basis function (building block, atomic unit) of the Fourier Transform, so results here will transfer to  any other $h(t)$ the more $h(t)$'s (D)FT is "well-behaved" (more specifically, $h(t)$'s "components" obey synchrosqueezing assumptions detailed in followup paper), which spans a large variety of real-world signals.

Deriving phase transform
The SSQ STFT paper considers $e^{j\omega t}$ directly, but omits a derivation; I will derive it here for both the 'standard' and 'modified' STFT, along their derivatives, from which the phase transform will arise naturally.
Let $f(t) = e^{j\xi_0 t}$, and assume a real, even-symmetric ($g=g^*,\ \hat{g} = \hat{g^*}$), window $g(t)$. We have:
$$
\begin{align}
S_gf(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi t} dt \\
& = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi t} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju\omega} d\omega \\
& = \hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\
\end{align}
$$
where Plancherel-Parseval's Theorem was applied, and similarly
$$
\begin{align}
S^{\text{mod}}_g f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi (t-u)} dt \\
& = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi (t-u)} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\
& = \hat{g}(\xi - \xi_0) e^{ju\xi_0} \\
\end{align}
$$
The derivatives are:
$$
\begin{align}
S_{g'} f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi t} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju\omega} d\omega \\
& = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\
\end{align}
$$
$$
\begin{align}
S_{g'}^{\text{mod}} f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi (t-u)} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\
& = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{ju\xi_0} \\
\end{align}
$$
Then, for both cases, to retrieve $\xi_0$, we consider the ratio
$$
\begin{align}
\frac{\partial/\partial t( S_{g}f(u, \xi) )}{S_gf(u, \xi)} = \frac{S_{g'}f(u, \xi)}{S_gf(u, \xi)} = j(\xi - \xi_0)
\end{align}
$$
which leads us to define
$$
\begin{align}
& \boxed{ w(u, \xi) = \xi - \frac{S_{g'}f(u, \xi)}{jS_gf(u, \xi)} } \\
&= \xi - \frac{j(\xi - \xi_0)}{j} = \xi_0 \\
\end{align}
$$
Note the phase transform being same for both is one consequence of the division "reducing influence of windowing".

Verifying derivation
$S_g f(u, \xi)$ expression can be verified with an example, comparing against Wolfram's output (set normalization=1, oscillatory factor=-1):
$$
\text{Let}\ f(t) = e^{jt},\ g(t) = e^{-\pi t^2}, \text{ and thus}\ \hat{g}(\xi) = e^{-\xi^2 / (4\pi)}.\ \text{Then}, \\
$$
$$
\begin{align}
S_gf(u, \xi)
&= \hat{g}(\xi - 1) e^{-ju(\xi - 1)} \\
&= \exp{\left( -\frac{1}{4\pi}(\xi - 1)^2 - ju(\xi - 1) \right)} \\
&= \exp{\left(\frac{1}{4\pi}( (\xi - 1)(-j4\pi u - \xi + 1 )\right)} \\
\end{align}
$$
More completely, ssqueezepy's ssq_stft on a variety of signals yields expected behavior.

Why the $2\pi$?
The TL;DR pretty much says it, which can be verified in the original MATLAB Toolbox's paper. In code, with frequency-domain differentiation, the $2\pi$ makes it so we define $j \omega$ to span $[0, ..., \pi]$ and $(-\pi, ..., 0)$, and without uses $0.5$ instead of $\pi$, reinterpreting the derivative's dimensionality as linear cyclic fraction of sampling frequency, rather than radian. See code.

Why the $\xi\ -$?
CWT's $\partial W_f$ takes $\partial / \partial b (\psi((t-b)/a))$, where $b = $ timeshift, whereas STFT's takes $\partial / \partial t (g(t))$. This is per how each is computed:

CWT takes DFT along rows; for every scale (row), we use frequency-domain convolution to correlate the wavelet with the signal at every timeshift $b$ (column).
STFT takes DFT along columns; for every timeshift (column), we apply DFT to correlate the windowed signal with cisoid bases at every frequency (row).

Following the derivations from motivating examples reveals the $\xi \ -$ stemming from this difference.
The definitions can be interchanged or made the same for both; why exactly they differ I'm unsure. My bet is, CWT is best computed for every timeshift (hop_length = 1), and the scales are already sparse, whereas in STFT we often apply large hop_length (small overlap), so convolution theorem is very compute-wasteful.


LMS adaptive filter relatively delayed signal and reference inputs

The adaptive filter tries to emulate the assumed filtering process between the noise reference signal and the actual noise in the noisy signal.

If $n(t)$ is the actual noise in the noisy signal, and $n_r(t)$ is the noise reference, it is assumed that there's a linear filtering relationship between the two:

$$n(t)=(n_r*h)(t)$$

where $*$ denotes convolution (filtering), and $h(t)$ is some unknown impulse response, which should be approximated by the adaptive filter.

The practical effect of a delay of the noise reference depends on the impulse response $h(t)$. A delay of $\tau$ in the noise reference means that the adaptive filter needs to model $h(t+\tau)$ instead of $h(t)$. This can be a good thing if $h(t)$ is almost zero in the interval $0<t<\tau$, because then optimal use is made of the length of the adaptive filter. However, if there is significant energy in $h(t)$ in that interval, then, after delaying the noise reference, the adaptive filter must approximate a non-causal filter, which it will have a very hard time to do.

So in general, a delay of the reference signal results in a shifting of the impulse response that needs to be approximated by the adaptive filter. If this shifting is helpful or not depends on the amount of shift and the properties of the impulse response.


How is the number of amplitude levels decided in modulation techniques?

There are $M$ symbols in an $M$-QAM constellation where two carriers I and Q are amplitude modulated, and there are $N=\sqrt M$ amplitudes for each I and Q carriers (for $M\ge4$). In $16$-QAM there are 16 symbols on a $4\times4$ constellation. That is, there are $4$ amplitude levels for the I and $4$ levels for the Q component. Of course, $\log_2M$ bits are conveyed by each symbol. Since $M$-QAM normally uses carriers with symmetric amplitudes, $M$ is usually a power of $4$. For instance, $M$ can be $4, 16, 64, 256,\cdots$.


These are very computationally expensive problems and it also depends if your polygons are convex. Also in your case you don't have any extra constraints and it can be any line within. I think that the best read for you is the following publication:

  Hall-Holt O., et al. - Finding large sticks and potatoes in polygons

You will definitely find few algorithms among ones they described (computational complexities are given). Here is some brief presentation concerning their work: click


EM-algorithm jointly estimates the state space model parameters as well as state estimates. The E-step is a Kalman filter, which uses the current estimates to predict the new states. The M-step uses this result in an MLE procedure to obtain the parameter estimates.

I would suggest that you read the following paper, which explains a generative model, from which all of the tools such as HMM, Kalman Filter, VQ and etc. are derived:

  Read Roweis, Ghahramani, "A Unifying Review of Linear
  Gaussian Models", Neural Computation, Vol. 11, No. 2, 1999.


The x axis is simply an index number, and this is showing your the correlation for leading and trailing shifts. The shift at index 80 is lag = 0, the origin you are interested in. Type help(xcorr) in Matlab to get further detail in how to return the lag. Following that, if you do the following you will get the result you are looking for:

r, lag = xcorr(frame, frame);
figure()
plot(lag, r)

I did it as shown above to highlight how to return lag from xcorr specifically while minimizing change to your code, but xcorr(frame) will return the same result as xcorr(frame,frame) so can be further simplified for the case of autocorrelation.


MATLAB: Issue with passing data through AWGN Channel

UPDATE: Scaling was not adjust for upsample/downsample.

data = randi([0 3], 1, 2^12); % Generate data
dataMod = qammod(data, 4);  % Modulate data
dataMod_u = 16*upsample(dataMod,256); % Upsample
b = rcosdesign(0.25,8,256,'sqrt');
dataMod_f = filter(b, 1, dataMod_u); % Tx Filter

channel2 = comm.AWGNChannel('EbNo',15,'BitsPerSymbol',2,'SamplesPerSymbol',256);
ch_out2 = step(channel2, dataMod_f); % Pass through channel
ch_out_f2 = filter(b,1,ch_out2);    % RX Filter
ch_out_f2 = (1/16)*downsample(ch_out_f2,256); % Downsample
scatterplot(ch_out_f2)


The standard approach would be to

Choose your feature vectors.
Partition your data set into training and testing.
Calculate the feature vectors for all your data (both training and test).
Use the training data set to train your classification algorithm.
Apply the trained classification algorithm to the test data set and calculate your correct classifications and incorrect classifications to get your precision / recall.

The trick will be step 1: choosing your feature vectors, given that you've only got Fourier coefficients (and values derived from them).


What is the correct length for obtaining a true linear convolution from DFT?

I think what you are asking is: what happens of I implement the convolution of two sequences of length M & N, using spectral multiplication with an FFT of length K? We defined "correct" here as "matching the result of a linear convolution". For simplicity we define the correct length $L = M + N -1$
Here is what you get

K > L: The first L samples will be "correct", followed by (K-L) zeros
K = L: Correct
max(M,N) < K < L: the first (L-K) samples of the result are wrong because of time domain aliasing. You are also missing (L-K) samples at the end or your result since your buffer is too short. Some samples towards the end of your result to match the correct answer somewhere the middle.
K < max(M,N): everything is wrong


How to check if Fourier components are in phase of out of phase?

Depends on the phase reference point in time.

Depending on the reference point, every DFT basis sinusoid can be either in-phase or out-of-phase (by 180 degrees) with every other basis, thus the same with any coefficient weighting of these basis vectors.

If you reference phase to the center of a DFT (using an fftshift), then the angle will tell you whether a particular component has even or odd symmetry relative to the center.


Is it appropriate to use Sobel operators to find the derivatives of an image?

Indeed, it adds smoothing in the $y$ direction. The Sobel filter is the separable combination of the centered derivative $[−1,\;0,\;1]$ along $x$, and the $3$-point binomial smoother $[1,\;2,\;1]$ along $y$.


Prove using superposition that an LTI system is linear

I'm going to assume the system has the initial condition $x(0) = 0$, as otherwise the system isn't linear.

Suppose the input $u_1(t)$ yields the output $x_1(t)$ and the input $u_2(t)$ yields the output $x_2(t)$.

We need to prove that the input $u(t) = \alpha u_1(t)+\beta u_2(t)$ yields the output $x(t) = \alpha x_1(t) + \beta x_2(t)$ for any scalars $\alpha$ and $\beta$.

The given input and output pairs must satisfy the differential equation and the initial condition. Hence, $x_1'(t) = Ax_1(t)+Bu_1(t)$ and $x_2'(t) = Ax_2(t)+Bu_2(t)$ for all times $t$ and $x_1(0) = x_2(0) = 0$.

Now, add $\alpha$ times the first equation to $\beta$ times the second equation to get: $\alpha x_1'(t) + \beta x_2'(t) = \alpha(Ax_1(t)+Bu_1(t)) + \beta(Ax_2(t)+Bu_2(t))$

Then, rearrange stuff to get:

$\dfrac{d}{dt}[\alpha x_1(t) + \beta x_2(t)] = A(\alpha x_1(t)+\beta x_2(t)) + B(\alpha u_1(t)+\beta u_2(t))$

$x'(t) = Ax(t)+Bu(t)$

Also, we have $x(0) = \alpha x_1(0)+\beta x_2(0) = \alpha \cdot 0 + \beta \cdot 0 = 0$.

Thus, the signals $u(t) = \alpha u_1(t)+\beta u_2(t)$ and $x(t) = \alpha x_1(t) + \beta x_2(t)$ satisfy the differential equation and the initial condition. Hence, the input $u(t) = \alpha u_1(t)+\beta u_2(t)$ yields the output $x(t) = \alpha x_1(t) + \beta x_2(t)$ for any inputs $u_1(t)$ and $u_2(t)$ and any scalars $\alpha$ and $\beta$.

Therefore, the system is linear, as desired.


Resample audio signal with a low factors

You can resample to any, including floating point approximations of irrational rates, using reconstruction by Sinc kernel interpolation (realistically, by using a windowed Sinc interpolator).


Autoregressive modeling (linear prediction) of electrical transmission lines?

This is a neat idea and could probably be made to work - it seems like a variant on time-domain reflectometry, which exploits correlations of this kind https://en.wikipedia.org/wiki/Time-domain_reflectometry. However, there's an important difference between your two use cases: we can normally inject a test signal onto an electrical transmission line, while we normally can't do that in the vocal tract. Measuring the cross-correlation of the injected and reflected signals is generally going to be more reliable than measuring the cross-correlation of noise.


Butterworth filter transfer function in time domain

A plot of the normalized impulse responses, for the n = 2 through 10 Butterworth low pass filters, are given by H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976, p. 113. This is shown below. They do not give the h(t) expressions in the book, at least where I have looked thus far. They say the figure is from Zverev's earlier book: A.I. Zverev, "Handbook of Filter Synthesis", Wiley, NY, ©1967.

EDIT: Looking back in my old computation notebooks, I found the impulse expression for the 3rd order Butterworth low pass filter. This pic is from 1983:

Hopefully, this is algebraically the same as your result.

EDIT 2: Gregg 1 gives $h(t)$ for the $M = 2$ Butterworth LPF:
$$h(t) = e^{-x}sin(x)U(t)$$
where $x = t/τ√2$, $τ = RC = 1/(2πB)$ and $B$ is the corner frequency in Hz. Butterworth LPFs are fine for many purposes, but not particularly good for others, e.g., they have relatively poor 1% settling times with step inputs 2.

References:

W.D. Gregg, “Analog and Digital Communication”, Wiley, NY, ©1977.
E. Voigtman, J.D. Winefordner, “Low-pass filters for signal averaging”, Rev. Sci. Instrum. 57 (1986) 957-966.

EDIT 3: Even though Matt has provided an excellent answer, for completeness I decided to test my 1983 $h(t)$ expression (see above) by comparing it with the numerically generated impulse response from a computer simulation $H(s)$ transfer function block. As per the OP, I assumed $2πB = 1$. The next figure shows that the results are the same.


DFT without knowing the exact value of sampling rate?

If you don’t know the sample spacing, you don’t know the sample rate.

If you know that samples are uniformly spaced, you can do a DFT and get relative frequencies. You won’t be able to map one dft bin to a specific frequency (in Hertz) unless there is some additional side info (perhaps a known reference hidden within the data).


How can you intutively figure what the impulse response will behave like?

Here's an idea...

In the first case, the impulse response has only positive coefficients. Groups of nearby samples will be averaged to yield the output signal. Thus, short term bumps/variations will be "averaged out", so the signal will be smoother. Furthermore, it looks like all coefficients sum up to 1 (there about 30 of them, in a triangle-like shape, the maximum is about 0.06, I suspect they sum to 1), so the overall amplitude / magnitude of the output waveform will be preserved. Last bit of information: the impulse response has its center of gravity at $n=15$, so the output signal will be delayed by 15 samples.

The second impulse response looks like a kronecker delta (at $n=15$) minus the impulse response from the previous example. Because convolution is a linear operation, the result will be the original signal shifted by 15 samples, minus the smoothed signal obtained in the previous example. In other words, the result will be whatever has been smoothed out in the previous example because of averaging - only the fast wiggles will remain.

Summing the output signals from example 1 and example 2 will yield the input signal delayed by 15 samples.


Artefacts from playing ultrasound through PA system

First of all, like Cal-linux mention in their comment, I would strongly suggest intermodulation distortion for the frequencies that are not multiple integers of the fundamental.
Now, regarding the sweep tone artefacts. Since you are talking about high frequencies we are talking about a high-frequency driver. Nowadays, most high-frequency drivers are domes with some kind of compression driver attached to them. It has been proposed by Hubbard that compression drivers generate subharmonics. By visual inspection (as you also state) part of the energy of the sweep seems to be present at the first subharmonic $\frac{1}{2} f$, which could potentially be somehow related to the results Hubbard presented.
Further investigation is needed to be able to say for sure but this is indeed a very interesting observation. I work as a P.A. engineer, so I will try to duplicate the effect and work more on it. I will post back if I manage to reach any conclusions.


In theory, M[0] doesn't exist. Your discrete derivative is only defined on [1,127].

As a practical solution, you can define M[0] as a null element of your signal. For example if you expect values between -1 and 1, take M[0] = 0.

In the case of real-time, there is two options:

Either you save the last value of your buffer so you can take M[0] = last for the next ones, and initialize last as 0, or whatever seems coherent with your expected signal;
Either you start computing the derivative at M[1] and your processing will have a latency of 1 sample.


How to label sound events in noisy .wav files

Automatic detection of birdsong is a well studied problem in machine-learning applied to ecoacoustics/bioacoustics. State of the art methods use Convolutional Neural Networks on spectrograms. In binary classification of bird/no-bird they reach AUC ROC of around 0.90 in noisy and mismatched conditions. The DCASE2018 Bird Audio Detection challenge is a good starting point. Open-source code is available for many of the systems, including the challenge winner. Several of the systems could be adapted to Audio Event Recognition to give time-information of predictions.

Note that to verify your approach is working well, you should always also have a manually labeled test-set.


Why is the second-order-section cascade approach giving "worse" results?

I don't think that has anything to do with SOS vs transfer function representation, but it appears to be an artifact how filtfilt()  handles initial conditions. The documentation is annoyingly vague about this and simply says it

  filtfilt minimizes start-up and ending transients by matching initial
  conditions.

If you implement filtfilt() manually you get the expected result: both results are basically correct, but SOS has less  numerical noise. See code snippet below.

In general, using filtfilt() is tricky: ideally you have an input signal that has lots of zeros at the beginning and the end, where "lots" means more than the length of the impulse response. Otherwise the results are highly dependent on the initial states of the filters and how you manage them.

Numerically, SOS is always better. Take a look at your "a" coefficients: the ratio of the largest to smallest coefficient is in the order of 1e10. Adding numbers with that different in magnitude creates numerically issues, even if you use double precision. For the SOS representation, that ratios is only about 14, that's about 9 orders of magnitude better!

%% signal
t = (0:.08:100)';
x = sin(2*pi*t);

%% original
n=15;
Wn=[0.5 3]/6.25;
[z,p,k]=butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
d1=filtfilt(sos,g,x);
[b,a]=butter(n,Wn);
d2=filtfilt(b,a,x);

%% manual implementation of filtfilt()
y = filter(b,a,x);
y1 = filter(b,a,flipud(y));
y = g*sosfilt(sos,x);
y2 = g*sosfilt(sos,flipud(y));
plot([y2 y2-y1]); legend('result','difference');


How does summation over time samples influences DFT?

If the length of your data is $NM$ long (regardless of $N$), your DFT should be $NM$ long, so the correct equation is (1).

The explanation is simple: how long is your data? Use that length for the DFT.

You may want to define:
$$
y'_n = \left \{ \begin{array}{cl}
y_n, &  0 \le n \le N-1\\
0, &  n \ge N\\
\end{array}
\right .
$$
and operate on that because $y_n$ for $n\ge N$ is not defined.


This really comes down to how you define LTI.
Far as I can tell, a time-discrete system being LTI means its in- and output are related by a countably large system of ordinary differential equations.
Since each of these transforms into a polynomial through the Laplace-, Fourier- or z-transform, the transfer function, in which ever domain it exists, is always a quotient of polynomials. Coefficients of polynomials are countable.
Since polynomials have, including their multiplicity, exactly the number of roots as their degree, and poles are just the roots of the denominator polynomial, we'll always have countably many roots – and thus, poles.

Not quite sure you'd want to relax the constraint that poles be countable, anyway: If you do that, transfer functions cease to be meromorphic, and that means the residue theorem doesn't apply, and as far as I can intuit, that leads to unsolvable integrals that you'd need to define the $z$-transform to begin with...


Correct color values on image taken with no infrared filter

Unfortunately, no practically known signal processing method could provide what you wanted yet. Instead, you must use a hardware optical IR filter that would block the incoming IR waves, before they reach into your sensor.
An optical image is an array of intensity values captured by the sensor pixels which are illuminated by the incoming electromagnetic wave during the exposure time. During the capturing process, the static image sampling is performed on the 2D sensor surface, unlike an audio sampling which captures the time based variations of the incoming acoustic pressure wave.
Now the static-image sampling does not capture the time variations of the incoming electromagnetic wave. So those array of pixels represent only the spatial distribution of the illumination intensity $I(x,y)$ accross the sensor surface, over which an image was formed by a focusing lens system. Therefore, time-based property of light's wavelength is not represented in the samples, and cannot be processed later for modification.
Note that the camera system also samples the wave along temporal (time) dimension, the motion-picture; however it's typically about a few hundred frames per second at max. Given the required bandwidth of the visible spectrum; about $10^{14}$ Hz (visible light in [400-700] nm, plus IR spectrum to at least 1000 nm for near IR) this would require a sampling rate on the order of $10^{14}$ Hz. Such high data-rate is not practical for commercial applications of typical imaging devices.
So, for your Sony NoIR cam, the temporal-spectral information of the incoming electromagnetic wave is not captured, not represented by the image pixels, and therefore cannot be processes by any further means.
You should therefore physically place an IR blocking filter in front of the camera sensor, or the lens. Old GSM phones with built-in cameras, may incorporate an IR filter..?


Can Principal Component Analysis (PCA) Solve the Cocktail Party Problem?

The Cocktail Party Problem is a Blind Source Separation (BSS) problem.
Given a linear mixture of signals:
$$ \boldsymbol{y} \left[ n \right] = A \boldsymbol{x} \left[ n \right] $$
We're trying to estimate the signal $ \boldsymbol{x} \left[ n \right] $.
The model can get even more complex with $ A $ being time varying:
$$ \boldsymbol{y} \left[ n \right] = A \left[ n \right] \boldsymbol{x} \left[ n \right] $$
We have 3 main approaches to this problem:

Probabilistic Approach
Looking at the signals as an ensemble of points of a distribution and find the linear coordinate transform to guarantee some property. The PCA approach tries to remove correlation (2nd moment information) while the ICA tries to remove correlation in higher moments (Basically statistical independence).
Time Signal Processing Approach
In case of 2 signals one of them being a reference we can use the adaptive decorrelation filter. Basically we're after removing any time correlation from the signals.
Spatial Signal Processing Approach
We can utilize the known location of the microphones in the room to create adaptive beamforming. The idea is that delayed adaptive summation of the data can change the spatial curve of the array and making matching a certain direction.

Of course in late years we can find work on the subject utilizing Deep Learning approaches. Their main advantage is being able to incorporate additional information (Like using the properties of the signal, be it in a certain language, incorporating visual data on the scene like images and videos [Who is moving the lips when?]).
This is a vast subject and the main idea is to tailor the solution to the specific case of yours.
Modern Robust ICA and IVA (Independent Vector Analysis) can be very effective.
I'd try them first unless you have the case which matches the Adaptive Filter (Which can be proven to be matching the Beam Forming solution under some conditions).


There are many approaches to achieve this. For example, check the following work by Pederson et. al.:

Evaluation of contrast measures in relation to observers perceived contrast Marius Pedersen1, Alessandro Rizzi2, Jon Yngve Hardeberg1 and Gabriele Simone1.


Synthetic Aperture Radar (SAR). How to do Range Compression?

What format is your input data in? It sounds like you have a phase history signal as it's often called. Assuming that the transmitted waveform was a linear FM chirp, has it already been deramped? Some receivers perform deramp-on-receive processing, which removes the chirped nature of the waveform. In that case, range compression is as simple as performing a 1-D FFT on each pulse return (which appears to be one per column in your example).

If you instead have a raw collection from the receiver that hasn't been deramped, then you expect to have an input signal that consists of many reflections of the transmitted waveform, separated in time according to the round-trip delay to the various reflectors in your scene. In this case, in order to perform range compression, you need to generate and apply a matched filter to each pulse return (again, a 1-D convolution of each column of your dataset with the matched filter).


just implement the Modified Bessel function.  it's easy.

i always like my window definitions centered about zero, since pretty much all of them are even symmetry.

i'll do this in discrete-time, but it's essentially the same thing in continuous-time.

Kaiser window:

$$ w[n] \triangleq \begin{cases} \frac{1}{I_0(\beta)} I_0\left(\beta \sqrt{1 - \left(\frac{n}{M/2}\right)^2} \right) \quad \quad & |n| \le M/2 \\
0 & |n|>M/2 \\
\end{cases} $$

$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $M+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $M$ and, in my centered and symmetrical case, must be even). $\beta$ is a "shape parameter" and O&S recommend this heuristic:

$$ \beta = \begin{cases} 0.1102 \cdot (A-8.7) & A>50 \\
0.5842 \cdot (A-21)^{2/5} + 0.07886 \cdot (A-21) \quad & 21 \le A \le 50 \\
0.0 & A<21 \\
\end{cases}$$

$$ M = 2 \left\lceil \frac{A-8}{4.57 \cdot \Delta\omega} \right\rceil $$

$A$ is the desired stopband attention in dB and $\Delta\omega$ is the desired width of the transition band in normalized angular frequency.

finally, the Bessel is evaluated as:

$$ I_0(x) = 1 \ + \ \lim_{K \to \infty} \ \sum\limits_{k=1}^{K} \left(\frac{x^2}{4}\right)^{k} (k!)^{-2} $$

when you evaluate this with a computer, pick a $K$ decently large (my guess is that $K=32$ is good enough) and evaluate the summation starting with $k=K$ and work it backwards to $k=1$ to keep numerical accuracy.  you might want to use Horner's method.

$$ I_0(x) \approx  1 + x^2\left( \tfrac{1}{(1!)^2 \, 4^1} + x^2\left(\tfrac{1}{(2!)^2 \, 4^2} + x^2\left(... + \, x^2\left(\tfrac{1}{((K-1)!)^2 \, 4^{K-1}} +  x^2 \tfrac{1}{(K!)^2 \, 4^K} \right) \right) \right) \right) $$

you can evaluate all of the $(k!)^{-2}$ in advance with a short table.

Alright, someone made me do some work.  This took about 45 minutes to code up and debug.  So here is my MATLAB code for implementing the 0th-order Bessel function of the first find (which is $I_0(x)$ above):

function y = mybessel(x)
%
%   Computes the 0th-order Modified Bessel function of the first kind
%

    K = 32;

    bessel_coef =  zeros(1,K);

    kfac = 1;
    two_to_the_k = 1;

    for k = 1:K
        kfac = kfac * k;
        two_to_the_k = two_to_the_k * 2;
        bessel_coef(k) = 1/(kfac*two_to_the_k)^2;   % compute power series coefficients in advance
    end

    x = x.^2;

        y = x .* bessel_coef(K);
    for k = K-1:-1:1
        y = x .* (bessel_coef(k) + y);              % Horner's method
    end

    y = 1 + y;

end

and here is the test code:

x = linspace(-16, 16, 32*4096+1);

I_0 = besseli(0, x);            % MATLAB's modified bessel

y = mybessel(x);                % my bessel

figure(1)
plot(x, I_0)
hold on
plot(x, y)
hold off

figure(2)
plot(x, y - I_0)                % plot error

with results:

and error:

for $|x| \le 16$, the relative error is less than $10^{-15}$.  the error increases with increasing $|x|$.


Convolution and Cross Correlation Blackbox


How do you model a rolling shutter in the frequency domain

A rolling shutter introduces a displacement of the lines of an image that is proportional to the relative velocity of the camera (or target) and the velocity of the shutter. In terms of "simulating" this through a video, it translates to "borrowing" a few lines of pixels from each frame to compose one final frame. Notice here how the shearing produced by the rolling shutter, which could be seen as a displacement of $x$ in 2 directions just got translated to a simple "vertical" motion in 3 dimensions.

So, all that we have to do is look at what exactly this operation does in the time domain and see how (if possible) we can replicate it in the frequency domain.

From your question, I suppose that you want to create one image (let's call it $I$) from the "integration" of the frames in $a$. Before we construct $I$, we have to make some assumptions about how $a$ was acquired. So, let's assume that $a$ was acquired with a "perfect" camera. One that does not suffer from rolling shutter itself and one that shoots at a very high frame rate. So the period of sampling $T_s = t_k - t_{k-1}$ of $a$ is $T_s \ll T_e$, where $T_e$ is the frame exposure time (the time it takes for the shutter to roll and create one image). Furthermore, let's assume that the "width" of the shutter, in terms of image lines, is $W_s$.

Then, just considering the shutter, this is indeed a "white" rectangle that moves parallel to axis $y$ across all time instances of $a$ and it looks like:

$$S(x,y,t) = rect\left(\frac{x-\frac{N_x}{2}}{\frac{N_x}{2}},\frac{y-t \times T_e}{\frac{W_s}{2}}\right)$$

This is a two dimensional rectangle. It is centred in the middle of the $X$ dimension of the frame and is as wide as the frame and in the $Y$ direction it is $W_s$ lines tall, moving downwards.

The time instances that make up $I$ are:

$$I(x,y,t) = a(x,y,t) \times S(x,y,t)$$

But if you "replay" this, all that you see is a band of pixels from each frame moving downwards.

To get rid of $t$ in $I$, you need to integrate this signal with respect to $t$. In other words, you have to sum all pixel time courses to generate just one image that is independent of $t$.

Now, this is a linear time invariant system of operations in three dimensions. To find its coefficients, you need to first transform $S_{FT} = \mathscr{F}(S(x,y,t))$.

To produce $I$ we multiplied $a$ with $S$. Multiplication in the time domain is convolution in the frequency domain. Therefore to express $I(x,y,t)$ (notice the $t$ still there):

$$I_{FT}(k,m,f) = \mathscr{F}(a(x,y,t)) \ast \mathscr{F}(S(x,y,t))$$

Here of course, notice the change of variables, $m,n,f$ represent their frequency domain $x,y,t$ counterparts.

Right, so now, we want to integrate the three dimensional complex $I_{FT}(k,m,f)$ with respect to time. We look up the properties of the Fourier Transform and realise that integration in time is equivalent to multiplication in frequency by a $\frac{1}{j \omega}$ factor plus a static DC component scaled by $\pi \delta(\omega)$. So, all we have to do now is multiply our $I_{FT}$ coefficients and offset them by these factors, in the direction of the time vector. Or in other words, in the dimension that corresponds to frames.

$$I_{FT}(k,m,f) = I_{FT}(k,m,f) \times \frac{1}{e^{2 \pi i f/N_f}} + \pi I_{FT}(0,0,0)$$

(And of course then, above, obtain the inverse to get back into "time" domain)

Of course, you can express the two as one operation and do the transformation in one step.

Hope this helps.


What is the purpose of the following code?

I think fft_max is a kind of equalization table so that the visualization would be "flat" over all frequencies for a typical input. It does not look carefully crafted as it is so spiky. Probably the author just used some example audio input and kept book of the largest encountered values of the magnitudes of the bins, hence the name of the table. If you need a different length fft_max, it would be fine to interpolate the original fft_max say by linear interpolation, or to roll your own (not a cigarette).

You could use the synthetic pink power spectral density $a^2/f$, where $a$ is a normalization constant and $f$ is frequency. In terms of magnitude which you need to put into code, that would be $a/\sqrt{f}$, the square root of power. Put some reasonable value to $f = 0$, for example copy the value at the first non-zero frequency.

Figure 1. Pink magnitude spectrum (red diamonds) compared to the original fft_max (blue dots).


How to choose a fixed adaptation step for decision feedback equalizer

I'm leaving the answer here, if somebody ever stumbles upon the same question.
About LMS: Both DFE and ARC employ Least Mean Squares (LMS) adaptive algorithms: DFE is an adaptive filter and ARC can be considered as an additional tap of that filter. Hence the adaptation step size should be selected accordingly.
The general update equation for the LMS algorithm is the following:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot e(n) \cdot \mathbf{x}(n),
$$
where $\mathbf{w}(n)$, $e(n)$ and $\mathbf{x}(n)$ are respectively filter weights, an error and tap values at $n^{th}$ sample and $\mu$ can be a constant scalar.
With $p$ being the filter order, the bounds for $\mu$ can be defined as
$$
0 < \mu < \frac{2}{(p+1)E\{|x(n)|^2\}},
$$
where $E\{|x(n)|^2\}$ represents the power in $x(n)$ and can be estimated using an average:
$$
\hat{E}\{|x(n)|^2\} = \frac{1}{N}\sum_{k=0}^{N-1} |x(n-k)|^2.
$$
Then $\mu$ value can be chosen as a fraction of $\mu_{max}$, e.g. $0.125 \cdot \mu_{max}$.
However, in the sign-error LMS, which is an LMS algorithm with reduced complexity with the following update equation:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot sign(e(n)) \cdot \mathbf{x}(n),
$$
replacing an error value with its sign changes the magnitude of the correction. If $e(n)>1$, then the magnitude of the correction is smaller compared to the standard LMS. Hence the $\mu$ value can be chosen to be larger.
To the questions:

For example in a question, the input signal is a sequence of evenly (presumably) distributed PAM4 symbols, i.e. "-3", "-1", "1" and "3", which gives a power equal 5. This gives $\mu_{max} = \frac{2}{(p+1)\cdot 5}$.

DFE and ARC employ different LMS algorithms: sign-error LMS and standard LMS respectively. This calls for separate $\mu$ values for DFE and ARC. Hence the recommended adjustment would be to use sign-error LMS for both functions. Not only does this simplifies the $\mu$ value selection, it also simplifies the further implementation.

Considering the need to simplify the hardware implementation of the algorithms, $\mu$ can be chosen to be a power of two. For this particular example, $\mu$ equal to $\mu_{max}$ rounded to one significant binary digit happened to be a good choice.
References: A good reference on the LMS algorithms and step size selection is Chapter 9 "Adaptive Filtering" in Statistical Digital Signal Processing and Modeling by M.H. Hayes. Another good option is a Least_mean_squares_filter entry on Wikipedia.


What are the spike removal techniques in matlab without using filters

based on your results, you can try a median filter, in Matlab - medfilt1; which takes n-length samples of your series x at a time and computes the median, and apply the median value to each value of x. It is a very common non-linear filtering technique.

try something like

result = medfilt1(signal,6)

you and alter the length of evaluation from 6 to whichever number suits the data.


Both real and complex signals can be passband.

Assuming $x(t)$ and $y(t)$ are baseband signals with bandwidth much less than $f_c$, then $f(t)$ is a passband signal, even if it is complex.


linear convolution toeplitz matrix vs circular convolution toeplitz matrix

They are in general different. For two signals of lengths $N$ and $M$, linear and circular convolution are equivalent if the output is specified to be of length $N + M - 1$ with the appropriate padding. Convolution via the DFT is inherently circular, which is why padding must be done before the inverse DFT to yield the linear convolution. So, this is a special case where they are the same.
If your goal is to always yield linear convolution, then don't worry about forming a circular Toeplitz matrix since the result will be the same when using the regular Toeplitz and is simpler to do so.
Below is some sample code and output where we form regular and circular Toeplitz matrices with a specified output of length $N + M - 1$:
%% Toeplitz Convolution

x = [1 8 3 2 5];
h = [3 4 1];

% Form the row and column vectors for the Toeplitz matrix
r = [h zeros(1, length(x) - 1)];
c = [h(1) zeros(1, length(x) - 1)];

% Toeplitz matrix
hConv = toeplitz(c,r)

% Compare the two types of convolutions
y1 = x*hConv
y2 = conv(x, h)

hConv =

 3     4     1     0     0     0     0
 0     3     4     1     0     0     0
 0     0     3     4     1     0     0
 0     0     0     3     4     1     0
 0     0     0     0     3     4     1

y1 =

 3    28    42    26    26    22     5

y2 =

 3    28    42    26    26    22     5

%% Toeplitz Circular Convolution

% Convolution length
n = length(x) + length(h) - 1;
numElementDiff = n - length(h);

% Set up the circular Toeplitz matrix
c = [h(1) fliplr([h(2:end) zeros(1, numElementDiff)])];
hConvCirc = toeplitz(c, [h zeros(1, numElementDiff)])

% Compare the two types of convolutions
y1 = [x zeros(1, length(c) - length(x))]*hConvCirc
y2 = cconv(x, h, n)

hConvCirc =

     3     4     1     0     0     0     0
     0     3     4     1     0     0     0
     0     0     3     4     1     0     0
     0     0     0     3     4     1     0
     0     0     0     0     3     4     1
     1     0     0     0     0     3     4
     4     1     0     0     0     0     3

y1 =

     3    28    42    26    26    22     5

y2 =

    3.0000   28.0000   42.0000   26.0000   26.0000   22.0000    5.0000

Here we're testing three things:

Linear convolution conv() is equivalent to performing the matrix multiplication with the appropriate Toeplitz matrix.
Circular convolution cconv() is equivalent to performing the matrix multiplication with the appropriate circular Toeplitz matrix.
Output length is specified as  $N + M - 1$, so we see that linear and circular convolution are equivalent.

If you are going to perform circular convolution of varying sizes, then you must form the Toeplitz matrix differently. This usually involves some type of padding with the matrix entries themselves or the signal(s) being operated on. Mathworks has a good summary of Toeplitz matrices here and linear vs circular convolution here.


Please see help conv (especially third parameter) and help fft. Valid code sample:

x = zeros(1,256); x(120:170)=1;
y = zeros(1,256); y(120:170)=1;
z = conv(x,y,'SAME'); size(z)
z1= ifft(fft(x).*fft(y));
z2 = real(fftshift(z1));
plot(z2);


The two methods are equivalent over discrete images and in both of them it is "...tricky to select the correct boundary radii to really count the pixels that correspond to a certain frequency".

Whether "selecting pixels by distance" or converting from rectangular to polar coordinates, there will always be the problem of what do you do when you have to "land" between pixels. That is, when you have to "land" at an image site, within the acquired data range but between known values.

The answer here is of course "interpolation". Whether you are trying to select all harmonics that are $\approx 14.48765\ldots$ pixels from the centre of the image or you are converting $(x,y)$ pixel $(15,5)$ to its $(\rho, \theta)$ corresponding location of $(\approx 14.41479\ldots , \approx 81.36745\ldots)$, you still have the problem of "landing" between pixels.

You could of course $round(\cdot)$ everything, which would result in a nearest neighbour interpolation but that doesn't look very attractive. You are likely to have to use more complex interpolation schemes.

By the way, by convention, in digital image processing, the two dimensional discrete fourier transform's output is shifted by $\frac{N}{2}$ in both dimensions in order to bring the low frequencies to the centre of the "image". Therefore, you are likely to have to use fftshift too.

As far as the last part of your question is concerned, there should be no need for "scaling" if float or double data types are used, both to calculate the spectrum and its sum over different frequencies. If you are asking what sort of scaling you should do in order to be able to display the image (?), then, please note that you may have to use a logarithmic transformation on the spectrum values before you convert the image to something like unit8 to prepare it for display. Happy to clarify this later part more if more information emerges.

Hope this helps.


What is the essential bandwidth with the given function?

Find the Fourier Transform $F(\omega)$ and perform the following integration.$\frac{1}{2\pi}\int_{-W}^{W}|F(\omega)|^2\mathrm{d}\omega=0.7E$, where $E$ is the total Energy of the function given by: $E=\frac{1}{2\pi}\int_{-\infty}^{+\infty}|F(\omega)|^{2}\mathrm{d}\omega$.The value of $W$ is the 70% Bandwidth.


Generating a time series given a transfer function

I think that procedure should be as follows:

1.Generate the white (Gaussian) noise  signal - randn function in MATLAB.

2.Having a transfer function of your filter divide it by $z^4$ (highest power in denominator) to get a proper representation for IIR filter:

$$H_d(z)=\dfrac{0.1584\cdot z^{-1} −0.3765\cdot z^{-2} + 0.2716\cdot z^{-3} −0.0534\cdot z^{-4}}{1 −2.9951\cdot z^{-1} + 3.0893\cdot z^{-2} −1.1930\cdot z^{-3} + 0.0988\cdot z^{-4}}$$

3.Now you can extract the coefficients according to equations that can be found here: IIR.

$$ H(z) = \dfrac{\sum_{i=0}^{P}b_i z^{-1}}{1+\sum_{j=1}^{Q}a_j z^{-j}}$$

This will give you following vectors:

a=[1 -2.9951 3.0893 -1.1930 0.0988] % denominator
b=[0 0.1584 -0.3765 0.2716 -0.0534] % numerator

4.These are coefficients of your filter that you can apply to the white noise signal. If you are using MATLAB/Octave, then use the filter function. Something like:

y=filter(b,a,x)

where x is your white noise signal generated previously.

The amplitude response of your filter is shown below - because you are going to apply it to white noise (having all frequencies equally probable), then spectrum of your filtered signal should be corresponding to amplitude response of your filter. There is a high gain in lower frequencies which makes somewhat sense (wind gusts are mostly energetic in lower range):

Judging from your article, now you only have to normalise the filtered white noise and scale it accordingly.
Good luck!


Creating a convolution kernel from filter description

The key point in the description on the linked-to page is:

  This function will filter the image by the nonlinear contra-harmonic mean method.

that the method in question is nonlinear.

That means convolution kernels are not really applicable, because they are generally only useful for linear, time-invariant systems.

Having said that, you may be able to use a simple linear convolution kernel in the process of calculating this mean.

You need to calculate:
$$
\sum_{(i,j) \in M} A(x+i,y+j)^Q
$$

where $Q$ is either $P$ or $P+1$ and $M$ is an $N\times N$ square mask ($\square$).

One way to do this would be to:

Calculate the image to the power $P+1 = I_1 =   A(x,y)^{P+1}$.
Calculate the image to the power $P =  I_2 =  A(x,y)^{P}$.
Form
$$
I_\mbox{NUM} = \square * I_1\\
I_\mbox{DEN} = \square * I_2\\
\mbox{Contra-harmonic mean} = \frac{I_\mbox{NUM}}{I_\mbox{DEN}}
$$

Here, the convolution kernel, $\square$ is just:
$$
\square(x,y) = \left\{
\begin{array}{cl}
1 & |x|,|y| < N/2\\
0 & \mbox{otherwise}
\end{array}
\right .
$$
i.e. it's just a square $N \times N$ kernel of all ones.


Are there paper about black box testing of optic flow algorithms?

There is Middlebury, which evaluates almost all the algorithms and their evaluation results are online:

http://vision.middlebury.edu/flow/eval/results/results-n1.php

Or simply click any of the links at: http://vision.middlebury.edu/flow/eval/

They also have a paper on their research and evaluation methodology.


Block code error control with $d \ge c+u+1$

You have to realize that $d$ is specifically the minimum distance in your code, i.e. the distance between the to code words that are closest together.

And $u$ is not exactly the number of "uncorrectable errors", but the maximum number of errors that can definitely be detected, but not corrected correctly.

Note that "distance" here is the hamming distance, i.e. the number of bits you'd have to flip until you end up with another codeword.

  We conclude that all these combinations below are possible: ...

False. The first three lines ($c\in\{6,5,4\}$) are impossible.

Also, your table doesn't make much sense. Given a code, you'd always correct as many errors as possible¹, so for a given code, $c$ is given, and $u$ is given. The fact that

$$\begin{align*}
c &=& \#(\text{detectable bit flips, which ARE correctable})\\
  &=& \#(\text{total detectable bit flips}) &- \#(\text{detectable bit flips, which AREN'T correctable})\\
  &=& (d-1) &- u\tag{*}
\end{align*}$$

is just your definition of $c$.

A $k=7$ repetition code can never correct 6, 5 or even 4 bit errors; I find "you can't flip more than $\left\lfloor \frac d2 \right\rfloor$ bits until it looks more like another codeword than like the original" highly intuitive, so the only way I can explain this is by example:

take your 0 code word, and flip 4 bits:

0 0 0 0 0 0 0
0 1 1 1 0 1 0

Now, that looks much more like the 1 keyword, doesn't it? In fact, 0 1 1 1 0 1 0 has a hamming distance of 4 to the 0 code word, and a hamming distance of 3 to the 1 code word, so it's closer to 1, and will correctly be identified as error, but falsely get corrected to 1.

Looking at this, any error is detectable, as long less than the minimum distance bits are flipped (for any code, pretty logical, isn't it?). So,

$$c_\text{max} = d -1 \text.$$

Now, your formula simply "happens" when you put in your definition of $c$ and $u$, see $(*)$ above.

There's not much to explain here: for repitition codes $d=k$, and you can always correct $\left\lfloor \frac d2 \right\rfloor=\left\lfloor \frac k2 \right\rfloor$ bit errors and detect $d-1=k-1$ errors.

One can't explain your other combinations, as they are basically wrong for any maximum likelihood decoder. You can of course work "worse" than possible, but especially for the  repitition codes, summing up the received bits and deciding whether the result is larger than $\frac k2$ gives you $c=\left\lfloor \frac k2 \right\rfloor$ and thus, $u=(d-1)-c=(k-1)-\left\lfloor \frac k2 \right\rfloor =\left\lceil \frac{k-1}{2}\right\rceil $, and, in your $k=7$ case, $c=u=3$.

¹: point is that there's many codes where you can't, with reasonable effort, correct as much correctable errors as the code would allow you too, that means that $c \le \left\lfloor \frac d2 \right\rfloor$ in general, and that's why $d\ge c+u +1$ and not $d = c+u +1$. However, repetition codes are so trivial, we can truly Maximum-Likelihood-decode them and correct $\left\lfloor \frac d2 \right\rfloor$ errors.


confused about time shifting property of Laplace Transform

If you have written the function correctly then its Laplace transform could be found very similary to your first example:

Given
$$x(t) = e^{-2 t} u(t-1)$$
its Laplace transform could be found as follows. First denote the signal

$$x_0(t) = e^{-2} e^{-2t} u(t) $$

then its obvious that
$$x(t) =x_0(t-1) $$

Using the tables and properties to conclude:

$$X(s) = e^{-s} X_0(s) $$
$$X(s) = e^{-s} \frac{ e^{-2} }{s + 2} = \frac{ e^{-(s+2)} }{s + 2}   $$


Gabor can be used as a filter bank in order to extract texutre
features from an image (or image patch).
To compare images (or patches from image) using gabor I would follow the following pipline:

compute a filter bank from the Gabor (different scale and orientation)

post process to create new feature. you can compute the energy, and first statistics measures on the filtered images.
now that you have a vector of features for each image you need some measure to compare tham. You can use simple L2, or fisher criterior as some papers offer Comparison of Texture Features
Based on Gabor Filters

Please also take a look at the following question that is simmilar to yours.


Why is this FFT code not working to extract frequencies from a wav audio?

Your DFT is working as expected at extracting frequency information.  One musical note at one pitch usually contains a lot of spectral frequencies: from harmonics and overtones plus broadband noise from any transients, such as note attack and decay.  Furthermore, the max magnitude frequency from an actual sound recording is often not the same as the heard note pitch frequency (one possibility is the magnitude peak is instead a strong overtone).

Instead, search for pitch detection/estimation algorithms.  The good (working, more reliable) ones will do a lot more than just an FFT peak magnitude search.


What you need to know is that the signal $$s(t)=e^{j2\pi f t}$$ has frequency $f$. At any time $t_0$, $s(t_0$) is a point on the unit circle. As time advances, this point moves on the circle, and it completes $f$ turns in one second.

Knowing this, you can easily see that the maximum frequency in your signal $x_c(t)$ is 9000 Hz. That is because, for $k=9$, you can write the exponent as $j2\pi 9\cdot 1000 t$.

The sampling frequency you mention has no bearing on the signal's maximum frequency. However, if you intend to sample this signal, you need to make sure that the sampling frequency is at least twice the maximum frequency, or $f_s=18000$.


Phase wrap when computing impulse response from frequency response

No problem at all. Inverse Fourier Transform is totally blind to phase wrapping. You can add any multiple integers of two pi to any phase and you will still get the exact same impulse response


Apply wavelet transform to analyse EEG signal

I'm not sure whether you're talking about Discrete Wavelet Transform (DWT) or Continuous Wavelet Transform (CWT). Both can be used on discrete signals similarly to DFT and DTFT, I'm not sure if anyone calls it the Discrete Time Wavelet Transform instead. In any case, the dyadic wavelet transform is the non-redundant one (number of samples in $=$ number of coefficients out), and that's what we usually assume by talking about DWT.

It looks like you need to do a bit more reading on the wavelet transform in general, because of your first question. Wavelet transform doesn't work like the Fourier Transform; rather, it's akin to the Short Time Fourier Transform, which is a function of two variables: frequency and time. Similarly, the wavelet transform is a function of scale and time. This means that you don't pick a scale and stick to it. You go through many scales (in dyadic case they increment in powers of two), and calculate the transform for each one.

To address your second question, edge effects are usually not that easy to deal with and there's a plethora of papers on the topic. If you simply want to know what portion of the transformed signal will by affected by them, and it looks like that's what you're asking for, this paper has a good discussion. One thing to keep in mind is that the DWT, like the Discrete Fourier Transform, has circular wrapping symmetry, so dealing with edges is undesirable if you're hoping to get a perfect inverse transform as well. You can look at this paper for a discussion of that, and also a way to eliminate edge artifacts.


Generation of a stream of data bits 0s and 1s at a certain bit or data rate

Multimedia signals such as video,picture, music, speech have bandwidths associated with their baseband analog source nature.

For digital processing this bandlimited analog source is passed through a analog to digital converter (into a sequence of samples), which operates at a clock-frequency determined by the signal sampling rate, $F_s$, which is also limited by the source bandwidth due to Nyquist-Shannon sampling theorem.

Furthermore for each sample, a number of bits are used to quantize (numerically resolve) it. The number of bits is dictated by the signal amplitude accuracy requirement. For digital video, it's mostly 8-bits per RGB channel. Or for audio it can be 16-bits per stereo L/R channel.

Hence, the multiplication of samlping rate $F_s$ and bit-per-sample gives you the raw source bit rate $R_b$. You may prefer to apply some statistical or perceptual data compression to reduce $R_b$ to a minimum.

Then you choose a digital modulation scheme such as FSK,PSK,GMSK,OFDM... and transmit those bits through your favorite medium. Before doing this, you may increse $R_b$ (by adding redundancies) for channel transmission error correction purposes. The result is your bit-stream (not necessarily logic signal levels of 0s and 1s though) that's traveling accross the universe...

I hope I clarified it.


Sub nyquist sampling, required number of samples for time sparse grouped signals

Compressed sensing does not assume any distribution of non-zero elements in the input vector (signal), so it makes no difference if your non-zero elements are near each other or located uniformly on the sampling grid. The number of measurement, as theory states for all sparse signals, depends on how many non-zero elements you have (sparsity degree), what is your noise in input and what recovery quality your need (exact or approximation) and etc.
therefore,

Does it make sense to perform compressed sampling if the non zero samples are grouped in time?

Yes!
I suggest the following video lecture of Baranuik:
Video Lecture On CS


Retrieve pixel values from the output jpeg image

Ask yourself: why do you think the TIFF and BMP files are so large while JPEG files are small?

Answer: JPEG is a lossy compression format and achieves small file size by removing details from your image imperceptible to the human eye. By design, JPEG cannot guarantee that the decompressed image will be pixel-by-pixel identical to the compressed image - it just preserves the overall content the human eye can perceive. You could try increasing the quality parameter, but this probably won't solve the problem in your case. Brightness levels will still be quantized, and color information downsampled.

Some "solutions":

Do not use JPEG, use a lossless compression format.
Do your image editing/modifications directly in the compressed domain. I doubt there are tools in Matlab to do that, and compressed-domain image modification is not as straightforward as writing into an array of pixels!


Spectrum analysis : how to represent FFT results

Actually, what do the squared magnitude FFT result represent? They represent the total energy in a given frequency bin. So to downscale your results to less bands, you can sum them. It's not even necessary to average the results - however it would have the same effect, as you map the results to an image, thus rescaling them.

However, note that if you need to rescale by something else than an power of two, you'll need to do it using supersampling, as the results won't map directly to your new bins.


How to Find the Probability of Kalman Filter States?

The idea of Kalman is it fuses the data with the optimal weighing to get results.
Optimal requires definition and in the case of Kalman it means that if your model is right this is the optimal.

Since we operate Kalman filter under the assumption the model is right no Kalman Filter will allow you doing what you want.

Now, it seems that in your case what you don't like is the transient phase when there are low number of measurements.
What I can suggest is that for this start phase you will built your own filter based on your weighing and after few iterations send your result to the Kalman Filter as a starting point.


For BIBO stability in the case of discrete time, there is a necessary and sufficient condition given by $\sum |h[n]| < \infty$ that is if the impulse response is absolute summable then the system is BIBO stable. Clearly
$$h[n] = \delta[n] - \delta[n-1]$$ and it has a finite support, the impulse response is absolutely summable and therefore BIBO stable.


Spatial Audio (Ambisonic) and calculating sound source position

I know nothing about Ambisonic decoding, but the Wikipedia page suggests:
$$
p(\theta_n) = w(t) + x(t) \cos(\theta_n) + y(t) \sin(\theta_n)
$$
Perhaps you could form that and vary $\theta_n$ to find the direction of highest power?

Note: I realize that's not what the page is really saying. It's saying that this is the signal to feed to a speaker at $\theta_n$... but that seems similar to finding the directionality of the source.  Also note that it does not find the distance of the source, just the $x-y$ plane angle.


Determining the noise floor of a signal in the frequency domain

You need to normalize your data based on the type of window you are using to obtain your frequency domain representation of the data.  The normalization differs depending on whether you are measuring a narrow band (your signal peak) or broadband (noise) signal.  Once you have the data properly normalized, the power of the narrow band signal can be read directly from the data.  The noise measurement must be estimated from the "noise floor" of the normalized frequency data.  Your noise power estimate will be 6dB less than the noise floor.  For a detailed discussion,

Go to this link: http://www.fhnw.ch/technik/ime/publikationen

Download the paper ""How to use the FFT for signal and noise simulations and measurements".


Linear Difference Equation and Method of Least Squares

The notation $ \hat{y} \left( k \mid k - 1 \right) $ usually means this is an estimated value of $ y \left( k \right) $ given all the available data up to time index $ k - 1 $.
So generally speaking, this is a prediction of one step in time of the data.

The case above also suggests linear estimation.
Namely, $ \hat{y} \left( k \right) $ is built using linear combination of all the given data.


Consider three points $A$ in image 1,$B$ in image 2 and $C$ in image 3.$H1$ is the homography matrix from $A$ to $B$ and $H2$ is the homography matrix from $B$ to $C$.
Then,

$$B = H1A$$

$$C = H2B$$

Substituting for B from first equation

$$C = H2(H1A)$$

Due to matrix associative property,

$$C = (H2H1)A$$

The homography matrix from $A$ to $C$ is given by $H2H1$.


I believe that your problem is using nextpow2 to set the fft length.  fft(x) returns a vector length(x). fft(x) is not constrained to a power of 2. If you want to set a particular size by zero padding, use fft(x,N). NFFT=length(Ax) is the actual size of DFT.

Using psd has advantages but if you use the proper size, your code will work (assuming I am correct) . Taking a single long DFT isn't an accepted way to estimate a psd, unless you are looking for very fine, high snr features.


Homography for (small baseline) Stereo Image Registration

For most cases (and I tend to say) yes, but this also depends on the geometry. If there is a perpendicular plane and two cameras observe different sides, the assumption will cease to hold.
Also, in the cases where objects are far and no crazy things are going on, you could also parameterize the intrinsics using an affine camera model.


At this moment i'm not hos i can't check my answer...
You say

  In the output array I get amplitudes of sin and cos functions.

Well that is you frequency array. Probably not the way want it. I suspect youre looking for the magnitude and maybe the phase.

the magnitude is quite easy.. it is the $abs(array (real + imag)$.
the phase $ \phi = \arctan \frac{imag}{real} $

in this way it is common to show the half the frequency array because of it's symmetry.
if you plot $array[0]$ to $array[N/2]$ you shoult get your frequecy array from $ \frac{1}{N}$ to $\frac{F_s}{2}$


Acoustic room simulation - Craven's hypothesis

Question 1: The equations (6.7), (6.8), and (6.9) are all consistent, so I suggest that the wording for the definition of $g$ is incorrect in the book.

Question 2: The brain just hears the time difference and, by experience, can figure out what that means in terms of distance. So the brain has effectively "calibrated" for whatever the usual speed of sound in air is.  I find that hearing things under water (where the speed of sound is about five times faster) makes things seem much closer.


Effective field in MRI when solving for Bloch equations

In the rotating reference frame, $B_z$ from the static magnetic field $B_0$ is transformed away, hence $B_z = 0$ if no gradient field is applied. All gradients change the $B_z$ value, as you described it with the slice selection gradient. This gradient causes $B_z$ to increase when you move along $z$. Note that it is irrelevant if you do this from the center of the magnet or somewhere close to the bore - you will have the same $z$-dependence at each position $(x,y)$.

For the readout-gradient (assuming it is along left-right direction that we'll call $x$), the same is true: $B_z = G_Rx$. If you go from left to right, $B_z$ changes. In this case, it is independent of the position along the bore, so moving along $z$ does not have an influence on the value, if you stay at the same $x$ coordinate.

A gradient field always acts along the $z$-component of the magnetic field. This is always a bit tricky, since gradients can be "applied along each direction".


The Fourier transform of sinusoids' products with possible other components

HINT:
It is easy to see that things can be simplified using the trigonometric product-to-sum identity in Equation $(1)$ below:
$$
\cos(\alpha)\cos(\beta) = \frac 12\big[\cos(\alpha+\beta) + \cos(\alpha-\beta)\big]\tag{1}
$$

In the first example
$$
\frac{1}{t^2}\cdot\cos(2πt)\cdot\cos(2πt) = \frac 1{2t^2}\big(\cos(4\pi t) + 1\big)
$$
From there you would want to visit the cosine modulation frequency-shift and differentiation properties of the Fourier transform, here in $(2)$ and $(3)$ respectively.
\begin{align}
\mathcal F\big\{x(t)\cos(2\pi f_0 t)\big\} &= \frac 12\big[X(f - f_0) + X(f + f_0)\big]\tag{2}\\
\mathcal F\left\{\frac{d^n x(t)}{dt^n}\right\}&=  \left(j2\pi f\right)^nX(f)\tag{3}
\end{align}

In the second example, using Equation $(1)$ you then have
\begin{align}
8\cos(20\pi t)\cos(40\pi t)\cos(80\pi t) & = 8\bigg(\frac 12\big(\cos(60\pi t) + \cos(20\pi t)\big)\cos(80\pi t)\bigg)\\
& = 4\big(\cos(60\pi t) + \cos(20\pi t)\big)\cos(80\pi t)\\
& = 2\big(\cos(140\pi t) + \cos(20\pi t)\big)\\
&\quad + 2\big(\cos(100\pi t) + \cos(60\pi t)\big)\\
\end{align}
With this you simply have sums and you don't have to think of convolutions, you have individual sinusoids at frequencies $10\ \rm Hz$, $30\ \rm Hz$, $50\ \rm Hz$, and $70\ \rm Hz$.


I would equalize both channels.  The most robust method is to cut both back to a lower bandwidth.  Then you only need to equalize the gains.
Assuming that your scanning by hand a 1-2 second sweep rate seems reasonable.  Unless your imaging, you can narrow the sweep down to smaller distances when a signal needs investigation.
A 10-50Hz bandwidth should accommodate human response times.  In this case you can use simple "RC" filtering to a fixed standard on both samples.
You mention construction variations; so first, you need to buffer yourself away from the coils.  Then apply a fixed RC filter.  Since you have isolated yourself both in time and impedance from the coil variations you have complete control of the balance between the two sides.
Of course, this might all be done in software but make sure your A/D has sufficient resolution (time and amplitude) and compliance.


What will the impulse response of a matched filter look like if the input is complex?

The impulse response doesn't change because it is the output of the filter when the input is an impulse.

When the complex input is applied to a complex filter the output is still determined by the convolution of the input and the impulse response of the filter. The addition/integration and multiplication operations that take place inside a convolution operation have to account for complex numbers.

I think it is difficult to gain intuition for the result of a matched filter in the time domain. Instead look at it from the frequency domain.
\begin{align}
 Y(\omega) &=&  H(\omega)S(\omega) \\
 & = &  e^{-j\omega t_0}S^*(\omega) S(\omega)\\
 & = &  e^{-j\omega t_0} ||S(\omega)||^2
\end{align}

In the case where the envelope of $S(\omega)$ is rectangular (as in your example) then $||S(\omega)||^2$ is equal to a rectangular pulse. The output then be a time delayed version of a sinc() pulse. The wider the frequency range of $S(\omega)$ then the narrower the main pulse of the sinc() pulse.

If you are trying to gain insight into how a matched filter is derived - one insight I can provide is that you have to understand that the Signal to Noise ratio (SNR) is defined differently at the input of the filter, than the SNR that we are trying to maximize at the output. At the output we are interested in the maximum value of the output signal, while at the input the SNR is measured using the signal's energy.


Does the hold part in sample and hold turns the signal from discrete to continuous?

The sample-and-hold circuit is a physical, analog electrical circuit and, by definition, its output is continuous.

The actual conversion to discrete time happens at the quantizer, not at the sampler. The output of the quantizer is a discrete sequence of numbers associated with time instants, so it can be said to be a discrete signal.


Is it possible to transmit and/or receive multiple types of signals at the same time?

Theoretically it's possible to do a frequency multiplex of two different protocols with one SDR transmitter. However, the practical main limitation is the speed of the digital-to-analog converter (DAC). E.g. for a WLAN and GSM band the channels that are closest on the frequency axis are at 1800 MHz and 2412 MHz, respectively. For generating such a signal you'd need a sampling rate of over 1200 MSa/s. Even if you find an SDR device that is equipped with an DAC that fast, it will be very challenging to filter the aliasing spectra generated by the DAC in order to avoid out-of-band radiation. So your sampling rate has to be even higher in practice.

If you want to transmit two independent data streams at bands very close to another, OFDM can be advantageous because of it's narrow-band characteristic. If the modulation is synchronized, the bands can even overlap. The latter is an case without practical relevance, though.


How does time shift correspond to phase change in a discrete signal?

If you have a signal

$$f[n]=\cos(\Omega_0n)$$

and you apply a time shift of $n_0$ you get

$$f[n+n_0]=\cos(\Omega_0(n+n_0))=\cos(\Omega_0n+\Omega_0n_0)=\cos(\Omega_0n+\phi)$$

where $\phi=\Omega_0n_0$ is the phase shift.

The other way around, if you have a phase shift of $\phi$, this is not always equivalent to a time shift of the original signal:

$$g[n]=\cos(\Omega_0n+\phi)=\cos(\Omega_0(n+\phi/\Omega_0))$$

which only corresponds to an integer time shift if $\phi/\Omega_0$ is integer, i.e. if $\phi/\Omega_0=n_0$ which results in

$$g[n]=f[n+n_0]$$

EDIT:
Re-reading your question, I think that the misunderstanding lies in the fact that you believe that the phase of a discrete-time signal must be integer. This is not the case. Imagine two continuous-time signals
$$x(t)=\cos(\omega_0t)\quad\textrm{and}\quad y(t)=\cos(\omega_0t+\phi),\quad\phi\in\mathbb{R}$$

Note that the following holds for any value of $\phi$:

$$y(t)=x(t+\phi/\omega_0)$$

So we can always express $y(t)$ as a shifted version of $x(t)$. Now imagine that we construct two discrete-time signals by sampling $x(t)$ and $y(t)$ at times $t_n=nT$ with some real-valued $T>0$:

$$f[n]=x(nT)=\cos(\omega_0nT)=\cos(\Omega_0n)\\
g[n]=y(nT)=\cos(\omega_0nT+\phi)=\cos(\Omega_0n+\phi)
$$
with $\Omega_0=\omega_0T$. Note that the phase $\phi$ is of course the same as before (and for this reason not necessarily integer). The difference between the continuous-time and the discrete-time case is now that $g[n]$ cannot in general be obtained from $f[n]$ by time-shifting because we can only shift by integers. The condition under which $g[n]$ can be obtained by shifting $f[n]$ is if $\phi/\Omega_0$ is integer because then we can write

$$g[n]=\cos(\Omega_0(n+\phi/\Omega_0))=f[n+\phi/\Omega_0]=f[n+n_0],\quad n_0\in\mathbb{Z}$$


1) The sampling period ($T_{s}$) is the time between each sample of the observed analog signal that is captured and stored in 'Vreal'

2) The total observation time ($T_{obs}$) is then equal to $N T{s}$ with N being the number of samples you captured

Very roughly, In order to properly monitor your signal, you need $T_s$ to be at the very least 2 times lower than the fastest observed phenomenon.
In other words, you need the sampling frequency ($F_s$) to be 2 times greater than the signal's maximum frequency (see https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem).

Also, you need $T_{obs}$ to be great enough to capture enough heartbeats to the spectral analysis done to determine the heart pulse rate actually makes sense.

What you need to know then, is :

1) How much heartbeats you need to capture so the heartbeat can be properly determined. This will be answered by people working in the field of medicine more than people from signal processing

2) How is a heart beat represented in terms of spectral content, and how large is its spectral bandwidth so you can set up your sampling frequency high enough to be able to capture them without missing critical components


How to measure/estimate symbol rate of a signal with unknown ( but probably QAM ) modulation?

The nonlinear trick is in fact often used to recover the symbol clock of $\frac \pi {2N}$ phase-shifted signals, e.g. BPSK, QPSK: Squaring samples $x[n]$ (or even taking $x^4$) leads to a removal of phase information, and you can then optimize timing with something that is pseudo-derivate based like Müller&Muller.

Now, the problem with QAM is usually that there's no two constellation points with the same phase; so an arbitrary amount of squaring can not reduce the number of phases to 1.

For such signals synchronization is complicated even with exact knowledge of what to expect; it's pretty common to see preamble-correlation based synchronization in receivers. Upside of that is that you can not only get timing, but also a fair amount of CSI that way, which leads to the possibility of equalization.

As a first step, however, some kind of autocorrelation-based detector surely sounds promising: Assuming that with a finite set of constellation points, there's also only finitely many sequences of symbols within the influence of a single pulse shaping filter impulse response, which also means that one might expect the signal to repeat after that amount of time. Hence, a peak on the autocorrelation¹ function of your signal probably signifies the length of $N,\, N\in\mathbb N$ symbols.

With that known, you'd still be on the lookout for $N$; remember that the spectrum of a data-modulated signal is the pulse shape convolved with a diraccomb with an inter-dirac distance of $\frac{1}{f_{symbol}}$. Hoping that symbol timing coincides with the pulse shape filter's time domain roots, you could just look at a sufficiently oversampled signal and count the spectral "peaks"; in reality, this will typically not be as clear.

A few other hints: If the spectral shape is wideband and rectangular, possibly even with a bit of "sincy" sidelobe on either side, it's probably OFDM. That'll make your life a little harder, but: OFDM typically uses cyclic prefixes, which means that you can find the OFDM symbol duration using autocorrelation; since the cyclic prefix needs to be a whole number of subsymbols (IDFT bins), you can narrow down the possible OFDM configurations quite a bit.

¹In fact, any periodicity-based estimator might work; have a look at thegr-specest toolbox. It has some parametric spectrum estimators, of which MTM might be very interesting, here.


Convolution Integral of Harmonic Signal (Cosine) with the Sinc Function

As $T\rightarrow \infty$, the integral becomes the convolution integral. You can use the fact that convolution in the time domain is multiplication in the frequency domain and then you have a single carrier being passed through a low pass filter. The term where $\omega=\omega_c$ just has to do with how you define the ideal low pass filter and in this case you'd say that the value at exactly $\omega=\omega_c$ is equal to $\frac{1}{2}$.


Is there an existing feedback elimination algorithm that can be implemented in iOS?

I've never implemented this myself, but here are my (possibly incorrect) thoughts:

The general method of feedback elimination is to

Detect the frequencies that are feeding back
Tune a filter to eliminate those frequencies

Step 2. is relatively easy, especially if you are using iOS and have an audio unit that can do it for you. There are some things to watch out for -- because you need a fairly steep filter you need to make sure it is stable even at extreme settings -- but hopefully the folks at Apple have taken care of that stuff for you.

Step 1. is trickier -- how do you know something is feedback and not a genuine part of the audio? There are some advanced techniques related to echo cancellation and adaptive filtering for doing this, but I'm not aware of anything widely available as code or libraries. There might be some shortcuts for your specific use-case, e.g. If you are dealing with speech only, you might be able to simply track pitch, and if a single pitch is strong enough assume that it's a result of feedback, but even doing that quickly enough may not be trivial.

Note that this entire process usually needs to be done continuously/adaptively, at least to some extent to make sure movements of the mike and speaker don't cause new feedback.

There are several patents associated with this, mostly because it's an important aspect of designing hearing-aids. Perhaps some of them are expired and you use them:

https://patents.google.com/patent/US4783818A/en?q=method&q=of&q=and&q=means&q=for&q=adaptively&q=filtering&q=screeching+noise&q=caused&q=by&q=acoustic+feedback
https://patents.google.com/patent/US4677677A/en?q=acoustic+feedback&q=cancelation&page=1

(Note, wikipedia mentions another method based on pitch shifting. I'm not familiar with this, and it sounds pretty silly, but it might be worth investigating.)


Find the difference equation and draw the simulation diagram

Your equation is:
$$y[n] = 0.2929 x[n] + 0.2929 x[n−1] + 0.4137 y[n−1]$$
which has input $x[n]$ and output $y[n]$.

You also need $x[n-1]$, so there is a delay block required for that and $y[n-1]$ so a delay block is required to remember the last output.

You also need constant multipliers for the coefficients 0.2929 and 0.4137.

Read up on signal flow graphs.


wavelet transform also serves to you what you want to do, all you have to do is to first apply a low-pass your signal so as to keep the frequency range you want to scan only, another approach is to transform Gabor where you have to define the size and shape of the analysis window, but I would recommend better  using wavelet transform applying a lowpass filter your signal first.


What kind of discrete convolution is this?

Fast (FFT/IFFT) convolution results in circular convolution.

(But if sufficient zero-padding is explicitly added before the initial FFTs, circular and linear convolution would then produce identical results.)


Meaning of a formula about transient frequency bins

without reading the paper,

Eqn 4 has the form of a central difference form of the second derivative. It is a numerical second derivative compared to a threshold. It differs from the standard form in that it centers on $n-1$.

https://en.m.wikipedia.org/wiki/Finite_difference

Whatever $\phi$ is, it is the positive change compared to a threshold.


If you don't know what you're doing, naive bayes is hard to screw up. The naive part of naive bayes is that it assumes that all of your features are independent. They probably aren't independent but in many cases it doesn't matter so much that this assumption is violated. Naive bayes is also good at handling different feature types such as using both categorical variables and continuous variables for the same classification task.


Gesture recognition from analog sensor value

You can use a matched filter. Average all the gesture 1's and that's your filter. Do the same for however many gestures you have and you will effectively have a bunch of filters. Run the data through those filters in parallel and then you will need a decision threshold that indicates yes/no if that gesture is detected.

It should look something like this:


Calculate impulse response when output contributes in input

It should be pretty clear that just computing the root of the characteristic equation can't be enough to find the impulse response. Because in that case the impulse responses of all systems described by
$$y[n]=y[n-1]+\sum_{k=0}^Nb_kx[n-k]\tag{1}$$
would be the same, irrespective of the values of the coefficients $b_k$.
For the given example, the most straightforward (and most instructive) way to determine the impulse response is simply to go through the given iteration until a pattern emerges:
$$\begin{align}h[0]&=1\\h[1]&=\frac32\\h[2]&=\frac32\cdot\frac12\\h[3]&=\frac32\left(\frac12\right)^2\\\vdots\end{align}$$
Can you see any pattern?


Apply LTI filter to incoming signal stream

As a quick draft of an answer, it turns out the things I was looking for was:

Digital Filter Realization
Digital Low Pass Filter Realization


What is the role of a LPF in oversampling?

@user29568: Your question is very far from stupid. It's a darned good question. To see why time-domain "zero stuffing" generates the additional spectral images, that must be filtered out as described by Fat32, you're welcome to see my explanation of this effect at: http://www.dsprelated.com/showarticle/761.php


The confusion comes from the fact that what is tagged as "transmitted" isn't the real transmitted waveform but its baseband representation, which are $\left\{+1,-1\right\}$ symbols in 1 dimension. For 2-FSK, that representation would be $\left\{\left(+1,0\right), \left(0,+1\right)\right\}$ symbols in 2 dimensions.

A common convenient way to handle such 2 dimensional symbols is to represent them with complex number, mapping the coordinate $\left(+1,0\right)$ to the purely real value $+1$ and the coordinate $\left(0,+1\right)$ to the purely imaginary value $j=\sqrt{-1}$. Using that mapping, the symbol generation code for 2-FSK becomes:

if strcmp(tx_type, "fsk")
  tx_bits = round(rand(1,N));
  s0 = zeros(N,1);
  s0M = zeros(N*M,1);

  for i=1:N
    if tx_bits(i)
      s0(i) = 1i;
    else
      s0(i) = 1;
    end
  end

  k = 1;
  for i=1:M:N*M
   s0M(i:i+M-1) = s0(k);
    k ++;
  end
end

The problem then is that the CMA weight update rule from the original post is given for a purely real signal. The adjusted weight update rule for complex signals is given by:

$$
\begin{align}
  y &= \vec w_{n}^H x\\
  e_n &= \left|y\right|^2 - 1 \\
  \vec w_{n+1} &= \vec w_{n} - \mu \vec \nabla e_n \\
           &= \vec w_{n} - \mu' e_n \Re\left\{y^\star x\right\}
\end{align}
$$

This can be implemented with:

for i = 1:N*M-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = abs(y)^2 - 1;
    w = w - mu * e(i) * real(conj(y) * x);
end

Plotting the results (showing the real part and imaginary parts in separate graphs) should give you something like the following:


Types of interpolation used for reconstruction in DSP?

Zero-order hold will result in a piecewise-constant waveform.

Linear interpolation will result in a piecewise-linear waveform.

If you want a piecewise-quadratic or piecewise-cubic or higher order polynomial interpolation, it will not appear much different from the original bandlimited waveform.


Polarity inversion in frequency domain (polar coordinates)

This is really basic, but here's how it works. If $X(\omega)$ is the Fourier transform of $x(t)$, then we obviously have

$$\mathcal{F}\{-x(t)\}=-\mathcal{F}\{x(t)\}=-X(\omega)$$

Furthermore,

$$X(\omega)=|X(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase. From this we have

$$-X(\omega)=-|X(\omega)|e^{j\phi(\omega)}=|X(\omega)|e^{j(\phi(\omega)\pm\pi)}$$

because $e^{\pm j\pi}=-1$. So Jazzmaniac's comment that the phase needs to be shifed by $\pm \pi$ (or any odd multiple of $\pi$ for that matter) is of course correct.


So you have a complex impulse response, and the real and imaginary parts of each tap are modeled as i.i.d. zero mean Gaussian variables with variance $\sigma^2/2$ (correct me if I'm wrong). The DFT of the impulse response is a weighted sum of the filter taps:

$$H[k]=\sum_{n=0}^{N-1}h[n]e^{-j2\pi nk/N}\tag{1}$$

where I assume $N\ge L$ (where $L$ is the number of taps), and in Eq. (1) $h[n]=0$ for $n\ge L$. From (1) it follows that $H[k]$ is also a complex Gaussian RV with zero mean and with variance $L\sigma^2$ (the variance of the real and imaginary part is $L\sigma^2/2$, respectively). From this it follows that the magnitude $|H[k]|$ has a Rayleigh distribution:

$$f_{|H|}(x)=\frac{2x}{L\sigma^2}e^{-x^2/(L\sigma^2)}\tag{2}$$

And, accordingly, the squared magnitude $|H[k]|^2$ has an exponential distribution:

$$f_{|H|^2}(x)=\lambda e^{-\lambda x},\quad x\ge 0\tag{3}$$

With $\lambda=1/(L\sigma^2)$, i.e. with mean $1/\lambda=L\sigma^2$, and variance $1/\lambda^2=L^2\sigma^4$.


The equality
$$\frac{1}{T}\sum_{k = -\infty}^{\infty}X(j(\omega - k\frac{2\pi}{T})) = \sum_{k = -\infty}^{\infty}x(kT)e^{-j\omega kT}\tag{1}$$
is an instance of Poisson's sum formula. The term on the right-hand side of $(1)$ is just the Fourier series representation of the periodic function on the left-hand side of $(1)$.
The samples $x(kT)$ of the time domain signal are basically the Fourier series coefficients of the periodic spectrum of the sampled signal.
The significance of $(1)$ is that it shows that the discrete-time Fourier transform (DTFT) $X_d(e^{j\Omega})$ of the sequence $x_d[k]=x(kT)$ equals a periodized version of the Fourier transform of the corresponding continuous-time signal $x(t)$:
$$X_d(e^{j\Omega})=\sum_{k=-\infty}^{\infty}x_d[k]e^{-jk\Omega}=\frac{1}{T}\sum_{k=-\infty}^{\infty}X\left(\frac{j(\Omega-2\pi k)}{T}\right),\qquad\Omega=\omega T\tag{2}$$
Clearly, if $X(j\omega)$ is band-limited with a maximum frequency $\omega_c<\pi/T$, then the shifted spectra in the sum on the right-hand side of $(2)$ won't overlap, i.e., there is no aliasing and the signal $x(t)$ can be reconstructed without error from its samples $x(kT)$. That means that the basic form of the sampling theorem is implicit in Eq. $(2)$.


What is difference between Signal frequency and data rate

There is a mistake. The units of frequency is cycles/second and not Mbit/second. I noticed this mistake in one of the online links on CAN.

I think what it means is  that the maximum theoretically possible data rate is 1Mbps and the typical real life maximum data rate is 320Kbps(40KBytes/sec). The standard says that for network lengths below 40m, 1Mbps can be achieved but for network lengths of 500m the data rate achieved is around 125kbps. So I think that's reason typical data rate is given as 40KBytes/sec.

References:

http://en.wikipedia.org/wiki/CAN_bus
http://www.computer-solutions.co.uk/info/Embedded_tutorials/can_tutorial.htm


generating white gaussian noise in matlab using two different functions

A white noise sequence is one for which each (random) element is uncorrelated from every other element:
$$
E[y[n]y[m]] = \left \{ \begin{array}{ll}
0 & \mbox{for } n\not=m\\
\sigma_y^2 & \mbox{for } n = m
\end{array}
\right . \\= \sigma_y^2 \delta[n-m]
$$
where $\sigma^2_y$ is the variance of $y$.

Note that I am assuming (because of the whiteness) that the signal is zero mean.

The whiteness of a signal says nothing about the distribution of its values.  To know something about that, Gaussianity or some other distribution needs to be invoked.

The functions wgn and randn both produce white, Gaussian noise sequences.

Calling the function rand would produce a white, uniformly distributed noise sequence.


Property for time invariance of system in terms of its impulse response

Time-invariance means that if $y[n]$ is the response to $x[n]$, then $y[n-k]$ is the response to $x[n-k]$. So if $h_0[n]$ is the response to the impulse $\delta[n]$, then - if the system is time-invariant - its response to $\delta[n-k]$ must be $h_0[n-k]$.

If $h_k[n]=h_0[n]$ were true, then the response to any input $\delta[n-k]$ would be $h_0[n]$, regardless of $k$, which basically means that the system output would be independent of the input signal.


Academic problem solving in Computer Vision

With all due respect, your question is rather vague...  For example, it is not clear to me what you mean by "structured methods of computer vision problem solving", or what you mean by a "graduate manner" of object detection or recognition.

If you are doing a Master's thesis in computer vision then presumably you have an adviser, who is supposed to be knowledgeable in this field, and who is supposed to help you define your problem, and point you to the relevant papers.

By the way, defining and constraining your problem is the key here. "general steps of solving an object detection (or recognition)" is not a very well defined problem.  Are you trying to detect a specific object, such as a particular book or your favorite coffee mug? Or are you trying to detect members of broad categories of objects, such as faces, cars, people, or chairs?

Does your method need to be invariant to in-plane rotation?  The answer is probably "no" for faces, which are rarely upside down, but "yes" for detecting cars in aerial images.  There are many more questions like this.

So, my point is that object detection, let alone recognition, is a broad area, and it is best to talk to your adviser about scoping it down.


transient and steady-state response of first order system

The region of convergence (ROC) of the $\mathcal{Z}$-transforms of a step-modulated sinusoid

$$x[n]=\sin(\omega_0n)u[n]\tag{1}$$

is the region $|z|>1$. The ROC of the transfer function

$$H(z)=\frac{b_0+b_1z^{-1}}{1- z^{-1}}\tag{2}$$

also equals $|z|>1$ (assuming a causal system). Consequently, multiplying the two transforms doesn't pose any problem, and the result also converges for $|z|>1$.

As an example, let's compute the response $y[n]$ of the system $(2)$ to the input $(1)$. For the sake of simplicity we assume $b_0=1$ and $b_1=0$ in $(2)$. First, we compute the response to

$$\tilde{x}[n]=e^{j\omega_0n}u[n]\tag{3}$$

The $\mathcal{Z}$-transform of $(3)$ is

$$\tilde{X}(z)=\frac{1}{1-e^{j\omega_0}z^{-1}}\tag{4}$$

The $\mathcal{Z}$-transform of the response $\tilde{y}[n]$ is given by

$$\begin{align}\tilde{Y}(z)&=\tilde{X}(z)H(z)\\&=\frac{1}{1-e^{j\omega_0}z^{-1}}\cdot \frac{1}{1-z^{-1}}\\&=\frac{A}{1-e^{j\omega_0}z^{-1}}+\frac{A^*}{1-z^{-1}}\tag{5}\end{align}$$

with

$$A=\frac{1}{1-e^{-j\omega_0}}=H(e^{j\omega_0})\tag{6}$$

The inverse $\mathcal{Z}$-transform of $(5)$ is

$$\tilde{y}[n]=H(e^{j\omega_0})e^{j\omega_0n}u[n]+H^*(e^{j\omega_0})u[n]\tag{7}$$

The response to the step-modulated sinusoid $(1)$ is easily obtained from $(7)$ by taking its imaginary part:

$$\begin{align}y[n]&=\textrm{Im}\big\{\tilde{y}[n]\big\}\\&=\big|H(e^{j\omega_0})\big|\sin\big(\omega_0n+\arg\left\{H(e^{j\omega_0})\right\}\big)u[n]\\&\qquad -\textrm{Im}\big\{H(e^{j\omega_0})\big\}u[n]\tag{8}\end{align}$$

The first term in $(8)$ is the steady-state response, and the second term is the transient response, which doesn't decay because of the system's pole at $z=1$.

So the output has a DC component due the imaginary part of $H(e^{j\omega_0})$. Note that this is no contradiction with the linearity of the system, because the input signal is not a single spectral line but it has a continuous spectrum extending down to DC due to the sinusoid being switched on at $n=0$. This DC value of the input's Fourier transform triggers the system's eigenfrequency, which is a DC component. Note that this always happens with LTI systems: the output may contain oscillations at frequencies that are only determined by the system, not by the input signal. However, unlike in the given example, usually these transients decay because most of the time we consider asymptotically stable systems.

I slightly modified your Matlab/Octave script and added the analytical result $(8)$ for comparison. The analytical result is identical (up to numerical accuracy) to the result obtained by filtering:

fs = 48000;
f0 = 1000;
w0 = 2*pi*f0/fs;
N = 200;
n = 0:N-1;
x = sin(w0*n);

bd = [1 0];
ad = [1 -1];
yi = filter(bd, ad, x);
figure;
plot(n,x,'b'); hold on;
plot(n,yi,'r')

% analytical computation
A = 1 / ( 1 - exp( -1i*w0 ) );
y2 = abs(A) * sin( w0*n + angle(A) ) - imag(A);

plot(n,y2,'k.'), hold off
legend('input','response by filtering','response by analytical computation')


Peaks are determined by energy in a given frequency range which can be more or less visible depending on the dimension of frequency bins.


variance in the time domain versus variance in frequency domain

Variance is never defined as power. For a wide-sense stationary random process $X(t)$ with zero mean

$$\mu_X=E\{X(t)\}=0\tag{1}$$

the variance of $X(t)$ equals its power.

The autocorrelation of $X(t)$ is defined by

$$R_X(\tau)=E\{X^*(t)X(t+\tau)\}\tag{2}$$

The power of $X(t)$ is

$$P_X=E\{|X(t)|^2\}=R_X(0)\tag{3}$$

The variance of $X(t)$ is

$$\sigma^2_X=E\{|X(t)-\mu_X|^2\}=E\{|X(t)|^2\}-|\mu_X|^2=R_X(0)-|\mu_X|^2\tag{4}$$

Clearly, if $\mu_X=0$, we have $\sigma^2_X=R_X(0)=P_X$.

Note that the autocovariance of $X(t)$ is defined by

$$C_X(\tau)=E\{[X(t)-\mu_X]^*[X(t+\tau)-\mu_X]\}=R_X(\tau)-|\mu_X|^2\tag{5}$$

Consequently,

$$\sigma^2_X=C_X(0)\tag{6}$$

But since the power spectral density $S_X(\omega)$ is the Fourier transform of $R_X(\tau)$, the integral over $S_X(\omega)$ equals the power of $X(t)$ and not its variance, unless $\mu_X=0$:

$$P_X=R_X(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_X(\omega)d\omega\tag{7}$$

Also take a look at this related answer.

For the difference between the definition of autocorrelation in statistics and in signal processing look here.


What data should be collected to determine the autocorrelation sequence?

All you need to do an autocorrelation is a number of evenly spaced samples at a given frequency, so you will need some dwell time at each frequency (which I assume that you are doing already).  That is all you need.


What is the distance metric in dynamic time warping?

The Euclidean distance between two vectors $\mathbf a$ and $\mathbf b$ is $\mathrm d(\mathbf a, \mathbf b) \equiv \sqrt {\sum_i (a_i - b_i)^2}$

In your example you forgot to square the differences. Remember the way you calculate hypotenuse of a triangle; it's the same idea.


Yes, that's right.

Even more generally: Usually, we call what you call baud rate the symbol rate, and that is simply how many symbols per fixed time we transmit. The unit is baud [bd], and that's why many people call it (imho, slightly incorrectly) baud rate.

These symbols don't have to be on or off, only. If you, for example, had -3V, -1V, +1V, +3V as the four possible values on your line, than you'd have a bit rate (in bit/s) that would be twice as high as your symbol rate (in bd, or Sym/s), because every symbol carries two bits of information (four values = 2 bits of info, because 2bits = number of options).

Going away from the purely theoretical point of view: for longer links, a simple binary voltage communication is often unreliable, and you might very possibly have to look for other kinds of links that are designed to be used in this scenario. I'd recommend finding a way of analyzing your signal at the receiving end, and making sure you properly terminate that.


Try this:

[argvalue, argmin] = min(x);
[argvalue, argmax] = max(x);

Also, have a look to this question.


Assuming you want to remove the DC bias from your signal, a few thoughts come to mind. Your signal sequence looks like a noisy periodic wave riding on a DC bias. And the bias is trending upward in value from the beginning to the end of your sequence. Computing the average value of your signal sequence, a single value, and subtracting that value from your sequence will produce a sequence whose average value is very close to zero. But the upward bias trend will remain in your new sequence.

A process called "DC blocking" will remove both the average bias and the upward bias trend from your signal. If you have Matlab software, then use Matlab's detrend() command to implement DC blocking. If you don't have Matlab available then have a look at the following DC blocking networks (copied from Chapter 13 of my DSP book).

Other than different DC gains those above networks have the same frequency response. (Variable $\alpha$ is in the range of just greater than zero -to- one.) Note, those networks have very nonlinear phase response near zero Hz. If linear phase DC blocking is what you want then have a look at:
http://www.dsprelated.com/showarticle/58.php


MU-MIMO Sensitivity to Line of Sight (LOS)

There is no mention of light-of-sight (LOS) on the time you mentioned. There is a mention of point-to-point, however. Point-to-point channels could be LOS or multipath.

MU-MIMO is better than point-to-point MIMO because users are usually well-separated in space. MU-MIMO is particularly better when the channels are correlated at the base station (BS). The reason for that is that the spatial separation of the users make the channel correlation matrices from different users to be different. In the ideal case where these correlation matrices eigenspaces are orthogonal, the data from different users can be perfectly separated without any interference. In point-to-point MIMO, on the other hand, the data from different antennas from the multiple-antenna transmitter will have the same correlation matrix at the BS, and thus there will be more interference from other data streams.

Read section 2.4 from the book Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency.

EDIT: In the paper, they are talking about separating users that have LOS channels with the BS in massive MIMO, that are difficult to separate in traditional MIMO.