What happens with signal in frequency spectrum when it is time shifted in time spectrum?

Each frequency in the FT of a time shifted waveform is rotated in phase by an amount proportional to the frequency and proportional to the amount of time shift.

If you delay a pure sinusoid by 25% of its period, it's phase referenced to any fixed point in time will change by pi/2 radians.  Delay a slightly higher frequency sinusoid by the same amount of absolute time, and it's phase will change more.  So the phase change will change with frequency.

If looking at a 3d plot of the FT of a time shifted signal, it looks like taking the FT before the time shift and twisting it.  The twist will be linear, e.g. proportional to frequency.  The more time shift, the greater the amount twist per unit of frequency (or twist revolutions per graph width).


Fourier Series representation of a signal

You're getting yourself into unnecessary trouble by choosing the integration limits exactly at those values of $t$ where you have Dirac impulses. Note that you can choose any integration limit as long as you integrate over one period of the given function.

E.g., if you choose some positive $\epsilon$ satisfying $0<\epsilon\le\frac16$ and you integrate from $-\epsilon$ to $\frac13+\epsilon$ then you have the relevant portion of the signal during one period inside the integral, and there's no question as to the scaling factor of the Dirac impulses. The result of the integration will be a real-valued constant for even $k$, and another real-valued constant for odd $k$. I'm sure you can derive the final result yourself.


The point in the article is that you can sample the signal at lower frequencies than 60Hz by subsampling. Consider the following two scenarios of sampling the 60Hz signal:

You sample like Nyquist wanted, at freq $\geq 120\mbox{Hz}$. In one period of your 60Hz signal, or $T=16.67\mbox{ms}$, you will have a correct representation of your signal via sampling.

You sample at a frequency even lower than the $60\mbox{Hz}$ signal freq. The trick is to sample at frequencies that do not divide the main frequency, which means you cannot use $2, 3, 6, 10, 12, 15, 20, \mbox{ or } 30 \mbox{Hz}$ frequencies, because then you'd always get the same point of the signal.

Let's say you sample at $19\mbox{Hz}$ like in the article. The whole point of this is to save the points on your graph, and get the desired signal point-by-point. So you start sampling your signal, get one point, then whole 3 periods of the signal will pass, and you get a second point, which you put on the graph, and so on until you have enough.

How many points do we need? The Nyquist theorem you stated tells us the number of points we need to be sure we got 100% of the signal. $$T_s = (M + \frac{P}{N})\frac{1}{F},$$  where $F$ is our signal frequency of $60\mbox{Hz}$, $T_s = 1/F_s = 1/19\mbox{s} = 52.63\mbox{ms}$ is our sampling period, and $M, N, P$ are integers where $N$ and $P$ must have no common factors. If we try to get such integers, we can see that for $N = 19, M = 3, \mbox{ and } P = 3$ our equation holds. So, we need 19 points for our graph, in order to have the whole signal. This means we need to wait 19 times our sampling period to get our total time of sampling, so we get $$ T_{tot} = 19\cdot 52.63\mbox{ms} = 999.97\mbox{ms} $$

So we see that sampling is possible at even lower frequencies, but at much longer time :) Hope this clarified your question! :)


How to express linear convolution using positive frequencies from channel and symbol DFT vectors?

I think already your fist equation is potentially wrong. The output of the channel is given by LINEAR convolution signal and channel impulse response, not by circular convolution. Circular convolution can be used but only if signal and impulse response of the channel are sufficiently small and enough zero padding is applied to avoid wrap-around errors.

Dumping the imaginary part of the spectrum is also problematic. This basically an "analytic" signal (at least if the negative frequencies are set to 0). This typically is a non-causal operation and it becomes very difficult to determine when you will be getting wrap around errors and when you don't. You will almost certainly end up with a non-causal impulse response of the channel.

What exactly are you trying to achieve and why?


This formula is for baseband with real symbols in an AWGN channel.
However, it is trivially easy to adjust it for quadrature modulation or complex symbols: just double the capacity.
OFDM and other wideband modulations operate over frequency-selective channels, which have different capacity. There are capacity formulas for many channels in the literature.


What is the DFT of DFT of discrete signal

let

$$\begin{align}
 X[k] &= \mathcal{DFT} \Big\{ x[n] \Big\} \\
      &\triangleq \sum\limits_{n=0}^{N-1} x[n] \, e^{-j2\pi nk/N}
\end{align} $$

and

$$ y[n] \triangleq X[n]  $$

(note the substitution of $n$ in for $k$.)  then

$$ Y[k] = \mathcal{DFT} \Big\{ y[n] \Big\} $$

then, if the DFT is defined the most common way (as above):

$$ Y[n] = N \cdot x[-n] $$

where periodicity is implied:  $x[n+N]=x[n]$ for all $n$.


Is sound the best way to transmit a signal under water?

Yes, it is called acoustic communications. Iterative Carrier Frequency Offset and Channel Estimation for Underwater Acoustic OFDM Systems is an example of a paper that uses orthogonal frequency division multiplexing (OFDM) in an underwater acoustic channel.
EDIT:
Note that you wouldn't call it a SONAR any more because SONAR stands for SOund Navigation And Ranging, whereas this is a communication system, just like you wouldn't call your cell phone a radar.
EDIT 2:
Yes, the signal could be converted to a radio signal at a repeater. The repeater would ideally be floating in the surface and have a transducer (which is basically a microphone or array of microphones) on the underwater side to receive the sound waves. The transducer converts the sound waves into an electrical signal, which can be filtered and/or amplified, and then sent to an antenna for radio transmission. The antenna would likely need to be above water because radio waves don't propagate very far through water.


How to evaluate a close-loop control during transient?

In order to evaluate the performance of a controller during transients, you can use different metrics. Two metrics you could use are

ITAE= $$\int_0^{\infty} t*|e| \,dt$$

ITSE = $$\int_0^{\infty} t*e^2 \,dt$$

where e is the difference between your reference and the measurement and t is the time since the reference change.

Basically, you want to minimize these metrics. The basic idea is that an error in the beginning won't impact the metric much while an error 1 second after the reference change will impact your metric much more.

There are probably other metrics you could use. You could also use stabilization time at 1% or 5%.

Edit : If your signal is cyclical you could limit the integral upper bound to a limited number of 2-Hz cycles.

ITAE= $$\int_0^{NT} t*|e| \,dt$$

ITSE = $$\int_0^{NT} t*e^2 \,dt$$

where N is a whole number and T is the period of your reference signal.


Confusions regarding meaning of "kernel" in signal processing and image processing?

What is the meaning of "kernel" in context of filters?

A "kernel" is the thing you convolve by.  So for a FIR filter that's defined by the taps $k = \begin{bmatrix}a_0 & a_1 & a_2 & \cdots\end{bmatrix}$, $k$ is the kernel.

For a 2D filtering problem, the kernel would, in general, be 2D (and you tend to see the term much more often in image processing, but I've also seen it in really highfalutin' stochastic signal processing literature).

  What is meant by "neighborhood of operation" here?

For an output pixel at position $\begin{bmatrix}x, y\end{bmatrix}$, all of the pixels in the input image that affect the output.  In general, if you have an $n \times m$ kernel, the neighborhood of operation is $n \times m$.

  Is it true even in case of FIR windows such as Hanning window,Blackman window etc? Those windows are also kernels?

Different disciplines, different terminologies.  No, because those don't define filter kernels -- the "windowing" operation is different from the "windowing" in 2D image processing.


In the code the OP is comparing the real and the imaginary components of the FFT in each to determine similarity. If a phase shift occurs then the real and imaginary components will also change. What the OP should be comparing is the magnitudes of the complex values in each sample, which should be invariant under such image shifts.
Confirming this; the OP had done the following to compare the FFT results:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(H1.imag-H1.imag)),np.mean(np.square(H1.imag-H2.imag)),np.mean(np.square(H1.imag-H3.imag))])

resulting in:

If one does the following instead:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(np.abs(H1)-np.abs(H1))),np.mean(np.square(np.abs(H1)-np.abs(H2))),np.mean(np.square(np.abs(H1)-np.abs(H3)))])

We get a result that we expect:

The equivalent magnitude is given by the translation property of the 2D DFT where:
$$x[k-m, l-n] \mapsto X[u,v]e^{-j2\pi (mu/M + nv/N)}$$
This is similar to the time shift property of the 1D DFT (which shows fundamentally what is occurring and can be extended to a 2D DFT):
Given 1D DFT as
$$x[n] \mapsto \sum_{n=0}^{N-1}x[n]W_N^{nk} = X[k]$$
Where $W_N^n=e^{-j2\pi n/N}$
From which:
$$x[n-m] \mapsto \sum_{n=0}^{N-1}x[n-m]W_N^{nk}$$
Substitute $u= n-m$
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{(u+m)k}  = \sum_{u = -m}^{N-1-m}x[u]W_N^{uk}W_N^{mk}$$
$$=W_N^{mk}\sum_{u = -m}^{N-1-m}x[u]W_N^{uk}$$
for any k:
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{uk} = \sum_{u=0}^{N-1}x[u]W_N^{uk}$$
Therefore
$$x[n-m] \mapsto W_N^{mk}X[k]$$
And $W_N^{mk}$ as $e^{-j2\pi mk/N}$ only modifies phase, not magnitude
Also the DFT is circularly periodic in the time domain, so any such rotation must done in circular fashion;  x[mod(n-m,M)] = x [n]


How to compare binary signals with different lengths

you need to elaborate on what you mean by similarity because there are many possibilities.

If you key sequence is embedded in any of your test sequences, the cross correlation is as was suggested by one of the commentators is a good approach, which forms a point by point product between the test and key sequence and sums the results of the product and provides a number, which if the key has $n$ ones the sum of products equal to $n$ indicates an exact match.  The key sequence is shifted relative to the test and the sum of products is calculated for all possible shifts.  This also works for situations where there are error bits in the tested sequence, when the sum product is near $n$.

One might presume that a sum product equal to zero is an indication of dissimilarity, and could be depending on your application but it is usually preferable that measures of similarity should be invariant under information preserving transformations of the data. The replacement of each one by a zero is an information preserving transformation, so a cross correlation value near zero in this case should be considered to be very similar.

A reverse ordering of test sequences is another information preserving transformations.

It is usual to consider dissimilarity as a situation where  a one in the reference would correspond to either a one or zero in the other. In this case dissimilarity to equivalent to statistical independence.

It is sometimes useful to transform $(0,1)$ to $(-1,1)$ when calculating cross correlation.

You can also compare sequences by the ratio of ones to zeros. Run length statistics are also useful.

Cross correlation is compatible with statistical hypothesis testing and confidence bounds. These calculations are covered in any good calculus level statistics text book (like Bhattacharyya and Johnson)


Want to do shot boundary detection via SVM, what are some good features?

Basically, Shot Boundary detection technique's output is affected by camera motion, object motion and illumination conditions. The feature you choose should be invariant to these issues. The problem in your case may be due to poor training dataset/ imbalanced dataset. threshold can be fixed experimentally and it is tedious.


Why is `fftfilt` (i.e. `fft` of both inputs, then element-wise multiplication, then `ifft`) faster than direct convolution?

Whether the direct convolution or the FFT/IFFT method is faster depends on the length of the impulse response, $N_\mathrm{i}$ and the signal length $N_\mathrm{s}$. With the formulas taken from here I've created a small Matlab script that calculates the required number of real multiplications and additions for the direct convolution and FFT/IFFT method, respectively. I've assumed real-valued impulse response and signal. Taking your example of $N_\mathrm{i}=32$ (or the other way round, this makes no difference), the direct convolution always requires less multiplications:

But if we increase $N_\mathrm{i}$ to 50, for example, the direct convolution requires less multiplications up to a signal length of about 60. For signal length greater than 60 the FFT/IFFT method requires less multiplications:

The actual computation time depends on some other parameters, especially the CPU architecture but the number of multplications is usually a good indicator.

The MATLAB code for the above figures for reference:

Ns = 1:1000; % length of signal
Ni = 50; % length of impulse response

K = Ni + Ns - 1;

% number of real multiplications for convolution:
M_R_conv = Ns*Ni;
% number of real multiplications for convolution:
A_R_conv = Ns*(Ni-1);

% number of complex multiplications for freq domain conv:
M_C_dft = 3/4 * K .* (log2(K) + 1) + K;
% number of complex additions for freq domain conv:
A_C_dft = 3/2 * K .* log2(K);

% number of real multiplications for freq domain conv:
M_R_dft = 4*M_C_dft;
% number of real additions for freq domain conv:
A_C_dft = 2 * A_C_dft + 2 * M_R_dft;

figure('Position', [0 0 400 300]);
plot(M_R_conv, 'b');
hold on;
plot(M_R_dft, 'r');
hold off;
legend('direct', 'fft', 'location', 'Southeast');
title(['number of real multiplications. Ni: ' num2str(Ni)]);
xlabel('length of signal');

A sidenote to your code: as it stands it is not running as you're creating NxN and MxM matrices for x and y. What you probably want to do is

x = rand(1,N);
y = ones(1, M)
% ...
z = ifft(fft(x, K).*fft(y, K));

p.s. for those of you interested in Python code instead of Matlab, here is the same code in Python: https://gitlab.com/snippets/1789085


Is there a way to measure each pixel's content in an image somehow?

I would want to start from e.g. making a measure of "color" that's a single real value, rather than a triplet. Then I could perhaps want to make a function that calculates the "change in" color between pixels adjacent to each other. So a derivative of some kind

Why would you need the color to be a single value for that? I'll do an analogy: To measure the distance between to cities, you don't need to map these two cities onto the real numbers; they still have coordinate vectors. All you do is map the difference vector onto a number. And for numbers, it's just the same: It's mathematically impossible to sensibly map them onto a single dimension and still preserve distance. So, you can't.

What you need is just some kind of a norm to map a color difference vector to a "perceptive difference". For many color systems, something like the squared sum of vector components would be a good start.


Given LTI system with transfer function :
$$H(s) = \frac{ -2 s}{(s+6)(s+2) } \tag{1}$$
which has the poles $s = \{ -2 , -6 \} $, and zeros $ s = 0 $.
Apply partial fraction expansion :
$$H(s) = \frac{1}{s+2} - \frac{3}{s+6} \tag{2}$$
Based on pole locations, there will be three ROCs,

ROC1 :  -$2 < \mathcal{Re}\{s\} < \infty $  ; stable
ROC2 : $-\infty < \mathcal{Re}\{s\} < -6$  ; unstable
ROC3 : $ -6 < \mathcal{Re}\{s\} < -2$   ; unstable

The impulse responses found from (2) and the chosen ROC , with the fact that
$$\mathcal{L}\{ e^{at}u(t) \} \longleftrightarrow \frac{1}{s-a} ~~,~~ \mathcal{Re}\{s\} >\mathcal{Re}\{a\} \tag{3}$$
$$\mathcal{L}\{ -e^{at}u(-t) \} \longleftrightarrow \frac{1}{s-a} ~~,~~ \mathcal{Re}\{s\} <\mathcal{Re}\{a\} \tag{4}\\\\$$

ROC1 :  -$2 < \mathcal{Re}\{s\} < \infty $  ;  right-sided  : $h(t) = (e^{-2t}-3 e^{-6t}) u(t)\\\\$
ROC2 : $-\infty < \mathcal{Re}\{s\} < -6$  ; left-sided : $h(t) = (-e^{-2t}+3 e^{-6t}) u(-t)\\\\$
ROC3 : $ -6 < \mathcal{Re}\{s\} < -2$   ; two-sided : $h(t) = -e^{-2t}u(-t)-3 e^{-6t}u(t)$.


Important differences between BCJR and SOVA algorithms?

For an AWGN channel, the exact computation of the log-likelihood ratio requires knowledge of the SNR. However, in Hagenauer's original paper on the SOVA there is an approximation which allows the implementation of the algorithm without estimation of the SNR.


window function between program and math

The fact that you are using a time limited sample of the waveform is the rectangular window that is referenced. Since the fft operates on a time limited sequence, it is mathematically identical to a infinitely long duration of that sequence multiplied the rectangular window.


Help with obtaining the power spectral density of a simple continuous cosine (using both forms of the definition for PSD)

There are several misconceptions in the question that have not been addressed in the existing answers. First of all, the signal $x(t)=A\cos(2\pi f_0t)$ is a deterministic power signal (unless $A$ or $f_0$ are modeled as random variables). For this reason several definitions in the question are inappropriate. First, the auto-correlation of a power signal is given by

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^Tx^*(t)x(t+\tau)dt\tag{1}$$

The integral given in the question (with infinite limits and without division by $T$) does not exist for the given $x(t)$. With definition $(1)$, the auto-correlation of $x(t)$ is indeed obtained as

$$R_x(\tau)=\frac{A^2}{2}\cos(2\pi f_0\tau)\tag{2}$$

The Fourier transform of $(2)$ results in the power spectrum of $x(t)$.

The power spectrum can also be computed directly from $x(t)$, but the formula given in the question only applies to random signals, but not to deterministic signals. For deterministic signals, the appropriate definition is

$$S_x(f)=\lim_{T\to\infty}\frac{1}{2T}\left|\int_{-T}^{T}x(t)e^{-j2\pi ft}dt\right|^2\tag{3}$$

The computation of $(3)$ for the given signal is discussed in this question.


The spectrum of your windowed sinusoid is just a shifted sinc function (assuming a rectangular window). If $T=0.015s$ is the length of the signal then the magnitude of the spectrum at positive frequencies is given by

$$S(f)=T|\text{sinc}((f-f_0)T)|,\quad f>0$$

with $f_0$ the frequency of the sinusoid, and where I use the definition $\text{sinc}(x)=\sin(\pi x)/(\pi x)$. So with $\Delta f=100Hz$ (the distance to the other sinusoid), you get a relative attenuation of

$$|\text{sinc}(\Delta f\cdot T)|=0.21221\doteq 13.5\,\text{dB}$$

So at the frequency of the neighboring sinusoid, the spectrum of the other (finite length) sinusoid is about $13.5\,\text{dB}$ lower.


You are right, but you have mistaken the DTFTs...

Considering the Parseval's relation for real discrete-time sequences to be:

$$ \sum_{n} x[n]y[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega})Y(e^{j\omega}) d\omega $$

We shall define $x_1[n]$, $x_2[n]$ and find their corresponding DTFT as:

$$ x_1[n] = \frac{\sin(\frac{1}{4} \pi n)}{2\pi n} \Leftrightarrow X_1(e^{j \omega}) = \begin{cases} \frac{1}{2} &,& \text{ for } |w| < \frac{\pi}{4} \\ 0 &,& \text{otherwise} \end{cases} $$

$$ x_2[n] = \frac{\sin(\frac{1}{6} \pi n)}{5\pi n} \Leftrightarrow X_1(e^{j \omega}) = \begin{cases} \frac{1}{5} &,& \text{ for } |w| < \frac{\pi}{6} \\ 0 &,& \text{otherwise} \end{cases} $$

Where we have used the DTFT pair;
$$ x[n] = A\frac{\sin(\omega_c  n)}{\pi n} \Leftrightarrow X(e^{j \omega}) = \begin{cases} A &,& \text{ for } |w| < \omega_c \\ 0 &,& \text{otherwise} \end{cases} $$

Hence applying Parseval's theorem we see that:

$$ S = \sum_{n} \frac{\sin(\frac{1}{4} \pi n)}{2\pi n} \frac{\sin(\frac{1}{6} \pi n)}{5\pi n} = \frac{1}{2\pi} \int_{-\pi/6}^{\pi/6} \frac{1}{2}\frac{1}{5} d\omega = \frac{1}{2\pi} \frac{2\pi}{6} \frac{1}{10}  = \frac{1}{60} = 0.01667$$

To check this result you can use the following OCTAVE/matlab code for getting the sum in the time domain:

L = 128;    % select large enough L for better approximation
n = -L:L;
x =  sin(pi*n/4).*sin(pi*n/6)./(2*5*pi*pi.*n.*n);
x(L+1)= (1/4)*(1/6)/(2*5);
S = sum(x)
>> ans = 0.01667

which verifies that our calculation was correct.


Hilbert transform too large to store (out of core processing)

I would use a linear phase FIR Hilbert transformer, and use block processing, such as the overlap-add method. That means that you partition the input signal into contiguous non-overlapping blocks and compute the convolution of each block with your filter impulse response. The results are then overlapped and added. Overlap occurs because the result of the convolution is longer than the input block. E.g., if the impulse response has length $N$ and  the blocks have length $L$, the result of the convolution is $N+L-1$, and the last $N-1$ output samples must be added to the first $N-1$ samples of the next output block (assuming that $N<L$, i.e., the output block only overlaps with the next block).

In short:

$$y[n]=\sum_ky_k[n-kL],\qquad y_k[n]=(x_k\star h)[n]$$

where $h[n]$ is the impulse response, $x_k[n]$ is the $k^{th}$ input block, $y_k[n]$ is the convolution of $h[n]$ with $x_k[n]$, and, finally, $y[n]$ is the output signal.


Is it customary to correct for the gain of a window?

Yes, it is customary to correct for the gain of a window, except for some cases I refer to later. (If you are interested only in the relative amplitude, of course you do not need to correct for the gain.)

Because the window reduces the gain of the original signal (time domain), the amplitude obtained through FFT need to be corrected. For example, if you use the Hanning window, you need to multiply all the amplitudes by 2 (the reciprocal of 0.5). As I understand it, most of the software packages for FFT automatically correct for the window used.

However, such correction is good only when all the frequencies of interest distribute throughout the time domain window. For example, suppose you have 1024 data with all signal levels of zero except for #512 point which has a value of 1 (impulse signal). Obviously, any windows do nothing to the data. So, if you correct the amplitudes for the window gain (multiply by 2), then you will end up with overestimation of the amplitude. If your 1024 data are all zero except for the very 1st point with a value of 1, then every point has a value of zero after windowing, and you lose the signal.

So, if you are dealing with random signals, with all the frequency components expected to lie almost evenly over the length of the signal, you need (or should) correct for the gain of the window you use.


Determine if $ y[n] = ny[n-1] + x[n]$ is linear time invariant and BIBO stable

I'll give you some hints that hopefully will allow you to do your homework yourself.

In order to check linearity, write down $3$ difference equations: one for the response $y_1[n]$ to an input signal $x_1[n]$, one for the response $y_2[n]$ to an input signal $x_2[n]$, and one for the response $y_3[n]$ to an input signal $ax_1[n]+bx_2[n]$. Multiply the first difference equation by $a$, the second one by $b$, and add them. If the resulting difference equation for the sequence $ay_1[n]+by_2[n]$ is the same as for $y_3[n]$, then you can conclude $y_3[n]=ay_1[n]+by_2[n]$, and, consequently, the system must be linear.

For checking time-invariance, write down the difference equation for the response $y_2[n]$ to the input $x[n-k]$. Then replace in the original difference equation the index $n$ by $n-k$ and check if the two difference equations are the same. If they are, the system is time-invariant.

To answer the question about BIBO stability, ask yourself what the output signal will be for $x[n]=\delta[n]$ and draw your conclusion.

Now that you've figured out the solution I would like to point out that you can also find an explicit (non-recursive) expression for $y[n]$:

$$\begin{align} y[0]&=x[0]\\
y[1]&=1\cdot x[0]+x[1]\\
y[2]&=2\cdot 1\cdot x[0]+2\cdot x[1]+x[2]\\
y[3]&=3\cdot 2\cdot 1\cdot x[0]+3\cdot 2\cdot x[1]+3\cdot x[2]+x[3]\\
&\vdots\\
y[n]&=\sum_{k=0}^n\frac{n!}{k!}x[k],\qquad n\ge 0\tag{1}
\end{align}$$

from which it is very straightforward to show that the system is linear and time-varying. Furthermore, from $(1)$ it should be clear that the system is not BIBO-stable.


Calculation of sampling rate and its effect on an Image?

Depends on how the image is created.  If by photography, the focus cone and Airy disk for the given lens have a somewhat low-pass filter effect.  Figure out the effective transition frequency at the imaging plane, and have a sampler with a bit over twice that resolution to diminish aliasing effects during image capture.

If you undersample (either during capture or with later processing) without low pass filtering, then you may get alias artifacts such as Moire patterns on high contrast fine detail.  Classic example is a pin striped shirt at SD TV resolution.


choosing an input signal for system identificaiton

conceptually, any driving signal with bandwidth that covers the maximum expected bandwidth of the system being identified is sufficient.  you drive input and synchronously measure input and output, use the FFT to compute the spectra of both input and output, divide the output spectrum by the input spectrum, inverse FFT the resulting spectrum and you have the impulse response.  at least if nothing bad happened, like the system went non-linear because the driving signal exceeded some bounds of saturation.  or if the amplitude of the driving signal was so low that the result of the division has measurement and numerical errors.

so the whole idea of a good driving function is one that does not exceed limits of linearity in the device, yet delivers a lot of energy to reduce measurement and numerical error.  a measure of efficacy for a driving signal is crest factor, which is the ratio of the maximum signal amplitude to the r.m.s. amplitude.  the lower the crest factor, the better.  it means more energy driving the device under test while remaining within maximum limits.

the impulse driving function is natural in that what is measured at the output is the impulse response.  making the nascent impulse thin enough to have wide bandwidth might mean that the input energy is small.  you would need to measure the impulse response several times by driving the system with a train of impulses, spaced apart sufficiently (longer than any expected impulse response) and synchronously average the results.  but the impulse train will be zero most of the time, and if you limit the impulse height to prevent saturation, then the driving signal has low power and your wait is longer in collecting data.

a single sine wave is not sufficiently broadbanded.  it will tell you only how that system behaves at the frequency of the sine.  a collection of added sinusoids at various frequencies will give you a better broadbanded result, but might add up to a spiky driving signal if all of the phases are the same.  scrambling the phases can help, but finding the best "random" phases can be ad hoc.

a swept-frequency sinusoid is a popular driving signal and there is some theory behind it, if the swept frequency varies linearly with time.  if log-frequency is swept linearly with time, the preceding theory no longer applies.  the crest factor is $\sqrt{2}$.

another possible driving signal is a bipolar maximum-length sequence, which is $(-1)^{a[n]}$ where $a[n]$ is a binary {0,1} maximum-length PN sequence.  with that signal the crest factor is as low as it can get (which is 1), but there can be unpredictable behavior if there are sufficient non-linearity in the system being tested.


O-QPSK constellation interpretation - GNU Radio

The phase recovery was finaly perform in Matlab. Because the signal we are studying (Zigbee) send a preamble before the communication start, we can compensate for the group delay and the channel's distortion with a anti-causal analysis of the expected signal and the received signal.

The GNU-Radio's blocks  weren't efficient at recovering O-QPSK's signal. However with fact that we don't real time treatment of the signal, Matlab was finaly most suitable way to resolve the O-QPSK demodulation of the Zigbee signal.


There are lots of FFT algorithms, and they don't all require power of 2 sizes.  For instance, you can do any size (including prime sizes like 7) with Bluestein's (Chirp-Z) algorithm, and the algorithm is still fast ($N \log(N)$).

I've heard that MATLAB uses FFTW under the hood, and FFTW includes many (all?) of these other algorithms.  So MATLAB can do a length 7 FFT without rounding it to a power of 2.


The Gradient Operator of a Vectorized Image in Matrix Form

It is pretty simple to create those Matrices.
The real issue with them is their size which is enormous for real world images.

For small kernels they are sparse which saves the day.
Indeed for the Derivative Operator, which has only 2 elements, they are highly sparse.

I built them in MATLAB using:

mI = im2double(imread(imageFileName));
mI = mI(11:410, 201:600, 1);

% mI = mI(11:20, 201:210, 1);

numRows     = size(mI, 1);
numCols     = size(mI, 2);
numPixels   = numRows * numCols;

mO = zeros([numRows, numCols, length(vParamLambda)]); %<! Output

vDx = [1, -1]; %<! Matrix is doing Correlation, this for Convolution
vDy = [1; -1];

% Sanity Check
mIxRef = conv2(mI, vDx, 'valid');
mIyRef = conv2(mI, vDy, 'valid');

% mIx = reshape(mDx * mI(:), numRows, numCols - 1);
% mIy = reshape(mDy * mI(:), numRows - 1, numCols);

mDh = sparse(numPixels - numCols, numPixels);
mDv = sparse(numPixels - numRows, numPixels);

% mDx = zeros(numPixels - numCols, numPixels);
% mDy = zeros(numPixels - numRows, numPixels);

tic();
colShift = 0;
for ii = 1:(numPixels - numRows)
    if(mod(ii + colShift, numRows) == 0)
        colShift = colShift + 1;
    end
    mDv(ii, ii + colShift) = -1;
    mDv(ii, ii + colShift + 1) = 1;
end

for ii = 1:(numPixels - numCols)
    mDh(ii, ii) = -1;
    mDh(ii, ii + numCols) = 1;
end
toc();

mIx = reshape(mDh * mI(:), numRows, numCols - 1);
mIy = reshape(mDv * mI(:), numRows - 1, numCols);

mE = mIy - mIyRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Vertical Derivative) - ', num2str(max(abs(mE(:))))]);

mE = mIx - mIxRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Horizontal Derivative) - ', num2str(max(abs(mE(:))))]);

If you load image into mI you see the result is identical to the convolution operator.

The full MATLAB code which implements the solution is at my Signal Processing StackExchange Question 50329 - GitHub Repository.

Remark
In practice usually those kind of equation are solved using Preconditioned Conjugate Gradient Method.
The nice trick is that method requires only the result of the operators applied on a vectorized image. This can be done using the convolution operator itself instead of using the matrix form.


Video stitching for static scene (infinite loop video)

Lets assume your camera axis has polar and azimuthal angles and you want to align your frames only regarding the azimuthal angle and set the polar angle free.

First you have to find the geometric transform between sequential frames of your video, then by multiplying these sequential geometric transforms between your reference frame and the current frame you could obtain the geometric transform between your current frame and the reference frame. To obtain geometric transform between sequential frames you could use automatic point matching algorithms like SIFT or SURF (considering all your frames don't have any intersections but sequential frames have intersection I suggested to use it sequentially).

After finding the geometric transform you have to modify the geometric transform to set the polar direction of camera free. To do this You have to decompose the geometric transform to find external camera parameters which contain the position and the pose of camera. Pose of camera is a Rotation matrix which rotate the axis of camera between your reference and current frame. After obtaining rotation matrix you have to decompose the rotation matrix to obtain the rotation angles which are yaw, pitch and roll and set Yaw angle to zero (I'm not sure which one is the polar angle of your camera but any of these angles which is polar angle of camera must be zero to set the polar direction free), then you have to obtain the modified geometric transform from the modified rotation matrix.

Now you could align all your frames using the inverse geometric transform.

All parts of suggested method have intrinsic MATLAB functions, also if your camera has nonlinear distortion like fish-eye you have to take care of that before point matching part.


Auto-correlation of C/A code with 30dB noise

GPS has a processing gain of 43 dB assuming you correlate over 20 consecutive symbols. (GPS sends 20 C/A code sequences for 1 data bit: Chip rate 1.023 MHz, data rate 50 bps).  It looks like you are correlating over just one symbol which is at 1 KHz in which case the processing gain would only be $10Log_{10}(1023)= 30 \text{ dB}$.

However I don't think your noise is 30 dB stronger than your pre-correlation signal. The standard deviation for a uniform distribution with a range R is $\sigma = R/ \sqrt{12}$. If your range is 30, then the standard deviation would be 8.66. The DC offset does not effect standard deviation. The standard deviation for the C/A code that is +1/-1 is 1, so the predicted pre-correlation SNR would be 18.75 dB. In which case we should see peaks given an expected post-correlation SNR of 11.25 dB when correlating over 1 ms (1023 chips).

Here is the result I got from Matlab with the same test parameters (uniform noise 0 to 30, 1023 chips GPS C/A code +1/-1, correlated over one 1023 chip PRN sequence), where I also shifted 0 lag to be in the center of the plot.

Comparing my plot to yours I think your results are fine: I see a larger peak on the left side of your plot that is reasonably the expected correlation peak: Given the mean of the added noise is 15, the correlation and noise will be shifted up by 15, such that the mean of the peak would be 1023+15 = 1038. The peaks will vary by the post-correlation noise, which would have a standard deviation of  $8.66\sqrt{1023}= 277$ and be now Gaussian distributed due to the law of large numbers. Therefore 68% of the peaks would be expected to be in the range of 1038 +/- 277.

Cross Correlation of SV6 with SV6 + Noise over 1 ms

Added Noise

To Ben's point in the comments, this is a LOT easier to prototype from Matlab/Octave or Python:

Matlab/Octave Code:
(C/A Code generator cacode.m and further details on its implementation are available on this post GPS Coarse Acquisition PRN Codes)

s = cacode(6);  % C/A for SV6
s1 = s*2 -1;    % map to +1/-1
noise = rand(1,1023)*30;    % uniform noise
% circular cross correlation:
xcorr = ifft(fft(s1).*conj(fft(s1+noise)));
figure;
plot(real(xcorr));


What are some typical lossless compression ratios?

It really depends on how much redundancy there is in the data. If all 250GB was just '0', then you could get fabulous levels of compression.

This page shows some results for compression of English text. It compresses 2,988,578 bytes worth of text using various techniques. The top three are: 330,571 (88.94%), 333,759 (88.83%), and 352,722 (88.20%).

Relating this back to your case: that means 250GB would compress (at best) to 27.65GB.

So perhaps there are lots of blank, repetitive patches in the data?

UPDATE:

Based on your edit, here is another paper that looks at XML compression rates.
The upshot (Figure 8; screenshot below) is that they appear to see compression rates ranging from 1 (none) to about 50 at a maximum.  That suggests that you are correct to be suspicious of your customer's 62.5 to 1 (250 to 4) compression rate.


Is a filter/control transfer function with positive phase "causal"?

In general, from the positiveness of the phase nothing can be concluded about the causality of the corresponding system. Note that the phase of the given system is

$$\phi(\omega)=\arctan(\alpha\omega\tau)-\arctan(\omega\tau)\tag{1}$$

which is positive for all $\omega>0$ if $\alpha>1$. Yet the system is causal, which can be easily seen by calculating its impulse response:

$$\begin{align}G(s)&=\frac{\alpha s+\frac{1}{\tau}}{s+\frac{1}{\tau}}=\alpha\frac{s+\frac{1}{\tau}-\frac{1}{\tau}}{s+\frac{1}{\tau}}+\frac{1}{\tau}\frac{1}{s+\frac{1}{\tau}}=\alpha+\frac{1-\alpha}{\tau}\frac{1}{s+\frac{1}{\tau}}\\\Longleftrightarrow g(t)&=\alpha\delta(t)+\frac{1-\alpha}{\tau}e^{-t/\tau}u(t)\tag{2}\end{align}$$

where $u(t)$ is the unit step function. Obviously, we have $g(t)=0$ for $t<0$, i.e., the system is causal.

Note that the system's phase is not only positive, but it also has a positive slope for $\omega<1/(\tau\sqrt{\alpha})$. This means that its group delay (the negative derivative of the phase) is negative in that frequency range, as shown in the figure ($\alpha=2$, $\tau=1$):

So even a negative group delay in some frequency range says nothing about the causality of a system.

A rational transfer function with all its poles in the left half-plane (as is the case in your example) and a region of convergence $\text{Re}\{s\}>0$ always corresponds to a causal and stable system, even if its phase is positive or its group delay is negative in a certain frequency range.


How does delta modulation work practically?

One key point to note about delta modulation schemes is that they quantize the change between the current value and the previously quantized value.

The other key point to note is that, in general, delta modulation schemes are run at sampling rates MUCH higher than the Nyquist rate would suggest. Early portable CD players used a scheme called MASH (multi-stage noise shaping) on their DACs to allow cheaper, more accurate reproduction of sound.

I can't find a reference quickly, but I seem to recall that for the 44.1kHz CD sampling rate, MASH 1-bit DACs ran at sampling rates in the order of 5-6 MHz.


Locking on to a square wave signal with minimum oversampling

One partial answer to my own question is to use the Goertzel algorithm instead of a PLL. The Goertzel algorithm makes it possible to trivially determine the power of a particular frequency; it won't provide the phase, however.

The Goertzel algorithm can be implemented quite trivially using a multiplication and two additions per sample; the sampling rate need only be a minimum of twice the frequency of the signal to be detected. Determining the amplitude of the signal from the Goertzel algorithm's output requires four more multiplies and two more additions.

Since the multiplication in the main part of the Goertzel algorithm is by a constant factor of 2 * cos(2 * pi * w), where w is the number of samples per cycle, one multiplication can be eliminated entirely by sampling at 6 times the target frequency, making 2 * cos(2 * pi * 1/6) = 1.


Normalized cross-correlation in frequency domain

I was thinking on the math behind it and concluded that calculating the NCC in the frequency domain is basically the same process as in the time domain:

$F[k]$ is the (discrete) Fourier transform of my sample, $G[k]$ is the gaussian curve in the frequency domain. The cross-correlation is
$$CC = \sum^{L-1}_{k=0}F[k]{\cdot}G[k]$$
Then, to get the normalized cross-correlation we first normalize the vectors $F$ and $G$, as would be done to any vector:
$$\hat{v} = \frac{v}{|v|} = \frac{v}{\sqrt{\sum^{L-1}_{k=0}v^2[k]}}$$
Thus, the normalized cross-correlation should be the dot-product of the normalized vectors $\hat{F}$ and $\hat{G}$:
$$NCC = \frac{\sum^{L-1}_{k=0}F[k]{\cdot}G[k]}{\sqrt{\sum^{L-1}_{k=0}F^2[k]\cdot\sum^{L-1}_{k=0}G^2[k]}}$$


Why wavelets based transmitter/receiver digital signal processing aren't common?

One key difference between the DFT and discrete wavelet transform (DWT) is that the DFT deals with frequency whereas the wavelet transform deals with "scales". (e.g. you estimate the width of a pulse or a series of pulses in a time series).

A radio signal is more or less a continuous signal when observed locally. A constant waveform signal is very sparse when decoded by FFT. So it is very easy and natural to represent an RF signal using the DFT/frequency domain.

One of the common issues one have to deal with in regards to the DFT is spectral leakage. Unfortunately, the wavelet kernels generally have a good amount of spectral leakage also so this issue cannot be improved by using the DWT.

The WT does have some useful properties. Aside from the examples you listed, there are some applications in time series analysis where the WT trumps the DFT. For example, detecting drift of atomic clocks or measuring changes in tides (see Wavelet Methods, by Percival and Walden). Also, a more limited application: the WT can be used for denoising a signal when you want to preserve the spikes as well as the baseline drift (e.g. to visually make a spectrum look cleaner).

In terms of time-frequency analysis, the Gabor wavelet is probably the closest you get to STFT.

Like the DFT/FFT the DWT is an orthogonal transform. However, the continuous wavelet transform (CWT), which is computed using convolution, is not. It's useful as an analysis tool only.


Apply FIR filters to Extract Stego from WAV data

I never worked with Python or SciPy, but I don't see the sense in the actions you apparently perform in your code.
firwin provides you with the time-domain filter coefficients, which are to be convolved with the time-domain signal you wish to filter.
What you do in your code, is convolve the time-domain filter coefficients with the frequency-domain signal - which makes absolutely no sense.

What is the point of converting the signal to freq-domain (doing FFT)? Is the stego analysis/decoding performed in time- or frequency-domain?

If your processing chain requires conversion of the signal to freq-domain, for some purpose, then you don't really need the FIR filter, since you can get rid of the undesired frequency elements once you have the signal in the frequency domain, by simply zeroing-out the undesired frequency bins of the signal.

If the only reason for calculating FFT of the signal is to filter it - then there is no need to do it, as you can simply convolve the time-domain filter coefficients with the time-domain signal, as I explained.


Should the order of a filter be proportional to the sampling rate?

Keeping the FIR filter the same order and doubling the sample rate will effectively halve the accuracy of frequency resolution which will be noticable at low frequencies (this is based on the assumption you are using something like FFT to create your FIR, otherwise the order of the FIR and it's resolution is a bit more complicated).


The typical method is to use the cross correlation function to recognize the start of each pulse sequence. Typically pulse patterns are chosen for good autocorrelation properties for this very purpose (likely the start of the packet as a header if there is a random data payload also involved). "Good autocorrelation properties" refer to an autocorrelation that is very strong when the sequence is aligned to itself but very weak for any other delays between the sequence and itself (white noise has this property). If the sequence correlates to delayed versions of itself, it will be difficult to differentiate in time when the sequence was received or to differentiate from separate sensors if they all use the same sequence.


It seemed to work relatively well where there was a strong lightness difference at the edge of the poster, compared to background.

This matches what squares.cpp seems to do: scale the image, do a canny edge detection, simplify the resulting contours and look for squares.

There's two things that, looking at your example images, go wrong here:

The contour detection / simplification doesn't seem to deal well with reflections crossing the edge (see Cpt America poster in last photo, and the batman vs superman poster in the 4th image)
The canny edge detector will fail (and so does my brain) on posters when the frame/background difference isn't clear. I had to look twice myself at the second picture to make sure where the chrome frame ends and where the poster begins

So, this points to it being helpful to pre-process your images in different way, try to detect the rectangle, change your preprocessing, repeat, and so on.

A few observations:

maybe you can just reduce the aggressiveness of the contour simplification
Try different edge detectors! Canny is not always the best choice.
the weaker reflections just seem to be relatively white light – maybe convert your image to a brightness/hue/saturation colorspace, and reduce the variance of the brightness channel

As a side note: this approach (pyramid scaling, edge detection, contour extraction, contour simplification, rectangle detection) sounds more complex than just a normal Feature extraction by means of a transform of the image. For example, there are several transforms in OpenCV, some of which are very usable to detect dominant straight lines in the picture. If you can find parallelograms in those, you should be one step ahead, especially if you can really assume the poster is the biggest rectangle.

I think the Hough transform is quite exactly what you need.


Noisy Frequency Estimation using Hilbert Transform

Digital Hilbert transform is tricky to do since an ideal Hilbert transformer requires an FIR filter with an impulse response from -infinity to +infinity. Obviously that can't be done, so numerical problems are unavoidable. They tend to be worst at low frequencies so my guess is you'd see a smaller error at, say, 100 Hz. Andrew Duncan's paper "The Analytic Impulse" is thorough analysis of what the digital Hilbert transform can and can't do. See http://www.andrewduncan.ws/air/index.html

Some alternative methods of measuring frequency have been discussed here: Is there an algorithm for finding a frequency without DFT or FFT?


Trying to understand WAVELET TRANSFORM Frequncy-Time diagram

I believe by the WT, you are talking about the discrete wavelet transform, DWT.

This can be thought of as a subsampling of the continuous wavelet transform, CWT. In the case of the DWT, we pick frequencies of the form $2^{j-1}$ for ($j=1,2,\dots$) and then pick times seperated by multiples of $2^j$.

You can see this in the diagram: As frequency increases, boxes double in height and half in width.

I believe what the author means by the statement

Note that boxes have a certain non-zero area, which implies that the value of a particular
point in the time-frequency plane cannot be known.

is that the DWT, is an 'averaging' of the CWT within the boxes i.e. the boxes represent only one value and that is the average of the the values within it from the CWT. Hence as the boxes have non-zero area, it is impossible to know any of the exact values making up this average, so it is impossible to know the exact value of any point in the time-frequency plane.


What Is the Intuition of Convolution in The Signal Processing World

I think this post gives a lot of intuition about Convolution: - What's convolution intuitively?
In the Signal Processing world, an LTI (Linear and Time Invariant) system basically scales and sums delayed versions of the input signal.
This can be mathematically described using the integral known as convolution.


We always want to apply some kind of a window function in order to minimize the effect of leakage. This makes rectangular window (lack of any windowing) case never used, this is why:
Any tapering function used is almost always decreasing to zero at boundaries.

This is why we are losing some data. In order to retrieve that somehow you will usually do 50% of overlap when processing. This will retrieve whatever was in between.

Another thing is that if you apply the Inverse STFT, you should use complementary window, that is summing to 1, i.e. Hanning with 50%.

Finalising - yes, you should pretty much always use windowing in your applications.

For more comprehensive informations please refer to great white-paper:

  Heinzel G. - Spectrum and spectral density estimation by the DFT,
  including a comprehensive list of window functions and some new
  flat-top windows


Filtering noise and identifying phones boundary in speech recongition

system using arduino

Arduino performance is not enough for any serious speech recognition. You should better play with your mobile phone instead. With mobile you can build at least some interesting applications. Or try Raspberry Pi.

  From what I have read so far the entire process is the following

You probably need to read more papers (including the source code of existing software like CMUSphinx) to understand the details because the description of the process you have is not very accurate either. Ideally you need to read a textbook like "Fundamentals of speech recognition" by L. Rabiner.

  i. input signal

This is two lines of code.

  ii. filter all noise, getting all signals that is high energy volume and low frequency. (How do I know what to compare against? What is consider high energy volume, and low freqency?)

This is not very accurate description, there is no such thing as "low frequency". There is a voice activity detection (VAD) task which is pretty complex and has many dedicated algorithms. The simple form is VAD based on energy of the signal. You compare background energy level with current energy level and if the difference is higher than threshold you count audio as speech. This is not very accurate VAD but it is ok for many cases (no big noise, clean audio). It can be extended to frequency subbands and mixed with noise suppression by spectral subtraction. Noise suppression combined with VAD is even more complex task to solve.

For simple start you can implement energy-based VAD (google can give you many examples of implementation)

  iii. divide the signal into phones using FFT (where do I get dataset for phones, what constitutes a boundary between two phones, and how similar does the detected phones have to match the dataset?)

For DTW comparison with the prerecorded audio you do not need to care about phones at all. For HMM-based recognition, HMM decoding algorithm finds boundaries between phones implicitly but actually it doesn't even need them. Speech is considered as a sequence of subphone states, not as phones. Phones are not very stable to model with HMM. You can find description of the algorithm in the Rabiner's paper linked above.

The database of phone properties is extracted from transcribed recordings. There are many public databases available (some are mentioned on cmusphinx website).

  iv. matching the list of phones to each word using HMM or DTW

There is no such thing neither in HMM nor in DTW. You match the audio to the word model, not taking individual phones into account. You will have to implement a separate training algorithm.

  For beginners like me I feel that steps 2-3 is the hardest part due to the difficulties I listed.

Since your process is not very accurate the difficulty estimation is not very meaningful too. HMM training and decoding are quite complex tasks, VAD or FFT are simple.

The topic is not very easy, it is complex. You can not expect to implement calculus algorithm without knowing the theory of calculus, like what is integration and so on, so the speech recognition. There are many tutorials around, this one is ok to get started too:

http://kastnerkyle.github.io/blog/2014/05/22/single-speaker-speech-recognition/

Maybe it will be helpful for you. However, take into account that any simple tutorial is usually not very professional and doesn't really cover things. Speech recognition can not really be explained in 10 pages of text. If you really want to understand a topic, read a textbook or attend a course on ASR.


If your application can use LGPL open-source code, you could reuse the code for the Remez exchange algorithm in the GNU Octave source.  I've got a standalone library that does that at www.ke5fx.com/dsplib.zip .  See if the function DSP_FIR_PM_Hilbert_transformer_design() does what you want.

The FIR filters implemented by my library are one-dimensional affairs, so if you actually need to do a 2D transform, that would be left as an exercise, so to speak.  Contributions back to dsplib are welcome!


Nyquist theorem vs sampling theorem vs shannon sampling theorem?

Sampling operation has its roots in the mathematical interpolation theory which was used to generate certain function values at specified points from the availabe set (the samples) of existing values. These kind of work is summarized as Whittaker interpolation. Lagrange interpolation is also another related concept.

Sampling theorem in electrical engineering, however, has its roots in the communication engineering literature and pioneered by communications engineer Harry Nyquist, who was among the first to consider digital transmission of analog signals. The Nyquist sampling theorem states the minimum number of uniformly taken samples to exactly represent a given bandlimited continuous-time signal so that it (the signal) can be transmitted using digital means and reconstructed (exactly) at the receiving end.

Later, C.E.Shannon also worked on communications theory (especially the digital one) more from a probabilitistic point of view and published his conclusions as Information Theory within his classic paper of The Mathematical Theory of Communications. Among the works shannon did also include a more in depth analysis of sampling operation and resulting Shannon Sampling theorem which is identical to Nyquist sampling theorem for bandlimited signals.

Hence the uniform sampling of bandlimited continuous-time signals is referred to as Nyquist-Shannon sampling theorem. You may omit the inventors' names and simply call it sampling theorem remembering that it's the uniform sampling of bandlimited signals.


The resulting FFT will be different: instead of calculating the FFT at the frequencies $2\pi n/1800\ \ $ for $n=0\ldots 1799\ \ $, you will be calculating them at $2\pi n/2048\ \ $ for $n=0\ldots 2047\ \ $. However, there is no degradation of the information.

Both approaches are correct: using 1800 or 2048. I'd use the "minimum energy solution" (i.e. the simplest, the laziest solution).  This would usually be using the 2048 length transform.

People tend to use radix-2 transforms because they don't know any better. There seems to be much mis-information about FFTs HAVING to be power-of-2. There is no such constraint.  Also, they probably don't know about about decent non-radix-2 algorithms, such as those available in FFTW and other libraries.

To see that the FFT of any length is information-preserving:

  Assume your data of length 1800 is in x.  Form X = fft(x,2048); (or some other  length different from 1800).
  Find xx = ifft(X);.
  See what sum(abs(x-xx(1:1800))); is.

See also this question and its answers.


necessity of patch size in edge detection algorithms

Let's start to given some information.

We now that the sobel operator perform a aproximmation of the image gradient using convolution ( $ * $ ).

Then,

Being $\mathbf{A}$ an image of $n \times m$ size.

The operator $\mathbf{K}$, commonly, uses two $3 \times 3$ kernels which are convolved with the original image. Thus we have

$\mathbf{G_x} = \mathbf{K_x} * \mathbf{A} $

$\mathbf{G_y} = \mathbf{K_y} * \mathbf{A} $

Being
$$
  \mathbf{K_x} =
  \left[ {\begin{array}{cc}
                -1 & 0 & +1\\
                -2 & 0 & +2\\
                -1 & 0 & +1\\
              \end{array} } \right]
$$

And,

$$
  \mathbf{K_y} =
  \left[ {\begin{array}{cc}
                -1 & -2 & -1\\
                 0 & 0 & 0\\
                +1 & +2 & +1\\
              \end{array} } \right]
$$

So, an extension of this kernel is a $5 \times 5$ matrix

$$
  \mathbf{K_y} =
  \left[ {\begin{array}{cc}
                -1 & -4 & -6  & -4 & -1\\
                -2 & -8 & -12 & -8 & -2\\
                 0 &  0 &  0  &  0 &  0\\
                +2 & +8 & +12 &  8 &  2\\
                +1 & +4 & +6  & +4 & +1\\
              \end{array} } \right]
$$

$$
  \mathbf{K_x} =
  \left[ {\begin{array}{cc}
                 +1 &  +2  & 0 & -2  & -1\\
                 +4 &  +8  & 0 & -8  & -4\\
                 +6 &  +12 & 0 & -12 & -6\\
                 +4 &  +8  & 0 & -8  & -4\\
                 +1 &  +2  & 0 & -2  & -1\\
              \end{array} } \right]
$$

This article demonstrates how to get Sobel gradient operators analytically. For that was using linear approximation of brightness in window 3x3.  Every pixel in window has its own weight. Special weight matrix selection gives Sobel gradient operator.


Time invariance from convolution integral

If your delayed input is $\hat{x}(t)=x(t-\alpha)$, then $\hat{x}(t-\tau)=x(t-\tau-\alpha)$, and not $x(t-\tau+\alpha)$. From here everything will work out as expected.


Is it possible to do further beamforming on already beamformed data?

For now I will focus on the conceptual, pending feedback from you:

  What I'm interested in to know if it is possible to take two adjacent
  beams, and do another step of beamforming,

The answer to this is yes, this is possible to do, and this is done in many beamforming applications. Think about it like this: All beamforming is really doing, is undoing the delays your signals exhibit (by virtue of different sensors at different locations), and summing the result. (Example: Delay-And-Sum beamformer). But in what order you undo the delays, before summing doesn't matter.

Thus, the beamformed result of a sub-group of sensors can be considered as one 'element', and the beamformed result of another sub-group of sensors, another 'element', and then both 'element's' outputs are beamformed.

In other words, first we remove the effects of the delays in one sub-array. Then we remove the effects of the delays in the other sub-array. At this point, both signal outputs might still have a delay relative to each other. So then again, both those inputs can be beamformed again. You can keep doing this ad infinitum.

  The objective is to reduce my side lobes even further, and to get
  narrower main beam.

If this is the objective, then the only way to reduce your main lobe width is to sample the spatial field with more sensors. This is the direct analogue to time-domain sampling. If we want greater frequency resolution (smaller main lobe width in the DFT analysis), we need to process a greater number of samples in time. Similarly, if we would like a smaller main-lobe width of our spatial spectra, we need to sample the spatial field with more sensors.

EDIT: Taking from the above, if you want to decrease your mainlobe width, you need more spatial samples. However if you want to reduce your sidelobes, then you would use a variety of windowing/shading functons. (Each of those functions will give you a slight increase in mainlobe width, but that is the price to pay). I am not sure what shading you have used for the above plots, but that will be good to know.


Can the intelligence message of an FM signal be cross-modulated to another center frequency due to non-linearities?

Yes and no,   it can't be intermodded onto another SIGNAL, but it can be intermodded onto another CHANNEL.  For example, if you are monitoring 144 MHz, an FM signal on 145 MHz and on 146 MHz CAN, due to non-linearity, end up being heard on 144 MHz.  But it can't be modulated onto an existing signal on 144 MHz.
Mark


What are the most common algorithms for adaptive thresholding?

I do not think mine will be a complete answer, but I'll offer what I know and since this is a community edited site, I hope somebody will give a complimentary answer soon :)

Adaptive thresholding methods are those that do not use the same threshold throughout the whole image.

But, for some simpler usages, it is sometimes enough to just pick a threshold with a method smarter than the most simple iterative method. Otsu's method is a popular thresholding method that assumes the image contains two classes of pixels - foreground and background, and has a bi-modal histogram. It then attempts to minimize their combined spread (intra-class variance).

The simplest algorithms that can be considered truly adaptive thresholding methods would be the ones that split the image into a grid of cells and then apply a simple thresholding method (e.g. iterative or Otsu's method) on each cell treating it as a separate image (and presuming a bi-modal histogram). If a sub-image can not be thresholded good the threshold from one of the neighboring cells can be used.

Alternative approach to finding the local threshold is to statistically examine the intensity values of the local neighborhood of each pixel. The threshold is different for each pixel and calculated from it's local neighborhood (a median, average, and other choices are possible). There is an implementation of this kind of methods included in OpenCV library in the cv::adaptiveThresholding function.

I found another similar method called Bradley Local Thresholding. It also examines the neighborhood of each pixel, setting the brightnes to black if the pixels brightness is t percent lower than the average brightness of surrounding pixels. The corresponding paper can be found here.

This stackoverflow answer mentiones a local (adaptive) thresholding method called Niblack but I have not heard of it before.

Lastly, there is a method I have used in one of my previous smaller projects, called Image Thresholding by Variational Minimax Optimization. It is an iterative method, based on optimizing an energy function that is a nonlinear combination of two components. One component aims to calculate the threshold based on the position of strongest intensity changes in the image. The other component aims to smooth the threshold at the (object)border areas. It has proven fairly good on images of analog instruments (various shading and reflection from glass/plastic present), but required a careful choice of the number of iterations.

Late edit:
Inspired by the comment to this answer. There is one more way I know of to work around uneven lighting conditions. I will write here about bright objects on a dark background, but the same reasoning can be applied if the situation is reverse. Threshold the white top-hat transform of the image with a constant threshold instead of the original image. A white top hat of an image is nothing but a difference between the image $f$ and it's opening $\gamma(f)$. As further explanation let me offer a quote from P. Soille: Morphological Image Analysis:

  An opening of the original image with a large square SE removes all relevant image structures but preserves the illumination function. The white top-hat of the original image or subtraction of the illumination function from the original image outputs an image with a homogeneous illumination.


How to set borders for image translation

What you observe in the image is the effect of the circular shift implied by the multiplication of the DFT of the image by the complex exponentials; i.e., the last (right most) sample becomes the first (left most) sample, hence circularly shift around the image, after a single pixel shift to the image is performed.

Theoretically you cannot prevent the shift phenomenon if you use the same method to translate the image; i.e., by using the inverse DFT of the exponentially multiplied forward image DFT.

However you can practically make it unobservable if you pad the original image by enough number of zeros at the ends so that the shifts of zeros will be unnoticed and only the linear translation of the image will be apparent.

To do so requires the amount of shift to be known in advance. Then given the shift amount and direction, you must pad that number of zeros in the correct end of the image so that when those zeros are circulated to the opposite end, they will not create a noticeable effect and the image will seem to be linearly shifted all.


How to Classify a Kernel as Low Pass Filter (LPF) or High Pass Filter (HPF)? How to Transform an LPF Kernel into HPF Kernel?

General

We assume 2 modes of filters: LPF or HPF.

Classifying Filter Type

Usually, if it is a well planned LPF and well Planned HPF a simple test will do.
Calculate the sum of all coefficients.
The sum of the coefficients is the first element of the DFT of the signal.
It means it is the DC gain and well behaved LPF has gain of 1 and well behaved HPF have DC gain of 0.

So sum the coefficients and:

If the sum is $ 1 $ then this is an LPF.
If the sum is $ 0 $ then this is an HPF.

Pay attention that those method can't differentiate between HPF and Band Pass or between LPF and Band Stop.

The more general method is to look on its DFT Spectrum to see its behavior.

Converting Filter Type

There is very intuitive formula - ALL_PASS_FILTER = LPF_FILTER + HPF_FILTER.
It's not accurate but this is the idea.

Be more accurate in the Frequency Domain All Pass Filter is a filter which has value of 1 for any frequency.
Assuming we have LPF Filter then HPF_FILTER = ALL_PASS_FILTER - LPF_FILTER.
In the spatial domain it means $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] $.

Where $ \delta \left[ n \right] $ is the filter with 1 in its origin and $ 0 $ anywhere else.

In your case $ \delta \left[ n \right] = [0, 0, 0, 0, 1, 0, 0, 0, 0] $.
So the HPF Filter becomes $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] = [-1/9, -1/9, -1/9, -1/9, 8/9, -1/9, -1/9, -1/9, -1/9] $.


Below is a block-diagram of your signal processing chain (ignoring multiplication by 2).

$$ x[n] \longrightarrow \boxed{ \text{ DFT-N } } \longrightarrow X[k] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow \boxed{ \uparrow 2 } \longrightarrow Y[k] \longrightarrow \boxed{ \text{ IDFT-N } } \longrightarrow y[n]$$

Where sequence $x[n]$ has a length of $N=16$ samples, $X[k]$ is its $N$-point DFT sequence. And $Y[k]$ is related to $X[k]$ as given by
$$ Y[k] = \begin{cases} { 2 X[k] ~~~ , ~~~ k = 0,2,...,N-2 \\ 0 ~~~ , ~~~ \text{ otherwise } \\} \end{cases} $$

I will solve this problem using the duality argument. Eventhough it's not necessary, it can reduce the mathematical burden via bypassing explicit computations.

The duality in DFT stems from the fact that the forward and inverse DFT transforms are (almost) identical, designated with the following:  ($n$ time , $k$ frequency)

$$ x_t[n] \leftrightarrow X_f[k] \iff X_f[n] \leftrightarrow N x_t[-k] $$

We can use the duality property to observe that if an operation  $T$ is performed on the time-domain signals $x[n] \xrightarrow{T} y[n]$ with a corresponding reflection on the frequency domain as $X[k] \xrightarrow{M} Y[k]$, where $M$ is the frequency-domain reflection operator of the time-domain operator $T$ , then the same operation $T$ performed on the freqency domain signals will have the same reflection operator $M$ on the corresponding time domain signals, with a possible scale factor if different sample rates are used on the left and right nodes.

Therefore, to use duality, in the above block diagram, I'll replace the frequency domain relation between $Y[k]$ and $X[k]$ with a time domain relation as:

$$v[n] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow w[n] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow z[n] $$

where in the time domain the operation is compression (downsampling) and expansion (upsampling) whose frequency domain DFT relation is readily available as:

$$ W[k] = 0.5 ( V[k] + V[k+N/2] ) ~~~, k = 0,1,..,N/2-1 $$

where $W[k]$ is $N/2$ point DFT of $N/2$ point compressed sequence $w[n]$, and $V[k]$ is $N$ point DFT of $N$ point sequence $v[n]$. Stated verbally, $W[k]$ is obtained from $V[n]$ by adding its second half onto its first half.

Then
$$Z[k] = W[ (k)_{N/2}] ~~~, k = 0,1,2,...,N-1 $$

Where $Z[k]$ is the $N$ point DFT of length $N$ point sequence $z[n]$. Again stated verbally, $Z[k]$ is merely a repetation of $W[k]$ two times. This is signified by the ramge of $k=0,1,2,...,N-1$ and the modulus operator inside the argument of $W[k]$.

Hence we can find the relation between $Z[k]$ and $V[k]$ as :
$$Z[k] = 0.5 \left( V[ (k)_{N/2}] + V[ (k)_{N/2} + N/2] \right) $$

This relation in the frequency domain was readily available. I didn't use any derivations to achieve it (eventhough at first it was derived by means of DFT manipulations anyway). Now using the duality property of DFT, I assert that, if te time domain signals are replaced by frequency domain signals then their effect on the time domain will similar to the above frequency domain relation, by simply replacing $k$ with $n$. Note that time reversal not needed, as it will be cancelled, but scaling by $2$ is required. So the result will be:

$$ \boxed{ y[n] =   x[ (n)_{N/2}] + x[ (n)_{N/2} + N/2] } $$

where $N=16$, $n=0,1,2,...,15$.

Conclusion. It takes more effort to describe a solution more easily obtained via the duality principle.


difference between conventional analog butterworth filter and that implemented through IIR

If you have designed the IIR Butterworth filter using the bilinear transform method, then due to the logarithmic frequency mapping (warping), the frequency and phase responses will be distorted...


Question regarding filter implementation - Audio EQ Cookbook

With the values you provided ($Q=1$, $f_0=5355\,\text{Hz}$, and $F_s=48\,\text{kHz}$) I get the following filter coefficients from the Cookbook formulas:

b =

   0.11789
   0.23578
   0.11789

a =

   1.32248
  -1.52844
   0.67752

Note that the formulas give you denominator coefficients that are not normalized. Many implementations use $a_0=1$, which can be easily achieved by dividing numerator and denominator coefficients by $a_0=1.32248$.

The maximum pole radius of this filter is $0.72$, which means that we're not at all close to the unit circle, so stability shouldn't be an issue, no matter which implementation you use.

You can add more information if you like, but from what I know up to now, the most likely problem is a bug in your code.

Just test your filtering routine by using as input sequence x=1,0,0,...,, which should give the following output (the impulse response):

ans =

   8.9144e-02
   2.8131e-01
   3.6860e-01
   2.8188e-01
   1.3694e-01
   1.3860e-02
  -5.4140e-02
  -6.9672e-02
  -5.2786e-02
  -2.5313e-02
  -2.2122e-03
   1.0411e-02
   1.3166e-02
   9.8827e-03
   4.6767e-03
   3.4194e-04
  -2.0007e-03
  -2.4875e-03
  -1.8499e-03
  -8.6361e-04
  -5.0389e-05


Does the FFT size matter in OFDM demodulation?

What is the impact of taking an arbitrary larger FFT size during the demodulation process?

It doesn't work.
This gets pretty intuitive if you take an arbitrary OFDM introduction and look at the spectrum of OFDM, consisting of overlapping sinc functions: you need to hit exactly that raster, or you don't see the individual channels, but an undisambiguable mishmash of carriers.
Using a larger FFT at the same sample rate can't work: don't forget that the FFT is just a DFT. A DFT maps vectors of a specific length N from time domain to frequency domain vectors of exactly the same size. So, a larger FFT length simply needs more time domain samples. But your OFDM symbol is as as many samples long as the transmit FFT was.
So, no, you can't use any other FFT size, sensibly.
You'd have to increase the sample rate by the same integer factor that you've increased the FFT size.
But that would only yield an FFT where all but the original bins are only noise, so you'd have won nothing.
(This is actually done in some OFDM systems to have more headroom for frequency synchronization, but it's not desirable if you don't have a frequency offset larger than your outer guard carriers.)


Response of a system to a step function (heaviside)

Fat32's derivation of the result via the Fourier transform is correct, but I think that your original question hasn't really been answered ("what am I doing wrong?"). The real reason why WolframAlpha gives you a wrong result is the different definitions of the Fourier transform used by WolframAlpha and the one used by you to derive $H(\omega)$. By default, WolframAlpha uses this definition of the Fourier transform:

$$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{j\omega t}dt\tag{1}$$

Apart from the factor $1/\sqrt{2\pi}$ (which makes the transform unitary), the actual problem lies in the sign of the exponent. You implicitly assume a negative sign, whereas WolramAlpha uses a positive sign.

There are two ways to solve the problem. Either you change Wolfram Alpha's default parameters using FourierParameters -> {1, -1} as explained here, or you change $H(\omega)$ accordingly, which means that you would need to use

$$H(\omega)=\frac{1}{1-j\omega RC}\tag{2}$$

As an example I tried the second method, and WolframAlpha indeed provides the correct result.

One more thing, for this kind of problems it's usually more convenient to use the Laplace transform because this avoids the use of generalized functions, such as the Dirac delta.


(Note: the paper pointed by hotpaw2's link is actually describing in more detail the algorthm I presented here)

Consider a data window length of $N$ samples from $n=0$ to $n=N-1$. Let your original data window be $x_1[n]$, whose first sample is $x_{old} = x_1[0]$. Now your new data set is denoted as $x_2[n]$ whose samples are actually one sample left shifted verisons of $x_1[n]$, i.e. $x_2[n] = x_1[n+1]$ for $n=0,1,...,N-2$ plus a new arrived data to position $n=N-1$ which is denoted as $x_{new} = x_2[N-1]$.

Then the following algorithm will compute the N-point DFT, $X_2[k]$ of the new data set $x_2[n]$ from that of the already computed and stored N-point DFT $X_1[k]$ of the old data set $x_1[n]$ as:

$$ X_2[k] = e^{j \frac{2\pi}{N}k }( X_1[k] + (x_{new} - x_{old}) ) $$

for each $k=0,1,...,N-1$.  This updated computation of $X_2[k]$ from pre-computed $X_1[k]$ requires N complex multiplications and N real additions. Compared to a direct N-point FFT which requires $N \log_2(N)$ complex MACs, this would therefore be an improvement by a factor of $\log_2(N)$ which for example would render to a factor of 10 when $N=1024$ using low level languages such as C.


Both Impulse and Step Response can help identify a system, however it might be hard to actually perform the identification. You'd have to fourier transform the response. Basically, this is only possible using digital signal processing. Using a sine sweep, the result is easily read, but the sweep will take a relatively long time.


Proof that first difference filter amplifies noise

Any filter whose frequency response magnitude is above unity at a given frequency $\omega$, amplifies signals at that frequency.

In discrete-time, the typical white noise is a wideband
signal whose power spectrum ranges from $-\pi$ to $\pi$.

The digital  first difference filter is an LTI system with the impulse response $h[n] = \delta[n] - \delta[n-1]$ and frequency response $H(\omega) = 1 - e^{-j\omega}$, whose magnitude is:

$$ |H(\omega)| = \sqrt{ 2 - 2 \cos(\omega) } = \sqrt{ 4 \sin^2(\omega/2) } = 2|\sin(\omega/2)|$$

as you can see as $\omega \to \pi$, then the magnitude goes to $2$ which is larger than unity and therefore amplifies any noise content found at those frequencies.

Note that, in any digital system, there will always be noise due to inherent quantization of numerical representation format, and therefore it will always be amplified using the given first difference filter.


Ideal high pass filter for discrete signal

HINT: Your original LPF is $-\omega_c \lt \omega \lt +\omega_c$. When you shift it by $\pi$, your integral limit for second method is $\int_{-\pi}^{-\pi + \omega_c}$ and $\int_{\pi - \omega_c}^{\pi}$.

For the first method, it should be
$$
h[n] = (-1)^{(n-\alpha)}\frac{sin(\omega_c(n-\alpha))}{\pi(n-\alpha)}
$$
because you have to apply shift of $\alpha$ on $n$. This will fix your plot of (1).

For (2), with the changed limits as I hinted at the beginning
$$
h[n] = \frac{1}{2\pi}\int_{-\pi}^{-\pi + \omega_c}e^{-j\alpha\omega}e^{j\omega n}d\omega + \frac{1}{2\pi}\int_{\pi-\omega_c}^{\pi}e^{-j\alpha\omega}e^{j\omega n}d\omega\\
= \frac{1}{2\pi j(n-\alpha)}[e^{-j(n-\alpha)(-\pi+\omega_c)} - e^{j(n-\alpha)(-\pi)}] +\frac{1}{2\pi j(n-\alpha)}[e^{-j(n-\alpha)(\pi)} - e^{-j(n-\alpha)(\pi-\omega_c)}] \\
= \frac{1}{2\pi j(n-\alpha)}e^{-j(n-\alpha)(-\pi)}[e^{-j(n-\alpha)\omega_c} - 1]+\frac{1}{2\pi j(n-\alpha)}e^{-j(n-\alpha)(\pi)}[1 - e^{j(n-\alpha)\omega_c}]\\
= \frac{1}{2\pi j(n-\alpha)}(-1)^{(n-\alpha)}[e^{-j(n-\alpha)\omega_c} - e^{j(n-\alpha)\omega_c}]\\
= (-1)^{(n-\alpha)}\frac{sin(\omega_c(n-\alpha)}{\pi(n-\alpha)}
$$


Running the BM3D Denoising Algorithm in Python

What I resorted to was using the PyPI package, which is advertised here: http://www.cs.tut.fi/~foi/GCF-BM3D/index.html#ref_software.

I dug a bit in the source code, and found that I could perform BM3D, in the following fashion:

import bm3d

denoised_image = bm3d.bm3d(image_noisy, sigma_psd=30/255, stage_arg=bm3d.BM3DStages.HARD_THRESHOLDING)

There are also some examples in the library's source code download.

I installed bm3d using pip (pip install bm3d) and needed OpenBlas (sudo apt-get install libopenblas-dev).


Downsampling lowpass filter for audio: FIR or IIR?

If linear phase is a requirement, that will probably steer you toward an FIR implementation. It is possible to build IIR filters that have approximate linear phase, but it is easy to design a linear-phase FIR.

If you're concerned about latency, forward-backward filtering as in filtfilt isn't really a good option. In general, it's really meant to be used an offline process, since to implement the technique exactly, you have to run the entire signal through forward, then do the same in reverse. That implies that you have access to the entire signal at once, which is not commensurate with low delay.

In general, an FIR filter will require a higher order for a given set of performance requierments. However, FIR filters bring some real advantages, such as guaranteed stability, lower susceptibility to roundoff errors (since the quantization error doesn't get fed back through the filter, although you can compensate for this with some increased complexity), and simply-achieved linear phase response. In addition, efficient FIR filter implementations are available for many processor architectures, mitigating the cost of the extra taps somewhat.

Another way to mitigate the extra-tap cost for FIR filters in your situation is to take advantage of efficient multirate signal processing techniques. Specifically, you can use a polyphase decimation approach to significantly reduce the number of computations that you perform in the process of decimating your signal. This has the effect of decreasing the number of effective taps (in terms of computational complexity) in the decimation filter. In addition, if you need to decimate by a large factor, then multistage approaches can help further reduce your load. Lyons' introductory DSP book has some good easy-to-read material on these topics.

Given more specific parameters of your system, one could make more pointed recommendations. What are your filter design requirements? What sort of computational capabilities does your platform have? What sample rate will you decimate to?


The result is expected as you are simply appending zeros to the product

fft(A1).*fft(A2)

Which does the time domain interpolation as you describe. However the ifft of the above product results in a series of constant values in the time domain since it is a circular convolution. Appending additional zeroes simply interpolates this for more constant values (in addition to a scaling factor).

To achieve the result you desire, which is the linear correlation of the two vectors, you must interpolate in the frequency domain prior to the multiplication (not insert zeros). This is exactly what occurs if you review the result of fft(A,M) which appends zeros in the time domain prior to taking the FFT, resulting in frequency domain interpolation.

Essentially, the process is to zero-pad the time domain samples prior to circular convolution, which results in linear convolution. Zero padding in frequency is not actually introduced in the process ifft(fft(A,M).*fft(B,M)).


To implement successive interference cancellation (SIC) process, should the receiver know the power allocated to each symbol?

Yes, one would require to know the power of other users. There are two options, the network can either signal the UEs this information in a control frame or the UEs need to estimate the power of the streams allocated for other users (knowing its own power allocated and some smart signal processing).
In any case interference cancellation techniques are suspetible to the "effectiveness" of the cancellation. Usually the SNR degardes rapdily with decrease in estimation accuracy or error in knowledge of the parameters of the other users, such as power


Loop Bandwidth for Carrier Demodulation loop

Your results are not showing that the data is being "filtered out" necessarily but that you are using a non-linear phase filter resulting in group delay variation over your passband, as evidenced by the spread in phase.

Don't use a Butterworth filter to model this. Instead use filters derived from either firls or firpm (remez in Octave) which are linear phase and have constant group delay.


Modelling propagation of sound wave by particle simulation

There is a huge body of research and commercial work on doing acoustic simulation. There are a two basic classes of simulation systems:

Finite element or finite boundary: these model the actual physics by dividing air and or/surfaces into small patches and modelling step by step the interaction between the patches by locally solving the wave equations given the boundary and initial conditions of the patches near by
Approximation methods: Ray tracing, particle tracing, mirror image method, etc. These are macroscopic models that use certain simplifying assumptions such as specular reflections, simplified diffusion models, no diffraction, etc.

This is a very complicated problem and it takes many person years to create a modelling system that actually produces useful results. There are commercial systems out there that you can buy. They vary greatly in complexity, capability, cost, and training required so it all depends on your specific application.


When I ran your code in Matlab I got the following warning message:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  6.912349e-21.
> In firls at 175

It turns out that the problem is the "/Fs/2" part.  You are trying to normalize the frequencies from $0$ to $\frac{F_s}{2}$ into the range $0$ to $1$.  The problem is that Matlab divides the numbers in order, so even though what you want is /(Fs/2) what you are getting is /(Fs*2).  Changing the line to the following-

b = firls(N,[0 fpass fstop Fs/2]/(Fs/2), [1 1 0 0], [1 1]);

fixes the problem.


It's due to the weighting function used in the Gaussian pyramid downsampling.

Lets take a 1D patch(1D patch because it makes our calculation easy and also the kernel is separable) with p10,p11,p12,p13,p14,p15 in the Level 1.

Assume maximum gradient possible is Max - Min between any two pixels.

The pixel p20 in level 2 is formed by influence of pixels p10,p11,p12,p13,p14 in level1 and pixel p21 if formed by influence of pixels p12,p13,p14,p15,p16. The corresponding image(I'm not not able to upload image).

The weights of the kernel  are

w(0) = a,

w(-1) = w(1) = 1/4

w(-2) = w(2) = 1/4 -a/2.

  For Gaussian pyramid a = 0.4(approximately)

  w(0) = 0.4

  w(-1) = w(1) = 0.25

  w(-2) = w(2) = 0.05

The equation for new pixel p20

  p21 = w(-2)*p10 +w(-1)*p11+w(0)*p12 +w(1)*p13+w(2)*p14

Substituting the values for the weights we get,

  p20 = 0.05*p10 + 0.25*p11 + 0.4*p12 + 0.25*p13 + 0.05*p14

Similarly value of pixel p21

  p21 = 0.05*p12 + 0.25*p13 +0.4*p14 +0.25*p15 + 0.05*p16

Now the difference between pixels p20 and p21 is,

  p20 - p21 = 0.05*p10 + 0.25*p11 +0.35*p12 -0.35*p14 -0.25*p15+0.05*p16

So for the difference to be maximum p10 = p11 =p12 = Max and p14=p15=p16 = Min

  p20 -p21 = 0.7(Max-Min)

  p20- p21 = (Max -Min)/√2

Since the filters are separable, applying the kernel in another direction over the image again will reduce the pixel difference by a ratio of √2.
The maximum gradient in x and y direction in next level of Gaussian pyramid is (Max - Min)/2

The maximum gradient possible is reduced by half or in frequency domain reduces the frequency by half(an octave).


Pre-emphasis: What is connection between parameter a and 6db/oct?

You're right that the value $\alpha$ will not influence the slope of the filter (approx. $6$ dB/octave). That slope is determined by the filter order, which is $1$ in this case (1 delay element in the filter implementation). The value $\alpha$ only influences the filter's cut-off frequency.


How to prove that the peak of the autocorrelation function is at zero lag?

The Cauchy Schwarz inequality states that:
$$
\left|\int_{-\infty}^{\infty}g_1(t)g_2(t) dt\right|^2 \leq \int_{-\infty}^{\infty}|g_1(t)|^2 dt \int_{-\infty}^{\infty}|g_2(t)|^2 dt
$$

I'm going to assume that $f(t)$ is real, just to make the math a little easier. From the above we can write:
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq \int_{-\infty}^{\infty}|f(t)|^2 dt \int_{-\infty}^{\infty}|f(t-\tau)|^2 dt
$$
For the second integral on the right, just use a variable substitution - let $u=t-\tau$ and then rewrite it with the dummy variable $t$ again - we have
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq \int_{-\infty}^{\infty}|f(t)|^2 dt \int_{-\infty}^{\infty}|f(t)|^2 dt
$$
but the integral on the right is just the autocorrelation function i.e. $ \int_{-\infty}^{\infty}|f(t)|^2 dt= R_{f}(0)$. So we can write
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq R^2_{f}(0)
$$
Note that from the definition of $R_f(0)$ it must be a positive quantity, so taking the square root of both sides and paying attention to the positive and negative quantities on the left hand side, we can write
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau) dt \leq R_{f}(0)
$$

Note - you don't have to use the Cauchy Schwarz inequality - another proof:
$$[f(t)-f(t-\tau)]^2 = f^2(t)+f^2(t-\tau) -2f(t)f(t-\tau)
$$
Integrating both sides from $-\infty$ to $\infty$ and using the same trick (variable substitution for the integral over $f^2(t-\tau)$) we can write
$$
\int_{-\infty}^{\infty} [f(t)-f(t-\tau)]^2dt= 2\int_{-\infty}^{\infty}f^2(t)dt-2\int_{-\infty}^{\infty}f(t)f(t-\tau)dt
$$

Rearranging we can write:
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau)dt = \int_{-\infty}^{\infty}f^2(t)dt -\frac{1}{2}\int_{-\infty}^{\infty} [f(t)-f(t-\tau)]^2dt
$$
Now - the second integral on the right has to be a positive quantity and it is subtracted so we have
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau)dt \leq \int_{-\infty}^{\infty}f^2(t)dt
$$


Help in analytical expression of BER for M-QAM

You need to read a good book to form a solid foundation first and then read research papers to understand advanced concepts. Eq (12) comes through a standard procedure for BER calculations. See Digital Communications by John Proakis, for example.


Possibilties of chebyshev polynomial waveshaping

For $-1 <= x <= 1$, let's compare Chebyshev polynomials of the first kind, $T_n(x)$, and the basis functions of the Fourier cosine series, $F_n(x)$:

 $F_n(x)=\cos(n \pi x)$

 $T_n(x)=\cos(n\ \text{acos}\ x)$

Writing $T_n(x_T) = F_n(x_F)$ and solving for $x_F$ gives $x_F = (\text{acos}\ x_T) / pi$, revealing that the Chebyshev polynomial series is simply an argument-warped version of the Fourier cosine series that maps the range $0 <= x_F <= 1$ to $-1 <= x_T <= 1$:

This exclusion of $x_F < 0$ makes irrelevant the constraint that the Fourier cosine series is only applicable to even functions, functions symmetric around $x_F = 0$, meaning that $f(x_F) = f(-x_F)$. The mapping function is well-behaved in the range of interest, so if the Fourier cosine series can represent any even function, then the Chebyshev polynomial series can represent any function.

Not knowing the amplitude of the input sinusoid, it is not possible to control the proportional amplitudes of the generated harmonics by the choice of the waveshaper. For example if the waveshaper is $f(x) = x^2 + x$:

$$f(\cos x) = \frac{1}{2} + \cos x + \frac{1}{2} \cos 2x\\
f(\frac{1}{2} \cos x) = \frac{1}{8} + \frac{1}{2} \cos x + \frac{1}{8} \cos 2x$$

The first one has the amplitude of the second harmonic at -6 dB compared to the fundamental, and the second one at -12 dB. Chebyshev polynomials are no miracle cure. The above polynomial $x^2 + x$ was a weighted sum of the first three Chebyshev polynomials. Any polynomial can be written as a weighted sum of Chebyshev polynomials, and vice versa.

The problem persists for pure Chebyshev polynomials, for example $f(x) = 4 x^3 - 3x$:

$$f(\cos x) = \cos 3x\\
f(\frac{1}{2} \cos x) = \frac{1}{8} \cos 3x - \frac{9}{8} \cos x$$

A Chebyshev polynomial series does have the advantage of better numerical stability than a simple power series, and if the amplitude of the input sinusoid is controlled somehow, it can be used to control the level of harmonics in the way you describe.


Why is the correlation of two random markov chains so large?

Hi: you can't cumulate two sets of random numbers and then calculate a correlation between them. The two samples are non-stationary ( since you are cumulating them )  and what you get is referred to as "spurious correlation" or "non-sensical correlation" in econometrics. ( some would say statistics but the econometricians really figured out what was happening ).  I wouldn't do it justice trying to explain it here but, essentially, what you are doing is calculating two series whose means are changing over time so the correlation calculation is illegal statistically. Granger and Newbold figured out the statistics behind why this is in 1976 so if you want all the details, just google for "spurious regression". The econometrics literature on this is enormous and, even though what you are doing has maybe zero to do with econometrics, you are basically generating two series that are each random walks and then correlating them which is a garbage in garbage type of scenario.  I hope this helps.


Compute step response from impulse response of continuous-time LTI system

Since your question is not about how to obtain the answer, but specifically about the reason why the unit step $u(t)$ appears in it, I assume that you know that the response to a step input is computed by convolving the impulse response with a step:

$$y(t)=\int_{-\infty}^{\infty}h(\tau)u(t-\tau)d\tau=\int_{-\infty}^{t}h(\tau)d\tau\tag{1}$$

Since $h(t)$ is causal, i.e., $h(t)=0$ for $t<0$, the lower integration limit can be changed to $0$. This implies that if in $(1)$ the upper integration limit $t$ is less than $0$, the result must be zero. Consequently, the response to a step input can be written as

$$y(t)=u(t)\cdot\int_0^{t}h(\tau)d\tau\tag{2}$$

Evaluating $(2)$ gives you the final answer, which includes the step function $u(t)$.

EDIT: Triggered by Dilip Sarwate's comment I add some extra explanation. Note that the step function in $(2)$ is actually redundant because for $t<0$ the integral in $(2)$ is zero because for $t<0$ we have

$$\int_0^{t}h(\tau)d\tau=-\int_{-|t|}^0h(\tau)d\tau=0,\quad t<0$$

because $h(t)=0$ for $t<0$. However, if $h(t)=f(t)u(t)$ with some function $f(t)$ that is not zero for $t<0$ (in the given example we have $f(t)=e^{-4t}$), we cannot just replace the lower integration limit in $(1)$ by $0$ and leave out the step function:

$$\int_{-\infty}^th(\tau)d\tau=\int_{-\infty}^tf(\tau)u(\tau)d\tau=\int_{0}^tf(\tau)u(\tau)d\tau\neq \int_{0}^tf(\tau)d\tau$$

because now we generally have for $t<0$

$$\int_{0}^tf(\tau)d\tau=-\int_{-|t|}^0f(\tau)d\tau\neq 0$$

The correct way to evaluate the integral for $h(t)=f(t)u(t)$ is

$$\int_{-\infty}^th(\tau)d\tau=u(t)\cdot\int_0^tf(\tau)d\tau$$


Your sequence simply is a sampled version of the continuous signal $x(t)$. So

$$x_d(n)=x(nT),\quad T=\frac{1}{50}$$

The relation between the CTFT of $x(t)$ and the DTFT of $x_d(n)$ is that the DTFT is the sum of shifted versions of the CTFT. The spectra are shifted by the sampling frequency. That's the reason for aliasing.


Determine maximum frequency of input signal to make system LTI

This question, and a number of similar others, from DTSP book can be a little tricky to recognise the fact that it's not actually asking alias free operation (which would require 8 khz cutoff signal as you expected), rather it is actually asking how much aliasing (due to initial sampling block) is tolerable: you can allow aliasing in those regions of the DT filter $H(e^{j\omega})$ whose stop band will help you get rid of the aliased spectrum, hence enable an LTI processing possible, but not with the samples of the complete original fullband signal $x_c(t)$ but with that of the reduced bandwidth input signal.

To see this you shall better consult to the spectral figure below:
The mathematical condition on the maximum allowable cutoff frequency $\Omega_c$ of the input signal, whose aliased spectra will be deleted by the ideal lowpass filter with cutoff frequency $\pi/2$, is given as:

$$\Omega_s - \Omega_c > \frac{\pi/2}{T_s}$$
$$\frac{2\pi}{T_s} - \Omega_c > \frac{\pi/2}{T_s}$$
$$\Omega_c < \frac{2\pi}{T_s} - \frac{\pi/2}{T_s}$$
$$\Omega_c < \frac{3\pi/2}{T_s}$$

Thats stated in hertz frequency as:
$$f_c < \frac{3}{4}f_s$$

which gives a maximum allowable input signal bandwidth as $f_c=12,000$ Hz when the sampling rate is $f_s=16,000$ Hz.


Fractional Bandwidth of a Gaussian Amplitude Modulated Signal

The fact that it's modulated with a sinusoid doesn't change the FWHM bandwidth of your pulse – the $e^{jx}$ function has $\left\lvert e^{jx}\right\rvert\equiv 1$ at every point. That doesn't change the amplitude, so the FWHM of a sinusoid-modulated gaussian is just the same as of the unmodulated gaussian.


If you look at:

https://en.wikipedia.org/wiki/Matched_filter

the derivation does not depend on the filter being either causal or have a finite duration. The derivation is fully satisfactory for convolution for time index $k$ over $-\infty,\infty$.  It's all about what filter will maximize the Schwarz inequality.

In this case
$$
x(t)=\exp(-\frac{(t-0)^2}{2}) = x(0-t)
$$

Functions that are symmetric are their own matched filter.

The $T$ parameter enters the discussion in terms of when does one sample the filter for the purpose of detection.  If $x(t)$ was delayed by some unknown time(nonsynchronous detection),  I would use the matched filter to estimate it's delay.


Inverse fourier transform by using a random phase and an original magnitude (matlab)

You image is not an $N \times M$ matrix. It's actually an $N \times M \times 3$ array of numbers. I assume that's why you're using fftn instead of fft2 to take its Fourier transform.

You have to either

convert your image to grayscale; or
process each component independently.

To work get grayscale, we can naïvely average over the RGB colors like so:

img = mean(double(img),3)

To work with all three components independently, your phase array needs to have an extra dimension:

rand_phase= 4*pi*(rand(s(1),s(2), s(3))-0.5);


I suggest you to investigate Image Restoration using deconvolution since it is more appropriate for your application. Kalman filtering is designed to remove or ignore noisy data in a signal but it is not suited for deblurring (as far as I know).

Noting $\epsilon$ the additive noise, $H$ the convolutional mask (your Point Spreading Function), $x$ the original image you want to retrieve and $y$ the blurred image you have:

$y = H*x + \epsilon$

Deconvolution allows you to retrieve $x$ knowing $y$ and $H$.

Using Matlab you can use the following function:

deconvreg(image, PSF)

where PSF is an estimate of the Point Spread filter and image is the blurred and noisy image.

If you want to know more about the theory, this is called an ill-posed problem that you solve by minimizing a Least Square Error through an iterative process:

$$x_{estimated} = \arg\min_{x}||y - Hx||^2$$

Under some approximations this can be computed in a fast way through FFT.

This does not answer your initial question about Kalman Filtering but I hope it will help you! If you are interested I have documents about image restoration using inverse and Wiener filtering.