A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering

Basically they cascade moving-average filters, or CAF.
a Moving-average filter has the following coefficients
$h_{av}(m) = 1/M, m = 0, 1.. M-1$
In Matlab it would look like this
h = ones(1,M)./M;
y = filter(h,1,x);

If you want to cascade three moving-average filters
h = ones(1,M)./M;
y1 = filter(h,1,x);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

M should be equal to the number of samples per period. So if your frequency is approximately 50 Hz and your sampling frequency is 12800 Hz then M = 256;.
Edit : I modified your script to make it work. There were 2 problems, first of all you only processed 256 samples, your couldn't measure anything significant because the transient response of your 3 cascaded filters is 768 samples long. Secondly, you filtered the reference and not your complex signal...
clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
ADCsin = dsp.SineWave(1,50.1);
ADCsin.ComplexOutput = false;
ADCsin.SampleRate  = 12800;
ADCsin.SamplesPerFrame = 12800;
ADC = ADCsin();
plot(ADC);
grid on
hold on

%ref sin wave with 50Hz complex value
Refsin = dsp.SineWave(1,50);
Refsin.ComplexOutput = true;
Refsin.SampleRate  = 12800;
Refsin.SamplesPerFrame = 12800;
ref = Refsin();
plot(imag(ref))

%calculating  xs
xs = ref .* ADC;

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,xs);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

%calculating the final result
Fest = 50 - (angle(y3(end)) - angle(y3(end-1))) * 12800/(2*pi)


Why is the range of frequency for discrete time Fourier transform $-\pi<\omega<\pi$?

Short intuitive explaination: In a time discrete system, the spectrum repeats itself at multiples of the sampling frequency (normalized sampling frequency is $2\pi$). Because of that, the limits of the fourier integral can be reduced from $\pm \infty$ to $\pm \pi$.


I understand your confusion. What we are looking for is a function of $z$ which equals the squared magnitude $|H(e^{j\omega})|^2$ of a Fourier transform $H(e^{j\omega})=\mathcal{F}\{h[n]\}$ when evaluated on the unit circle $|z|=1$. Note that the inverse discrete-time Fourier transform of $|H(e^{j\omega})|^2$ is given by the convolution $h[n]\star h^*[-n]$, where $^*$ denotes complex conjugation. Since the $\mathcal{Z}$-transform of $h^*[-n]$ is given by $H^*(1/z^*)$ we see that the $\mathcal{Z}$-transform of $h[n]\star h^*[-n]$ is given by
$$\mathcal{Z}\big\{h[n]\star h^*[-n]\big\}=H(z)H^*\left(\frac{1}{z^*}\right)\tag{1}$$
If $h[n]$ is real-valued we obviously have $h^*[-n]=h[-n]$, and $(1)$ is equivalent to
$$\mathcal{Z}\big\{h[n]\star h[-n]\big\}=H(z)H\left(\frac{1}{z}\right)\tag{2}$$
It's important to understand that the equations
$$\big|H(z)\big|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{3}$$
and, for real-valued $h[n]$,
$$\big|H(z)\big|^2=H(z)H\left(\frac{1}{z}\right)\tag{4}$$
are only valid on the unit circle $|z|=1$.
In general we of course have
$$\big|H(z)\big|^2=H(z)H^*(z)\tag{5}$$
but $(5)$ is no valid $\mathcal{Z}$-transform of any time domain sequence, because it is not analytic. On the other hand, $(3)$ and $(4)$ are valid $\mathcal{Z}$-transforms, namely of the sequence $h[n]\star h^*[-n]$. On the unit circle they equal the squared magnitude of the Fourier transform of $h[n]$, i.e.,
$$H(z)H^*\left(\frac{1}{z^*}\right)=\big|H(e^{j\omega})\big|^2,\qquad |z|=1\tag{6}$$
Let's look at a simple example. With $h[n]=\delta[n]+\delta[n-1]$, we have $H(z)=1+z^{-1}$. Since $h[n]$ is real-valued we can use Eq. $(4)$:
$$H(z)H\left(\frac{1}{z}\right)=(1+z^{-1})(1+z)=z+2+z^{-1}\tag{7}$$
The actual squared magnitude of $H(z)$ is
$$\big|H(z)\big|^2=|1+z^{-1}|^2=1+2\textrm{Re}\{z^{-1}\}+|z|^{-2}\tag{8}$$
On the unit circle $z=e^{j\omega}$ both $(7)$ and $(8)$ equal
$$\big|H(e^{j\omega})\big|^2=2+2\cos(\omega)\tag{9}$$
However, the inverse $\mathcal{Z}$-transform of $(7)$ equals the inverse Fourier transform of $|H(e^{j\omega})|^2$, whereas $(8)$ is no $\mathcal{Z}$-transform.


I think you should first construct a smooth function of $r$ (where $r=\sqrt{x^2+y^2}$). For example, your original solution was based on the following function of one variable:

$$f_1(r) = \begin{cases} 1 & r < R \\
0 & \text{otherwise} \end{cases}$$

but you might try using a "softer" step function such as:

$$f_2(r) = \begin{cases} 1 & r<R-\frac{s}{2} \\
0 & r\ge R+\frac{s}{2} \\
\frac{1}{2} - \frac{r-R}{s} & otherwise
 \end{cases}$$

and you can create smoother variants, such as $$f_3(r) = 3f_2(r)^2-2f_2(r)^3$$ since $f_3$ is simply formed by mapping the linear ramp function of $f_2$ to smoother cubic shape.

The transition width of the ramp, $s$, should be set to $\frac{1}{2}$ or $\frac{1}{\sqrt2}$ if the linear ramp, $f_2$, is used, and you may need a larger value of $s$ (say, $s=1$) if the cubic function, $f_3$ is used.

Note that, strictly speaking, what you are trying to do is sample the original pillar-box function at discrete sample-points, so if you think about the function prior to discrete sampling (let's call this the pre-sampling function, $f(x,y)$ being defined over the Cartesian plane, for all real $x$,$y$)  the problem with your original pre-sampling function was that it contained a step-discontinuity (a $0^{th}$order discontinuity), which contains large amounts of high-frequency energy (in 2D spatial frequency), so you need to remove some of this high-frequency energy.

The radial function, $f_2(r)$, removes the $0^{th}$order discontinuity, resulting in far less high-(spatial)-frequency energy in the pre-sampling signal (but the resulting signal still has a $1^{st}$order disconinuity, since it's first derivative is not continuous). Likewise, the cubic function, $f_3(r)$, has even more reduced high-frequency energy, since it has no $0^{th}$ or $1^{st}$order disconinuities (but it's second derivative is discontinuous).

In general, if you make sure your pre-sampling function has it's first $n-1$ derivatives being continuous (and it's $n^{th}$ is disconinuous), then the amplitude of it's high frequencies will fall of like $\left(\frac{1}{f}\right)^{n+1}$

I hope this helps...


Why does my signal magnitude increase after adaptive filtering?

I don’t have the specific details for your filter but with digital filters in general it is typical for the filter to grow the signal in band in contrast to analog filters that shrink the signal out of band. It is all just a matter of scaling. Consider the simple case of a moving average FIR filter consisting of the summation of the previous N samples; such a filter will grow any low frequency signal within its passband, specifically growing a DC component by a factor of N. (We could then divide the result by N to get the actual average of try samples, but the primary filter itself here is growing the signal).

This approach of considering digital filters as "growing a signal in band" is particularly important for noise considerations in fixed point design, where you should typically avoid scaling the signal prior to the filter or the filters coefficients in order to normalize the result, but always allow the filter to grow the signal and then scale afterward. This is the reason for extended precision accumulators, and the reason is rather straightforward: if you scale the signal prior to the filter in a fixed point design you are effectively increasing the quantization noise. This quantization noise is typically modeled as independent white noise from sample to sample so in the filter will increase at a rate of $10Log_{10}(N)$ for the $N$ taps in the filter. If you scale after you are only adding this higher quantization noise level once.  Scaling the coefficients can be shown to accumulate quantization noise contributions in a similar way. This is also a reason why equalizers and adaptive filters should not be any longer than they need to be (to avoid noise enhancement).


In contrast to common misunderstanding, aliasing does not affect the process of sampling, rather it affects the process of reconstruction. i.e., the samples themselves do not contain error, but their interpretation at the reconstruction (or anything related with it) will be in error.
So, if you sample a sine wave $$x(t) = A \sin( 2 \pi f_0 t + \theta) $$ at a sampling frequency of $F_s = f_0$, (which is undersampled according to Nyquist-Shannon sampling theorem for bandlimited signals), then your samples will be $$ x[n] = A \sin(\theta) $$ a constant signal! When reconstructed back into continuous waveform at the interpolator output, you will get a DC signal, $x(t) = A \sin(\theta)$, as opposed to the original sine wave.
The error is due to the false assumption at the interpolation stage that those samples $x[n]$ came from an original continuous-time signal which was band-limited to $f_0/2$, Nyquist frequency of sampling frequency $f_0$, which is not the case in reality.


Is there a name for the procedure of taking the FT over separate consecutive small time-blocks?

The Short Term Fourier Transform (STFT)

https://en.m.wikipedia.org/wiki/Short-time_Fourier_transform

.

For your question, the windows are finite and non overlapping.

The continuous time case is probably more of a mathematical device that could be approximated by diffraction over some periodic, channelized structure.

The discrete time case is very commonly used, although as RJB pointed out, usually with overlap.

If your signal consisted of flat noise and a set of sinusoids, for large DFT segments, the bins would be approximately statistically independent in time and frequency for a non overlapping window, so there are applications where overlapping is not used.


how to deal with signal to make it suitable for kalman filter

Kalman filter have basic requirements that the random processes like input signal, measurement signal,process noise and measurement noises all should be white Gaussian in nature.what Kalman filter actually do is that it uses the system model and input signal and generates its own output which is actually a predicted value now using the observation from the system it calculate the correction to be applied to the predicted output and hence final results are better than observations obtained from the system. Hence very first step- is that you need a system model. i am totally disagree with @Deniz notion that in non-linear case you should use particle filter. particle filter is used when the random process as mentioned above are not white Gaussian. Even if your system is non linear you can use Extended kalman filter or unscented kalman filter. As ur case is concerned i think you should apply extended or unscented kalman filter as velocities for navigation will keep on changing.


minimum requirement of SNR for fft algorithm

Some assumptions are made in this paragraph. The signal is sinusoidal. It lands exactly at a frequency domain bin, and has magnitude $a$. No frequency domain interpolation is needed. Frequency estimation will look for and report the bin with the largest magnitude. This can be the signal bin or a noise-only bin. The noise is additive, but its contribution to the signal bin will be ignored as an approximation. The real $\operatorname{Re}(Z)$ and imaginary $\operatorname{Im}(Z)$ components of the noise $Z$ in the noise-only bins are independent and identically distributed (IID) Gaussian random variables of zero expected value.

The probability that the magnitude of the noise, at a given noise-only bin, is lower than $a$ is given by the cumulative distribution function of the Rayleigh distribution:

$$P(|Z| < a) = 1-\operatorname{exp}\left(\frac{-a^{2}}{2\sigma^2}\right),$$

where $\sigma^2 = E\left(\operatorname{Re}(Z)^2\right) = E\left(\operatorname{Im}(Z)^2\right)$, the expected value of both the squared real and imaginary components of $Z$. If there are $N$ bins to look in, then the probability that none of the noise-only bins has higher magnitude than the signal bin is given by:

$$P(|Z| < a)^{N-1}.$$

A conclusion is that the correct bin will be chosen at a higher probability if you can narrow down the frequency range where to look for the signal, making $N$ smaller. It is up to you/someone to decide how much risk of reporting the wrong bin is acceptable. There is more to say about this I'm sure.


How to find out if a transform matrix is separable?

I admit I did not really thought about it before. I hope my notations won't be too sloppy.

I assume that given an operator matrix $A(u,v)$, you can apply this operator as a transform on an image $I$, to obtain an image in a novel domain $J(u,v)$. For instance a Fourier kernel would give $$a_{m,n}(u,v) = \exp^{-2\pi \imath \left(um/M+vn/N\right)}\,.$$

For each choice of a $(u,v)$ pair, you get a specific kernel, for instance a basis element.

The transform is separable, in the common sense, when you can write it as a element-wise product, for instance:
$$a_{m,n}(u,v) = b_{m}(u) . c_{n}(v)  $$
where $b(u)$ and $c(v)$ are 1D  vectors:
$$A(u,v) = b(u) . c^T(v)\,,  $$
one acting on  rows, the other on columns.

Thus, the vectors being of rank at most one (the rank of a matrix is lower than its minimum size), $A$ is at most of rank $1$ as well. I take it as sufficient condition. If I rule out uninteresting matrices and vectors of rank $0$, that should be a necessary condition. More details for instance at A rank-one matrix is the product of two vectors.

So:

if $A$ is of rank $>1$, it does not appear to be separable (at least with this definition)
if $A$ is of rank $1$, take a non-zero vector for the span, and then multiply it by appropriate coordinates to get the matrix back.

For $A_3$, the determinant is zero, so the matrix is rank $1$. You will get vector $[1 -1]^T$ in the image, and quickly:

$$2 A_3 = [1 -1]^T [1 -1]\,.$$

Note, the separation is not unique.


Why is the received signal=information+RF signal+noise,not just information+noise?

Well, it seems the authors assume a setup where there are separate data receivers (DR) and energy receivers (ER) ("To achieve a balanced user experience for separately located datareceivers (DRs) and energy receivers (ERs) in the network, joint transmit beamforming vectors are optimized [...]"). In that setup, $s_k^E$ would be symbols that are used for sending energy to the ERs whereas $s_k^D$ are data symbols for the DRs.

Therefore, if you are a DR and want to receive data, you care about two types of interferences: data symbols for other DRs and energy symbols for the ERs. Note that you didn't copy equation (1) exactly: the second term should be $s_i^D$ (data symbol interference). Then you have the third term with interference of energy symbols for the ERs and the fourth term for the noise.

As an ER, you wouldn't care about interference, it's actually additional power, if you like.

Other papers use other models where one node may want to receive data and harvest energy at the same time. As they say: "all models are wrong, but some are useful..."


Assuming your 512 samples of the signal are taken at a sampling freqeuncy $f_s$, then the resulting 512 FFT coefficients correspond to frequencies:

0, $f_s/512$, $2 f_s/512$, $\ldots$, $511 f_s/512$

Since you are dealing with discrete-time signals, Fourier transforms are periodic, and FFT is no exception.

Therefore you can think of your last coefficient as corresponding also to frequency $511 f_s/512 = (511-512) f_s/512 = -1 f_s/512$.

The same applies to the second to last coefficient, and so on. This is the mirroring commented by Daniel Hicks.

Also, if you are transforming a real signal, then all your information is contained in the first 256 FFT coefficients. The rest are simply complex conjugates of the first coefficients.


Detecting expected patterns in repeated measurements

So far I learnt how to use basic filters to get rid of noise (moving average, exponential smoothing)

Those are all low-pass filters!

  checking sudden changes (difference between two successive filtered values)

That is a high pass filter !

To be honest, I don't see "patterns" in your signal, the peaks are all pretty different, but you don't care for "patterns", either:

You really are just interested in the signal that first passes through a low-pass filter, and then a high-pass filter.

Combined, that is a band-pass filter.

So if I'd recommend something: Look into filtering; understand the spectrum of signals, and understand how to build filters that extract specific signals from the original signal.


Direct downconversion (IQ downconverter) with ASK modulation-- what kind of Balun do you need?

You could make your intentions clearer about "using direct downconversion for carrier interference cancellation" to "move the carrier interference to DC then DC block" in your scheme of "[d]irect downconversion ... with ASK modulation", had you attached the circuit drawing.
But, for a prescribed application of MAX2021 as GSM/CDMA signaling demodulator, you need a power combiner -- in fact, three power combiners: two for differential components in the $I$ and $Q$ paths and one to combine the single-ended $I,Q$ signals into one baseband signal. The "power" determinant should not confuse you: the name is inherited from a nomenclature of microwave power devices (e.g., microwave power amplifiers, attenuators, etc.), but, when used in the context of modulators/demodulators, designates exactly the component you need to combine $I_{+}/Q_{+}, I_{-}/Q_{-}$ outputs into single-ended signals. Besides, microwaves101.com directly indicates a correspondence between baluns and power combiners in the section 180 degree combiners.
Also, the baseband output of MAX2021 is not "square pulses" with zero rise/fall times; after output filtering according to a circuit shown in Figure 4. Baseband Port Typical Filtering and DC Return Network, baseband outputs look more like an "analog" signal with a frequency range defined by the symbol data rate.
In Figure 5. Demodulator Combining Diagram (page 15), the MAX2021 datasheet of your reference directs you to use Mini-Circuits' 180° power combiner ZFSCJ-2-1 and 0° power combiner ZFSC-2-1W-S+. ZFSCJ-2-1 features wideband operations, 1MHz to 3GHz. If your symbol data rate requires frequency below 1MHz, there are available on the market lumped-element/resistive power dividers/combiners of other suppliers, with a lower frequency value down to 4KHz and even DC.
Read also the MiniCircuits blog article and the application note AN10-006.


Fast method to decompress huffman and arithmetic encoded data

For actual implementations take a look at (C) programming codes for decoders.

Basically the process is simple, you create an abstraction such as

bits = get_bits(n)

Where, n is the n bits you want to read. the var bits is left aligned. It extracts it from some byte_buffer - say a 32 bit or 64 bit symbols.

  NOTE: All, practical codecs put major symbols at 'Byte aligned' even though >internal symbols can be bit by bit.

Now for the implementation efficiency part, basically, irrespective of what huffman code will come out to be, (and hence whatever be the required n).

this is how simple get_bits can be implemented:

// At start
frame_buffer = 0x4444
bits_delivered = 0

n = 1
mask = 0x8000
bits = (mask & frame_buffer) >> 32-n
mask = mask >> n
bits_delivered += n
if(bits_delivered == 32)
    { load_framebuffer( &frame_buffer) }

(Actual code needs to take care of few more stuff boundary conditions)
The actual data is first fetch into frame_buffer which is a 32 bit register inside ALU (actually 16/32/64 bit depending on type of CPU) and then rest of the processing is only inside ALU's registers only.

Only when you exhaust the frame_buffer, do you really flush that bit-buffer and get another integer from memory.

This is how it gets its efficient implementation. Hope this answers your question.


FFT doesn't always have frequency in all frames

You might be looking only at the cosine component (sometimes called the even or real component) of the FFT result, instead of the magnitude of each entire complex FFT result vector.  The real component will vary as your test signal goes in and out of phase with the cosine basis vectors (always even functions) of each frame, even if the complex bin magnitude stays constant (or varies only slightly due to windowing scalloping).

mag[i] = sqrt(real[i]*real[i] * imag[i]*imag[i]);


Quick question about difference between two implementations of MFCC

At a glance - the number of filters is different, the first algorithm has it equal to 48, while the second one - to 32. But even these small things aside, MFCC computation is basically a pipeline of algorithms with 4 or 5 stages. And each stage has some input parameters varying which you can get slightly different output. Given that each stage is implemented properly these differences have negligible effect on the recognition result.


State-Space Representation of Forward and Backward Filters

The backwards filter, as the name suggests filters the input signal with the same filter $(F,G,H,D)$ backwards in time

$$
\begin{align}
x(t-1) &= F\,x(t) + G\,u(t) \\
y(t) &= H\,x(t) + D\,u(t)
\end{align} \tag{1}
$$

But in order to define a filter using a state space model, the update law does have to go forwards in time, which can be done by solving the first equation for $x(t)$

$$
x(t) = F^{-1}\,\left(x(t-1) - G\,u(t)\right) \tag{2}
$$

however $x(t)$ can only a function of $x(t-1)$ and $u(t-1)$. This can be solved by using a new definition for the state $z(t) = x(t) + F^{-1}\,G\,u(t)$, so

$$
z(t) = x(t) + F^{-1}\,G\,u(t) = F^{-1}\,x(t-1) = F^{-1}\,\left(z(t-1) - F^{-1}\,G\,u(t-1)\right) \tag{3}
$$

which can be rewritten in the following more standard form

$$
z(t+1) = F^{-1}\,z(t) - F^{-2}\,G\,u(t). \tag{4}
$$

Substituting the definition of $z(t)$ in to the expression for $y(t)$ gives

$$
y(t) = H\,x(t) + D\,u(t) = H\,\left(z(t) - F^{-1}\,G\,u(t) \right) + D\,u(t) \tag{5}
$$

which can be simplified down to

$$
y(t) = H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t). \tag{6}
$$

Combining equation $(4)$ and $(6)$ gives the following (forwards) state space model, which is equivalent to the (backwards) state space model in equation $(1)$

$$
\begin{align}
z(t+1) &= F^{-1}\,z(t) - F^{-2}\,G\,u(t) \\
y(t) &= H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t)
\end{align} \tag{7}
$$

So the backwards filter based on the forwards filter $(F,G,H,D)$ can be written as $(F^{-1},-F^{-2}\,G,H,D - H\,F^{-1}\,G)$.


A question regarding z transform and its magnitude response

The frequency response of a given filter is calculated evaluating $H(z)$ in the unit circle, i.e. $z=e^{j\omega}$.

Hint 1:

Note that you can express your transfer functions as follows:
$$H(e^{j\omega})=\frac{1-a^*e^{j\omega}}{e^{j\omega}-a}=\frac{1}{e^{j\omega}}\cdot\frac{1-a^*e^{j\omega}}{1-ae^{-j\omega}}$$

Hint 2:

There is a property that states that given any complex number $z$, the following equality holds:

$$|z|=|z^{*}|$$

Can you solve it now?


Discrete time inverse fourier transform of cosine squared

The only thing that went wrong is your last statement:

They're only zero if $k\ne 0$.

Thus, you get three Dirac delta distributions.

Let's look at your last equation, and work with that:

$$
\begin{align}
 x(n) &= \frac{1}{2\pi}\int\limits_{0}^{2\pi} \frac{e^{j\omega(n+2)} + e^{j\omega(n-2)} + 2e^{jn\omega}}{4}\,d\omega \tag1\label{orig}\\
&=\frac1{8\pi}\left(\int\limits_{0}^{2\pi}e^{j(n+2)\omega}\,d\omega + \int\limits_{0}^{2\pi}e^{j(n-2)\omega}\,d\omega + \int\limits_{0}^{2\pi}e^{jn\omega}\,d\omega \right)\tag2\label{linearint}
\end{align}$$

From $\eqref{orig}$ to $\eqref{linearint}$ I just went ahead and used the fact that integration is a linear operation to decompose that single integral into a sum of three, very similar, integrals.

Now, let's look at the last one, as you (sadly, with a mistake), already did:
$$
i(n):=\int\limits_0^{2\pi}e^{jn\omega}\,d\omega =
\begin{cases}\frac1{jn}e^{jn\omega}\Big|_{\omega=0}^{2\pi} \overset{e^x\text{ is }j2\pi\text{-periodic}}= 0& \text{for } n\omega \ne 0\\
\int\limits_0^{2\pi}e^{0}\,d\omega= \int\limits_0^{2\pi}1\,d\omega= 2\pi& \text{else.}
\end{cases}
$$

So, this integral yields zero for all values of $n$ but for $n=0$, in which case it yields a fixed, finite, non-zero value.

Your other two integrals are just the same, but shifted w.r.t. $n$ by $2$ to the left or the right.

For more info on the Dirac delta distribution (looks like a function, isn't quite exactly a function), I'd recommend your favourite textbook or wikipedia. It's a thing that, when integrated, has an integral of $1$, but in an arbitrarily narrow support around $0$, outside of which its value is $0$.


Which filters are the most commonly used and easier to implement for ECG noise removal?

Different filters are used for different purposes. Removing powerline noise, which usually is concentrated at a single frequency, is usually performed with a notch filter. Removing baseline wandering is usually done with a DC blocking filter, which is a very narrow high pass filter, and is essentially another notch filter with 0 Hz center frequency. Removing EMG is usually more difficult, but you can try doing it with a low pass filter, since EMG has a lot of content in high frequencies.

In order to process EEG signals for simple things like heartbeat rate detection, look at the Pan-Tompkinks algorithm.

Questions 2 and 3 require some effort on your part, since they're very broad. I've used this book in a graduate class on the topic, and I found it to be very insightful.


Can I set a constraint on the first tap of an FIR filter such that its inverse is stable?

The Invertible FIR Filter

A constraint based on the first coefficient alone is developed as follows:

From Cauchy's argument principle any FIR filter that meets the following constraint will be invertible (including marginal invertibility, change $\le$ to $<$ otherwise):

$$\max\left(\arg \left( H(e^{j\omega}) \right)\right)-\min\left(\arg \left( H(e^{j\omega}) \right)\right) \le \pi\tag{1}\label{1}, \space\space \omega \in [\omega_o, \omega_o+\pi) $$

Where:

$H(e^{j\omega}) = \sum_{n=0}^{N-1}h_ne^{-j\omega n}$: Frequency response of FIR filter

$\arg()$ : unwrapped phase of $H(\omega)$

Note: In an earlier version of this post, I had a simpler subset of the above that constrained the absolute value of the phase to not exceed $\pi$, but consider given any such solution that we can rotate the filter response by a fixed angle without effecting invertibility. Thus it is more generally constrained to be that for a FIR with real coefficients, the difference of the unwrapped phase for $\omega \in [0, \pi)$ cannot exceed $\pi$. Extending this to an FIR with complex coefficients, and given that we rotate the filter response (which shifts the frequency response by a fixed amount), results in the generalized constraint above applicable to any interval over $\pi$ in frequency.

Thus to constrain the first coefficient $h_o$ we can derive from $\ref{1}$:

$$\arg( h_0) + \max\left(\arg\left(\sum_{n=1}^{N-1}h_ne^{-j\omega n} \right)\right)- \min\left(\arg\left(\sum_{n=1}^{N-1}h_ne^{-j\omega n} \right)\right) \le \pi\tag{2}\label{2}$$

Which shows the complexity of a constraint based on the first coefficient alone but that it can exist. It is equivalent from this and simpler to state that if all the zeros are inside the unit circle (an invertible FIR filter), then the plot of the frequency response on a complex plane as we sweep $\omega$ from $0$ to $2\pi$) cannot encircle the origin.

Further details below:

An invertible FIR filter is a minimum phase filter, since all zeros must be inside the unit circle (or on the unit circle for marginal stability). However a subset of all possible minimum phase filters given a first tap would be the set of decreasing coefficients as the OP has hypothesized as the would meet the constraint given below under “by inspection”.  However other minimum phase filters also exist where subsequent taps are larger than the first.

The easiest constraint I can think of that would be on all coefficients specifically beyond solving for the roots is given by Cauchy's argument principle: The frequency response as given by the coefficients cannot encircle the origin for an invertible FIR filter. With the frequency response given as:

$$H(z) \bigg|_{z=j\omega} = \sum_{n=0}^{N-1}h_n e^{-j\omega n}$$

Further details below:

The invertible filter must have all zeros inside the unit circle, since zeros become poles once the filter is inverted, and any poles outside the unit circle means instability for causal systems. (For marginal stability consideration; meaning a system that neither grows nor decays, the zero can be on the unit circle.)

An FIR filter with all zeros inside the unit circle is a minimum phase filter Any other FIR filters would not be invertible. This includes maximum phase filters, which have all zero's outside the unit circle, and mixed phased filters with both minimum and maximum phase components (linear phase filters are mixed phase).

So the constraint is the filter must be a minimum phase filter to be invertible. Below I list 4 tests for detecting a minimum phase filter, with the Cauchy's argument principle the closest to providing a simple coefficient constraint rule.

By Inspection: A tell-tale sign of a minimum phase filter is the concentration of coefficients toward the start of the filter. Given this, you can rule out many filters by simple inspection if the coefficients are either concentrated in the center or end of the filter. Specifically when considering all filters with the same magnitude response, the coefficients for the minimum phase filter (which is the filter's impulse response), will decay the fastest in time. A detailed proof of this property specific to minimum phase polynomials is given in Oppenheim and Shafer's Digital Signal Processing book and is summarized as:

$$\sum_{n=0}^N|h[n]|^2\ge \sum_{n=0}^N|g[n]|^2$$

Where $h[n]$ is a minimum phase filter and $g[n]$ is any other filter with the same magnitude response and $N$ can be any positive integer. This does not mean all coefficients are in decreasing order, for example [5 6 3 2 1] is linear phase, while [5 8 3 2 1] is not, so this is not necessarily a simple constraint that can be applied but can certainly identify obvious non-linear phase solutions.

Cauchy's Argument Principle: A very simple approach to test for this condition applicable to FIR filters is to use Cauchy's argument principle (see Nyquist's Stability Criterion and Cauchy's Argument Principle) by plotting the frequency response on a complex plane. For a causal FIR filter, the number of clock-wise encirclements of the origin will be equal to the number of zeros outside the unit circle. If all the zeros are inside the unit circle, there will be no encirclements of the origin (I show an example below). For an FIR filter, there cannot be any counter-clockwise encirclements (as all poles are at the origin), so if any encirclements do happen, they will be clock-wise only.

Solving for the Roots: Confirm the magnitude of all roots of the polynomial given by the coefficients of the filter are all $|z|\le 1$

Hilbert Transform: Compare the magnitude and phase response for the filter. Since there is a unique relationship between the magnitude response and phase response for minimum phase filters, the two can be compared to determine if the filter in question is indeed the minimum phase solution for that magnitude response. This is detailed further by PeterK in this post: Derive minimum phase from magnitude with relationship copied below:

$$
\theta(\omega) = - {\scr H}\left[ \ln(G(\omega)) \right]
$$
where ${\scr H}$ is the Hilbert transform.

Every magnitude response has a minimum phase solution, therefore every FIR filter can be decomposed into a minimum phase filter (invertible) cascaded with a all-pass filter (constant magnitude response with phase change only over frequency and not invertible).

These concepts are demonstrated with a simple example of a 2 tap FIR filter with coefficients [1 0.5] and it's reverse [0.5 1]. In the first case the filter is minimum phase with a transfer function $1+0.5z^{-1}$ and the second filter is the reverse, which is a maximum phase filter with the transfer function $0.5+z^{-1}$. The magnitude response of both filters is identical but the phase response is very different, as given by the vector diagram showing phase versus frequency for both filters. (This diagram is created by simply replacing $z^{-1}$ with the phasor $e^{-j\omega}$ which is exactly how to get the frequency response from the z-transform).

It is a bit of an optical illusion but the resulting phasor given by the sum $\sum_n h[n]e^{-jn \omega}$ does have the exact same magnitude for all frequencies as we sweep the frequency $\omega$ from $0$ to $2\pi$. However notice how constrained the angle will be in the diagram on the left, compared to the one on the right! Minimum phase versus Maximum phase.  This plot also demonstrates Cauchy's argument principle showing how the frequency response of a minimum phase filter when plotted on the complex plane will not encircle the origin.

Below is a plot of the magnitude and phase response for the two example filters above (with the same magnitude response for either). Since group delay is $d\phi/d\omega$ the minimum phase filter will have the lowest delay while the maximum phase filter will have the largest delay, which makes sense when you consider the placement of the largest tap in the FIR filter which is a series of summed and weighted delays (the energy will emerge from the minimum phase filter sooner).

Another example demonstrating the Cauchy argument principle is below, with a plot of the frequency response on the complex plane for filters [5 6 3 2 1] and [5 8 3 2 1]. The filter [5 8 3 2 1] is proven to not be minimum phase since the frequency response encircles the origin:

And here is another example for a maximum phase filter, which is also confirmed using Cauchy's Argument Principle. This is for the fourth order filter with coefficients [1 -3 -3 2 5] where we see all four zeros are outside the unit circle as we have four encirclements of the origin. (The easy way to count encirclements is to note a direction on the frequency response with a forward direction consistent with increasing $\omega$, and then draw a vector from the origin out toward infinity at any angle and count how many crossings of the frequency response take place: if the cross is of a forward direction the count increases, and if of a negative direction the count decreases).

And here is another simple example with the FIR filter given by coefficients [1, 1,5, 0.6] showing how the frequency response of an FIR with a real positive first coefficient can enter the LHP and still be a minimum phase filter. Specifically, it is the origin that is not encircled consistent with Cauchy's Argument Principle. With that constraint the phase response can never exceed $\pm \pi$. Below the plot is the standard magnitude and phase frequency response as two separate plots.


Finding the frequency response of a filter defined as a Z-domain transfer function

The analytic way is to substitute the variable $z$ by $e^{j\omega}$ to get the frequency response $H(\omega)$ (with $\omega = \frac{2 \pi f}{F_s}$) - that is to say, the frequency response is the $z$ transform evaluated on the unit circle.

Note that matlab has a built-in function for plotting the frequency response straight from filter coefficients (freqz), whose source code can be seen. Have a look! The implementation does not even involve solving the difference equation or filtering an impulse - you can simply divide the FFT of the numerator sequence by the FFT of the denominator sequence - padded up to the desired number of frequency points.


How to extract traffic signs from a photograph?

Did you try something simple like correlation?

(EDIT). The idea behind correlation is to use a template (in your case a trained road sign sample), and compare it to every position in the test image. The comparison operation I've used to generate the images below is called normalized cross-correlation. Roughly speaking, you standardize (mean=0, standard deviation=1) the pixels in the template and the image part you want to match, multiply them pixel by pixel, and calculate the mean value of the products. This way you get a "match score", i.e. a measure of similarity between the template and the test image at every position in the test image. The position with the best match (highest correlation) is the most probably candidate for the position of the road sign.  (Actually, I've used the Mathematica function CorrelationDistance to generate the image below, which is 1 - (normalized correlation). So the darkest spot in the match image corresponds to the best match).

I don't have any other templates, so I simple cropped the sign from the second picture you posted:

Even though the template is rotated slightly, cross correlation still looks usable

and the best match is found at the right position:

(You'd need multiple scaled versions of each template to detect signs at any size, of course)


The selection of the amplitude of the signal, consisting of rectangular pulses, at different frequencies

The following is what I believe to be an optimized approach for performance in the presence of additive white noise when no other information is known about quantity of pulses or their amplitude distributions, beyond that they are 100 or 200 us long rectangular pulses and repeat at the 1KHz and 1.4KHz rates. This can be even further improved if any other further information is known (including learning from the data to see if anything else about the waveform is repetitive and therefore can be predicted once identified).

Note that the composite waveform consisting of the two pulse types combined  (pulse type meaning 1KHz and 1.4KHz pulses) repeats every 5 ms.

Create a mask to establish the timing of the pulse types. If the starting time between the two pulses is fixed and known then you can create a single reference waveform that is 0 everywhere except where a pulse could be located, and 1 in those positions using 100 us assumed durations. The longer the mask the better the result up to the maximum expected duration of a burst of pulses but keep it in 5 ms durations to take advantage of the cyclical property for using circular cross correlation.

If the 100 and 200 us pulses are equiprobable then in would actually be slightly better to use 200 us pulses in the mask (the break-even is when the ratio of 200 us pulses to 100 us pulses is $\sqrt{2}$, since 200 us masks would have 3 dB more noise which is offset by the signal increase if the mask allows the accumulation of $\sqrt{2}$ more pulses.)

Perform a circular cross correlation of the mask with a capture of the data. The result will have cyclical peaks showing the starting position of the 5ms composite pattern.

The circular cross correlation is given by the following which could be used as a quick test of the procedure:

$$XCORR = ifft(fft(b)conj(fft(a))$$

Where fft(a) is the fft of your mask and fft(b) is the fft of your data.

Here is a very simple example of that to demonstrate how this is done properly and the expected result since the OP is having trouble with this step:

>> mask = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];
>> // create three cases of test data for example
>> dataAlligned = mask;
>> dataShiftOne = circshift(mask,1,2);
>> dataShiftTwo = circshift(mask,2,2);

Thus each would be:

mask = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];

dataAlligned = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];

dataShiftOne = [0 1 0 0 0 0 1 0 0 0 0 1 0 0 0];

dataShiftTwo = [0 0 1 0 0 0 0 1 0 0 0 0 1 0 0];

The circular cross correlation is done for each case above as:

>> corr = @(data,mask) ifft(fft(data).*conj(fft(mask)));

Below shows the result of this for each case:

>> corr(dataAllgined, mask)
ans =
   3   0   0   0   0   3   0   0   0   0   3   0   0   0   0
>> corr(dataShiftOne, mask)
ans =
   0   3  -0   0   0   0   3  -0   0   0   0   3  -0   0   0
>> corr(dataShiftTwo, mask)
ans =
   0   0   3   0  -0   0   0   3   0  -0   0   0   3   0  -0

Given the sparse nature of the mask, the implementation for the circular cross correlation would likely be more efficient in the time domain by repeated multiply and sums at all circular rotations of the mask. Given the mask only has values of 0 and 1 this would be quite simple to compute even in a microcontroller, as the multiply operation would just be a selection of which samples to sum.

If the starting delay between the two pulse types is not known, then do the above with two masks; one for the 1KHz pulses and one for the 1.4KHz pulses.

Once timing is established for the above techniques, then we know when a pulse of each type should appear if a pulse is to be present in that interval. So at this point to detect each pulse use an integrate and dump that would integrate over 100us in two successive durations at all the expected start times. (This can be implemented with a moving average filter that has a width equal to the pulse duration, and taking the samples of the moving average output at each expected end of pulse location.)  A decision metric is established using a noise optimized threshold on the integrate and dump output to decide pulse present=1 or pulse not present = 0. The use of two successive durations is to determine if the pulse was 100us long or 200 is long (decoding 1,1 or 1,0 for that interval).

The magnitude of each pulse is the magnitude of the integrate and dump output scaled by the number of samples—- keep a running average of the detected magnitude for each pulse type and adjust the threshold based on the averaged magnitude to optimize the decision for pulse present or not present.


Difference between power spectral density, spectral power and power ratios

I have no idea what your calculation of power spectral density gives since
I cannot understand it.

If a signal $x(t)$ has Fourier transform $X(f)$, its power spectral density
is $|X(f)|^2 = S_X(f)$.  The absolute spectral power in the band of frequencies from $f_0$ Hz to $f_1$ Hz is
the total power in that band of frequencies, that is, the total power
delivered at the output of an ideal (unit gain) bandpass filter
that passes all frequencies from $f_0$ Hz to $f_1$ Hz and stops everything
else.  Thus,
$$\text{Absolute Spectral Power in Band}
= \int_{-f_1}^{-f_0} S_X(f)\,\mathrm df + \int_{f_0}^{f_1} S_X(f)\,\mathrm df.$$
The relative spectral power measures the ratio of the total power in the band
(i.e., absolute spectral power)
to the total power in the signal.  Thus,
$$\text{Relative Spectral Power in Band}
=  \frac{\displaystyle\int_{-f_1}^{-f_0} S_X(f)\,\mathrm df
+ \int_{f_0}^{f_1} S_X(f)\,\mathrm df}{\displaystyle\int_{-\infty}^{\infty} S_X(f)\,\mathrm df}.$$


Radon Transform reducing the number of parallel beams (MATLAB)

Unfortunately you can't do this with the radon function as it comes.  However, the function is written in MATLAB form (as opposed to being a compiled function), so you could extend the function yourself.  This would also be a good way to learn more about the process.

To reconstruct properly, you would also need to perform edits on the iradon function.

An alternative way to get what you're looking for would be to get projections as they come from radon, then reduce the array size by averaging the projection rays in groups of 4.  You would probably want the reconstruction to be the same size, so you'd need to repeat the averaged values 4 times in the argument to iradon.  So maybe something like this:

p = phantom ();
proj = radon (p, 0:179);
proj_avg = proj (1:4:end,:) + proj (2:4:end,:) ...
         + proj (3:4:end,:) + proj (4:4:end,:);
proj_rep = reshape (repmat (proj_avg',1,4)', [], length (proj(:,1)))';
recon_few_beams = iradon (proj_rep, 0:179);


Do you need extract the spectral envelope from original and pitch shifted signal and then warp the shifted spectral to match with the original!

You can use cepstrum, LPC or a channel vocoder to extract the envelope!

The real cepstrum can be computed using this equation:

Rcepstrum=real(ifft(log(abs(fft(windowed_signal)))))

Now do you need low pass Rcepstrum

And finally apply FFT again to get the Spectral Envelope, you will get a smoothed version of the spectrum!

At this point you will need get the diff spectral shapes from the (Shifted-Original), this is your correction factor, now apply the correction factor in the shifted spectral.

Edit

Let me try show how it is done:

FFT1 = Original
FFT2 = Pitch shifted

Get the Real Ceptrum from both signals:

cepstrum1 = real(ifft(log(abs(FFT1))));
cepstrum2 = real(ifft(log(abs(FFT2))));

Low pass window from cepstrum:

Choose some index value to cut the cepstrum (50, 100, 150), get just the firsts indexes (0 to 50 for example) and put zero at all.

Extract the envelope from the slashed window

envelope1 = exp(2*real(fft(cep_cut1))); % spectral shape
envelope2 = exp(2*real(fft(cep_cut2))); % spectral shape

Find the correct factor from envelopes:

warp_factor=(envelope1/envelope2)

Apply the correction:

new=(FFT2.*warp_factor)

Back you signal to time domain again (not forget, do you need overlap and add)

New_Formant_signal = real(ifft(new))

My small tests some times give me poor results, maybe you can use some peak picking in the FFT and interpolate the results to give a nice envelope, I will try one day, bu I never have time lol

Of course if you are pitch shifting monophonic singnals you can use PSOLA to keep the formants Intact...


How to recover $f(t)$ from Fourier Transform of its absolute value $\mathcal{F}|f(t)|$?

Matt L.'s comment said it. Once you take the absolute value of the function, you lose its sign. Whatever you do to the result, like taking the Fourier transform, doesn't help. The information is lost.

It is not a part of your question, but the function could have known properties that enable recovery. If it is band-limited, like your cosine, then recovery is possible up to a sign ambiguity $f(t)$ vs. $-f(t)$. In other words, if you know the absolute value of a band-limited function at all times $t \in \mathbb{R}$, and if you know the sign of the function at time $t_1$, then you also know the sign of the function at time $t_2$. For an algorithm, see Gaurav Thakur (2010) "Reconstruction of Bandlimited Functions from Unsigned Samples".


After measuring a NFC tag with a USRP and making a couple of simulations in octave and python, it seems like I was incorrect regarding the point-wise multiplication. As it stands, I believe the best simulation is to do an addition of the carrier wave (13.56 MHz) and the subcarrier wave (847 kHz), with the subcarrier wave's amplitude set at 0.1x the amplitude of the carrier wave.


DC and Nyquist values are higher than all other bins

For a mean centered (DC removed, pure AC) signal the first bin of the FFT output should be an exact zero. In your plot, this zero output sample is possibly obscured by the squeezing of the graph display; i.e., if you zoom into the left side of your plot, you should see that zero valued first bin sample...

Furthermore, for a typical FFT computation of real data, the second half of the FFT bins, is a conjugate symmetric mirror of the first half (excluding the DC and / or the Nyquist bins). Note that the Nyquist bin is (about) at the center of the FFT output, and to the right of this Nyquist bin is the mirrored content. Note also that exact Nyquist frequency bin exists only for even length FFTs and not for odd length FFTs.

Looking into your output, the last and the first bins seem like mirrors; but that's wrong as I have indicated, your fist FFT bin must be a zero, and even if it's not zero, it won't be mirrored to the last bin. What's mirrored to the last bin is actually the second bin of the FFT output.


Is maximizing hamming distance the same as minimizing correlation?

Two codewords $c_1$ and $c_2$ of length $n$, with elements in $\lbrace +1, -1 \rbrace$, and Hamming distance $d$, have a cross-correlation given by $$(n-d) -d = n-2d.$$ The reason is that there are $n-d$ bits that are equal and their product is $1$, and $d$ bits that are different and their product is $-1$.

Note that:

The larger the distance $d$, the smaller the correlation.
The correlation by itself does not tell you anything about the distance, since you also need to know $n$. For example, $1,1,1$ and $-1,-1,-1$ have correlation $-3$ but $1,1,1,1,1,1$ and $1,1,1,-1,-1,-1$ have correlation $0$.


Why the requirement of the GCD of the lengths of all circuits in the graph being one?

This is my understanding: the statistics of the source described in the paper depend on which character is produced first. If the first character is $a$, then one of the source properties is that letters in odd positions are always $a$. However, if the first letter is $b$ (in other words, a shift of $1$ in the circuit), then letters in odd positions can be either $b$ or $c$. This is what makes the source non-ergodic.

Note that the statistics when choosing $b$ as first letter and $c$ as first letter are exactly the same. In other words, there are $d=2$ different classes of sequences, and those in each class have the same statistics.

Regarding the importance of these sources: for one, Shannon's theorems depend on the source being ergodic. From a more practical standpoint (and this is purely my interpretation), non-ergodic sources are not good models for real-world processes; referring back to the example, no practical language produces repeated letters in a fixed periodic pattern.

I recommend reading "An Introduction to Information Theory: Symbols, Signals and Noise" by John R. Pierce, along with Shannon's paper. It explains many of the concepts in the paper in a very didactic manner. For example, it takes four pages to very clearly explain ergoidc sources, with nice examples.


Why FFT does not retrieve original amplitude when increasing signal length

Your code is bit unclear, especially generation of your signal. Python allows for vectorized operations so it is good to use it. What's more, it is good to clearly specify the sampling frequency of your signal and use it then. Also please remember to normalize your FFT by length of your signal (in this particular case) and multiply by 2 (half of spectrum is removed so energy must be preserved).

Phenomena you are facing is connected with fact that your frequency is not matching exactly frequency bin for it. Because of that energy is leaking into other frequency bins. This is due to fact that you don't have integer amount of cycles of your sinusoid. Please search for spectral leakage. Easy way to check that is to use the equation for frequency spacing in frequency domain, that is:

$$\Delta f=\dfrac{f_s}{N} $$

Where $f_s $ is your sampling frequency. Then your frequency vector is:

$$f=[0, \; \Delta f, \; 2\Delta f, \; ...  \;,(N-2)\Delta f, \;(N-1)\Delta f] $$

If frequency you are analysing matches exactly one in your frequency vector, then amplitude will be correctly estimated, otherwise - leakage. In order to minimize this effect, you might want to apply window to your signal. It will reduce amount of energy leaking to other bins and in slight way will improve estimation of original amplitude, although you won't be really able to retrieve it in 100%.

Here is small modification of your code to help you visualize that. The red curve is ideal spectrum and blue one is for given number of frequency bins (strictly connected to your signal length).

import numpy as np
import matplotlib.pyplot as plt
import time

def get_fft(y, fs):
    """ Get the FFT of a given signal and corresponding frequency bins.

    Parameters:
        y  - signal
        fs - sampling frequency
    Returns:
        (mag, freq) - tuple of spectrum magitude and corresponding frequencies
    """
    n  = len(y)      # Get the signal length
    dt = 1/float(fs) # Get time resolution

    fft_output = np.fft.rfft(y)     # Perform real fft
    rfreqs = np.fft.rfftfreq(n, dt) # Calculatel frequency bins
    fft_mag = np.abs(fft_output)    # Take only magnitude of spectrum

    # Normalize the amplitude by number of bins and multiply by 2
    # because we removed second half of spectrum above the Nyqist frequency
    # and energy must be preserved
    fft_mag = fft_mag * 2 / n

    return np.array(fft_mag), np.array(rfreqs)

def build_signal(A, f, N, fs):
    """ Generate cosinusoidal signal of a given frequency and amplitude.

    Parameters:
        A   - amplitude of a signal
        f   - fundamental frequency
        N   - number of signal samples
        fs  - sampling frequency
    Returns:
        t   - time vector
        y   - generated signal
    """

    dt = 1/float(fs)            # Get the time domain spacing
    T  = N*dt               # Upper time limit
    t  = np.arange(0, T, dt)    # Construct time vector

    y = A * np.cos(2*np.pi*f*t) # Generate the signal

    return t, y

if __name__ == "__main__":

    # Define the parameters
    fs = 100   # Sampling frequency (must satisfy the Nyquist theorem)
    f  = 10    # Signal fundamental frequency
    A  = 0.5   # Signal amplitude

    # Prepare figure for plotting
    plt.figure()
    plt.ion()
    plt.show(block=False)

    f_off = 1                     # Frequency offset in Hz for zoom

    # Precomute the spectrum without any leakage
    (t, y) = build_signal(A, f, 1000, fs)
    (MAG, FREQ) = get_fft(y, fs)

    for N in np.arange(1000, 10010, 1):

        # Clear and scale
        plt.clf()

        # Generate the signal
        (t, y) = build_signal(A, f, N, fs)

        # Get the spectrum and frequencies
        (mag, freq) = get_fft(y, fs)

        # Plot the reference spectrum
        plt.plot(FREQ, MAG, '-or')
        plt.hold('on')
        # Plot the spectrum for given number of bins
        plt.plot(freq, mag, '-o')
        plt.title('Frequency spectrum for %d bins' % N)
        plt.xlabel('Frequency')
        plt.ylabel('Amplitude')
        plt.grid('on')
        plt.hold('off')
        plt.legend(['Ideal spectrum', 'Spectrum for %d bins' % N])
        plt.xlim((f-f_off, f+f_off))
        plt.ylim((0, A))

        plt.draw()
        plt.pause(0.5)

Two intermediate steps. Signal is 10 Hz sinusoid with an amplitude of 0.5, sampling frequency is 100 Hz.

Exact match:

Frequency component you are analysis falls between two frequency bins:


Line coding for a modulation order that is not power of 2

Another reason to encode to $3$ level is that the occupied bandwidth is less.   This is often driven by an emission mask constraint  for transmitter certification. To satisfy this constraint,  the tradeoff is made on symbol distance.   It is however a costly tradeoff affecting receiver sensitivity because the multi level scheme is done post detection where every factor of $2$ increase in throughput costs $6$ dB in sensitivity.   Other schemes like QAM cost only $3$ dB in sensitivity when doubling the number of bits per signaling element.  With the multilevel FSK, you do however retain some of the advantages of an FM strategy in  terms of the CIR performance of the system.


I'll try to guide your research a bit rather than giving you a direct answer. One-click recipes aren't too good for obtaining optimal results, and this requires some careful insight. After all, you are dealing with biological signals.

To develop a real-time filter you need to process each sample/window before the next one comes up. You need to take into account the CPU time available for your calculations (depending on the platform) and the number of signals you must filter at once (with the subsequent impact in used memory, of course). After all, a digital filter is no more than a bunch of sums and products.

Now you should study the kind of signal you need to filter. You seem to be handling frequencies from delta to beta waves. You need an insight in the characteristics of each signal so that you keep them after the filtering (e.g. if studying brain activity during sleep, you should make sure that you don't miss out k-complexes and spindles during Phase II). This way you can look into some values: the maximum band pass ripple you are allowed, the width of the transition band, etc.

Having this information, after some research (e.g. Digital filter design by Burrus & Parks; Discrete-time signal processing by Oppenheim & Schafer) you can choose the kind of filter you need, and go into further details about implementation.

Fortunately, with a simple google query you can find out lots of things. There are Matlab toolboxes available to aid your filter design:

FDATool, Matlab filter design toolbox.
you can also take a look on ERPLab, which seems a good starting point.

And lots of free papers: Google query for "EEG filtering".

The higher the order of the filter, the more accurate it is, but the more operations per second it needs. So the main constraint here is that you need to keep it real-time for a bunch of signals.


How do you map indices to twiddle factors in FFT

Here is C++ code snippet that shows a relatively efficient way of doing it. The idea here is that the step size is always half the previous step and that the index wraps around if it exceeds N/2 which we can implement very efficiently with a bit masking operation (if it's a power of two).
constant int N = 8; // FFT size
constant int moduloMask = (N/2)-1;  // implements modulo N/2 adding

// loop over all stages
int stepSize = N;  // initial step size if FFT length
for (int i = 0; i < numStages; i++)
{
  int twiddleIndex = 0;
  stepSize >>= 1;   // divide by two, it's half of the previous stage
  for (int b = 0; b < numButterflies)
  {
    currentTwiddle = W[twiddleIndex];  // current twiddle factor from table
    // update the index: add step size and warp it around if it exceeds N/2
    twiddleIndex= (twiddleIndex+ stepSize) & moduloMask ;
  }
}


Are there any good open-source (well, free) handwriting OCR programs?

tesseract is probbaly the best and most widespread OCR library.

It has been tested with handwriting and isn't too bad - although handwriting isn't an easy thing to read. http://arxiv.org/ftp/arxiv/papers/1003/1003.5893.pdf


Symbol Synchronization for GFSK/FSK signals in GNU Radio

The Symbol Synchronizer block is a PLL-based synchronizer that is trying to estimate the symbol clock period and symbol clock phase (aka timing offset) based on the samples coming in that represent the data symbols.

Being a PLL configured with static parameters, there is a fundamental trade off between acquisition speed and tracking stability of the symbol clock estimate.  The more reactive the loop configuration, the faster it will acquire, but the more likely it may loose track later and cause an error.  The more sluggish the loop configuration, the longer it will take to acquire, but it will track much better when locked.

You can speed acquisition by using a correlator out front (e.g. the corr_est block) that detects your preamble and marks the center of a symbol with a "time_est" tag.  The symbol sync block will notice the "time_est" tag (or a "clock_est" tag) and reset it's tracking estimates right there to rapidly acquire lock.

But that is not your immediate problem.

For the symbol sync block to work well, you must understand the requirements, assumptions, and behaviors of the Timing Error Detector you have selected.

The M&M TED, with the default "bpsk" decision constellation, is expecting:

symbols that are peaked in the center
symbols that ideally have peak(trough) levels of $\pm 1$
you to set an expected TED gain that comes from modeling and simulation of the TED in this particular situation (expected SNR, samples symbol in use, matched filter in use, etc.) in MatLab, Octave, Python or other M&S tool of your choice.

When you model and simulate the M&M TED, you may notice something peculiar about its behavior:

If you have a good estimate of the symbol clock period, then the M&M TED S-curve, and hence TED gain, collapses to $0$ on sequences of alternating symbols: $..., -1, 1, -1, 1, -1, 1, ...$

This is likely why you are failing to acquire lock at the beginning of some packets when starting with some random symbol clock phase (aka timing offset) estimate.

Either use a correlator that marks a preamble symbol to jump-start the timing synchronization, or use a different TED.  The M&M TED is old technology.

When using a different TED, make sure you understand its derivation and the assumptions that went in to its derivation.

Also design and modeling and simulation in MatLab, Octave, R, Python, or whatever, is absolutely essential if you want something other than suboptimal results.  GNURadio is a fine prototype tool, but is not a very good up front design tool.

Please read this brief, if you haven't already:

https://www.gnuradio.org/grcon/grcon17/presentations/symbol_clock_recovery_and_improved_symbol_synchronization_blocks/Andy-Walls-Samples-to-Digital-Symbols.pdf


What is the most efficient way to perform Gaussian blurs for varying sizes of Gaussian kernels in an image?

I'm not really used to thinking in texture fetches as I usually code high level non-hardware-specific signal processing, but low pass filtering also removes information you can couple it with sub-sampling which can give further gains. It is not clear to me if you have already done this.

In practice this means that after the first low-pass filtering, you can disregard $k$ out of $N$ of the samples in the subsequent low-pass filtering.

Simple example is the $2^n \times 2^n$ averaging-filter, first we can do 2x2 mean values and store. Then for the next pass, we only need to use the output of one of each $2 \times 2$ blocks. So we only need consider $1$ in $4^k$ data points at step $k$.

A 8x8 block this would mean:

$$I_1 = \left[\begin{array}{cccccccc}
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1
\end{array}\right]$$
First iteration, all are involved in 2x2 filtering. But they are only stored at the pixels which are ones here (so we only need to calculate 1/4 of them):
$$I_2 = \left[\begin{array}{cccccccc}
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0
\end{array}\right]$$
Next step only involves the ones in the previous one, now our image is 4x4 instead of 8x8, so calculations are reduced a factor of 4. But we only need to store the results at pixels:
$$I_3 = \left[\begin{array}{cccccccc}
1&0&0&0&1&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
1&0&0&0&1&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0
\end{array}\right]$$
We will end up with only one out-value, in this case upper left. And here's where the interpolations come into play. We can interpolate between these upper left results of neighboring "blocks". The interpolation kernel can be dependent on the low-pass filter, but I assume that this is where the hardware shines with it's strengths, having really fast ways to do this for us?

If you are interested in theory filter-banks, sub-band coding and maybe wavelets are interesting to read up on. The example I brought up here is the implementation low-pass filter of the basic but popular discrete Haar-wavelet. Gaussian pyramids in computer vision may be even more relevant. There's much to read about it so it's difficult to cram it all into just one answer.


Concepts / Algorithms for recognising words in the Time-Domain

I highly doubt that something as complex as recognizing a word can be done directly in the time domain using features as basic as zero-crossing rate.


Pose from Homography matrix - Zhang method

The coordinate system you choose is completely arbitrary, as no information about real-world coordinates can be inferred. From an image of a table there is no reason to know that one leg is located at any particular $(X, Y, Z)$, or that it is any particular size (you can't tell if it's a doll's table or a giant's table).

Normally you would choose one of your cameras to be located at the origin, looking down the $z$ axis, defined by the matrix:

$$[R|t] = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1&0 \end{bmatrix}$$

Then due to the scale ambiguity you would have to choose an arbitrary scale, for example if you are using a stereo camera you could set the distance between the cameras to be one unit distance.


What are the Segmentation algorithms that can be used to process Satellite images (obtained from Google Earth) to segment forest and barren areas?

You could try image saliency or entropy filtering. Entropy filtering calculates the entropy of gray values in the image Image within a rectangular mask. The result is written to the output image. The process is repeated in a sliding window fashion.

In the resulting image, the forest areas should appear more whitish while the barren areas appear black, because they don't contain enough information.

Of course I assume the roads would make a problem. However you could get rid of them using morphological edge finding operations.


How to find fundamental frequency of Fourier Series

Your solution is correct. Obviously, any periodic function with period $T$ is also periodic with period $kT$, $k\in\mathbb{Z}^+$. In the given example, the smallest possible period is $T=1/5$, which is equivalent to $\omega_0=10\pi$. There is no obvious reason why one should choose $T^\prime=kT$ and $\omega_0^{\prime}=\omega_0/k$, with $k>1$. Nevertheless, such a choice, even though unmotivated, is equally correct, because the resulting Fourier series are identical:

$$\begin{align}x(t)&=\sum_nc_ne^{jn\omega_0t}\\&=\sum_nc_n^{\prime}e^{jn\omega_0^{\prime}t}\\&=\sum_nc_n^{\prime}e^{j\frac{n}{k}\omega_0t}\\&=\sum_mc_{km}^{\prime}e^{jm\omega_0t}\end{align}\tag{1}$$

With $\omega_0^{\prime}=\omega_0/k$ we get

$$c^{\prime}_n=\begin{cases}c_{n/k},&n=mk\\0,&\textrm{otherwise}\end{cases}$$


differences between two L1 norm minization schemes

The first equation you have is often called the Quadratic Problem, which through the use of Duality can be shown to be equivalent to the Basis Pursuit De-Noising (BPDN) given as:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{1}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}< \delta .$$
In your 2nd problem, the L2 norm is replaced with the L1 norm in the constraint. When you constrain the L2 error, you tend to end up with an error vector that is quite dense, i.e. a lot of non-zero elements, and each error is small with respect to $\delta$.
By switching the error to the L1 norm, then they are trying to make the error vector sparse (as opposed to dense). So in this formulation:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{1}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{1}< \delta $$
They are trying to find:

An $\boldsymbol{x}$ that is sparse under the constraint that
$\boldsymbol{e}=A \boldsymbol{x} - \boldsymbol{y}$ is also sparse.

In this sense you might think that the problem they really want to solve is:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{0}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{0}< \delta $$
but because that is a very difficult problem, they are using the convex formulation they have given.


Is it possible to just remove the noisy part of signal?

The short answer is no, unless you know exactly what the noise is, you can not in any practical sense remove the noise. There are things you can do to maximize SNR but no technique exists to remove all noise.


Is there a way to obtain the transfer function from a bode plot on Python? (I know that it is possible on Matlab)

Yes, use the control toolbox that is available in python. You can install it using pip install control

From here you can implement transfer functions, bode plots, step responses etc similar to what you would be accustomed to doing with the Matlab control toolbox.

Here are further details:

https://github.com/python-control/python-control


What is the relation between windowing and hopping in audio DSP

What many people do not seen to realize is that using a finite length transform (such as one iteration of an STFT FFT) on a longer signal already alters the longer signal by windowing (e.g. ignoring the portions of the signal that do not fit in the current finite length FFT).  Not using a tapered window means that one is using a rectangular window.  Using a rectangular window may alter the signal spectrum even more than using a tapered window (from, say the infinite length FT of an infinitely long stationary signal, especially if the signal is not perfectly integer periodic in the STFT FFT aperture width).

The FFT spectrum is altered the least near the center of the window (where a tapered window has a gain of unity or a maximum gain, and non-tapered signal is farthest from both potential edge transients).  Therefore the closer together the window centers of individual STFT FFTs (the shorter the hops), the less the alteration or distortion of the total STFT spectrum.

Lots of overlap also provides greater time localization or locality (e.g. the STFT won't miss or slice up a short event that is right between two windows).
 However the greater the overlap, the greater the computation required.  So a compromise overlap of 50% or 75% window overlap is common.


Depending on the context, the use of the complex form could be for mathematical convenience or for a no-kidding need for both real and imaginary parts.
When you factor the expression, you get
$$u(t) = e^{{\eta}t^2}e^{j{\beta}t^2}$$
Where the first exponential is a generic magnitude envelope, in this case Gaussian. The second exponential is the chirp itself and is where all the action is.
So to simplify things, lets assume that the envelope is ideal so that the signal is just the chirp
$$u(t) = e^{j{\beta}t^2}$$
To view the chirp, you can take either the real or imaginary part, and it looks something like this

In a system that does not use I/Q, the real part is what you would expect to transmit as your waveform. The type of system will determine whether using a real signal or some type of I/Q is best. I'm going to use a radar example here.
In frequency-modulated continuous wave (FMCW) radars, a real chirp like the one above is transmitted and has the form
$$x(t) = cos({{\beta}t^2})$$
Which is just the real part of the complex form. It is received after a delay and mixed with itself, and without going into the mixing process, produces a single frequency sinusoid that can be used to determine range. Here, using a real part only is practical. Using I/Q in FMCW is also beneficial (SNR improvement), but is not usually necessary and many systems do not use it.
Another type of radar, pulse-Doppler, benefits greatly from using the complex form. The same chirp is considered, except now the imaginary version is used. This is important because pulsed-Doppler radars usually operate on performing pulse compression, which is just correlating the transmitted waveform with the received one.
The autocorrelation of a complex chirp looks like

Using a complex waveform allows us to mix our signals to baseband, which give the classic autocorrelation responses we expect without additional mixing and filtering.


what is the reason for the Change in DFT analysis after downsampling

Because you have resampled the signal to a sample rate below twice the highest frequency in the signal, you are seeing the signal plus an aliased image of the signal both summed together.  Two partially overlapped rectangles will sum to produce the "head-and-shoulders" shape you see in the aliased spectrum.


How to gather similar objects in an image?

First, I personally wouldn't recommend you using background subtraction in general for tracking applications.
There are way better methods like L1-minimzation of lasso problems, Kenalized Correlation Filters, Multiple Instance learnng, Struck, and whole bunch of stuff.

But if you must use background subtraction,, here are my thoughts:

You could obviously at the first frame apply some very basic background subtraction algorithm like using median filtering.
Now once you get these blobs, you can start generating a dictionary out of each blob.
That means you can now treat each blob as an individually part where you learn a dictionary for each. Usually in dictionary learning from first frame, people usually perturb the bounding box covering the blob a little (slight translations)(Maybe by gaussian sampled translations with zero mean) and then collect all these samples and vectorize them into your dictionary matrix.
You do this for all the other parts, or blobs if you will.

Now at the first frame you already have the bounding box governing the object ,which is very typical in any tracking problem. Also, you have the bounding box governing all these blobs individually.

What you do next in the next frame is to try track each of these parts independently say using simple cross correlation.
At last, you only need to learn one affine transformation that transforms the big box from the previous frame to cover all the blobs in the next frame.

In this approach you divided your main task of tracking from 1 bounding box to smaller ones on which you can use simple methods of tracking. Then you combine 1 bounding box by either learning an affine transformation or simply plotting a minimum bounding box that covers all the other parts.


How to reconstruct a time-discrete signal

Consider things you know about the sampling theorem:

If it is met, the original signal can always be reconstructed perfectly!
The signal has to be bandlimited for the theorem to be met!

Further, think about the spectral and time domains and their relation to each other:

If a signal is bandlimited in spectral domain, it has to be infinite in time domain.

Then, thinking about the given signal:

The given sampled signal has just zeros before and after the one none-zero value.
It is not periodic.
One can conclude, that in the original signal, there have to be none zero values between all (or at least some) sampling points, otherwise it would not be infinite.

So we have to ask ourselves the question:

Which signal is bandlimited, has zero crossings at regular intervals, a non zero value in the middle and is not periodic?

If you had any lectures on signal theory, the answer should be pretty obvious.


Nyquist plot interpretation when curve hits the origin

First to clear up the OP's misunderstanding: the Nyquist Stability Criteria involves clockwise encirclements of -1, not the origin, and this would be the polar plot for the open-loop gain specifically.

I've included some details below for those that are more interested.

First a review of the basic equation relating Open Loop gain and Closed Loop Gain for a control loop:

$$G_{CL}= \frac{G_F}{1+G_{OL}}$$

The denominator of the closed loop gain is known as the "Characteristic Equation".  Given that all physical processes that are linear time-invariant have transfer functions that are proper (the degree of the numerator cannot exceed the degree of the denominator), we are able to determine stability from the roots of the characteristic equation alone. The roots (zeros) of the characteristic equation are the poles in the closed loop transfer function (since the characteristic equation is the denominator of the closed loop transfer function).

Therefore all stability assessments are methods to determine if any of the roots of the characteristic equation are in the right half plane (and therefore unstable).

What Nyquist did in creating the "Nyquist Stability Criteria" was make use of Cauchy's Argument Principle as a method to determine this. The significance of his approach was that instead of using the characteristic equation, he was able to do this with just the Open Loop Gain alone ($G_{OL}$). And what you will observe when you see the approach is that you do not actually need an equation for the open loop gain; we can go in the lab (or simulate a complex digital design) and measure the transfer function as a magnitude and phase versus frequency, and from that measurement directly determine stability without having to solve for poles and zeros and polynomials etc.  This is very significant so I will stress that again with examples: We could have a system that has time delays that would not easily translate to a polynomial, or a digital system with decimation resulting in fractional powers of z, or a complex system where the effort to measure the magnitude and phase versus frequency would be a fraction of the effort to create an actual equation with accuracy. Often digital systems will have parasitic $z^{-1}$ delays that may easily be missed in the effort to create a transfer function, but cannot hide from the visible phase delay that can be measured.

Cauchy's Argument Principle

The essence of Cauchy's Argument Principle as it applies to Nyquist's method is that if we map an arbitrary contour that encircles poles and zeros of a function H(s) on the complex plane, the mapping will encircle the origin in the same direction Z-P times, where Z is the number of zeros and P is the number of poles.  This is demonstrated in the examples below showing the effect of each zero and pole; For each zero inside the contour, the mapping will encircle the origin in the same direction, for each pole inside the contour, the mapping with encircle the origin in the opposite direction. Anything outside does not contribute to an encirclement. An encirclement means the phase completes 0 to $2\pi$, for multiple poles and zeros, their individual phase (argument) contributions add and hence the name "Cauchy's Argument Principle":

Nyquist Plot

So using Cauchy's Argument Principle, Nyquist made the contour the entire right half plane to determine if any roots of the characteristic equation are inside the contour (and therefore an unstable system). Now if the Nyquist method did the contour for the characteristic equation itself, the encirclements would indeed be around the origin N= Z-P times as described above in the development of Cauchy's Argument Principle. Since the characteristic equation is $1+G_{OL}(s)$ and the roots are the solution to $1+G_{OL}(s)=0$, the Nyquist method instead does the contour around only $G_{OL}(s)$ and from $G_{OL}(s)= -1$ the encirclements if they exist due to roots in the right half plane will be around -1 (By subtracting 1 from both sides, we have simply shifted the origin but more importantly isolated the function for the open loop gain alone, as mentioned above it is a system parameter we can easily measure.

So, if there are clockwise encirclements of -1, the system is unstable as the number of clockwise encirclements is N=Z-P, and P cannot be negative. However from this we also see that if there are no clockwise encirclements, that does not guarantee stability unless we know that there are no poles in the right half plane for the open loop gain (P = 0). If our open loop system is stable, then we know that is true and we can make conclusions with regards to closed loop stability. If not, then in that condition we would need to determine the number of poles and their location (right or left half plane) in the open loop gain before concluding stability.

NOTE: I described the Nyquist plot in units of s implying an analog system, but everything described applies to digital systems as well in units of z (it is just the frequency magnitude and phase is done by sweeping the unit circle instead of the $j\omega$ axis.) It is interesting to observe the example below of an analog system that I mapped to digital using the method of impulse invariance, showing that the Nyquist plot for both systems is identical (at this scale). Point is the result of counting clockwise encirclements of -1 is still the same with the digital system:


Shortest geometric distance from surface in 3d dataset?

Do you just need the distance, or do you need the closest point?

For the closest point, the FLANN library can help, and it has Matlab bindings.

If you only need the distance, you can also use a distance transform. Try googling for "distance transform 3d matlab" for implementations.

Which one is faster depends on the number of "seeds" and "skin voxels".


(this should possibly be a comment but I don't have the rep yet, so I'm submitting it as an answer) Restricting your view of the signal to a finite time window will always introduce artifacts, unless the signal itself is zero everywhere outside of that window. Using a "square" window (simple sampling the signal for a finite duration) will have the effect of convolving the signal's true spectrum with a Dirichlet Sinc function (Dirichlet Kernel). See https://en.wikipedia.org/wiki/Window_function#Spectral_analysis


Is there an equation that will graph the shape or predict the maximum amplitude of a single sample impulse fed through a One-Pole LPF?

As mentioned in the comments, an impulse response with the shape shown in your question can only be obtained by a system with two real-valued poles. In continuous time, with two distinct poles, the total impulse response is
$$h(t)=\frac{e^{-\alpha_1t}-e^{-\alpha_2t}}{\alpha_2-\alpha_1}u(t)\tag{1}$$
It's a basic exercise to determine the location of the maximum of $(1)$, and the result is
$$t_{max}=\frac{\log\left(\frac{\alpha_1}{\alpha_2}\right)}{\alpha_1-\alpha_2},\qquad\alpha_1\neq\alpha_2\tag{2}$$
For $\alpha_1=\alpha_2=\alpha$, the total impulse response is
$$h(t)=te^{-\alpha t}u(t)\tag{3}$$
and the location of the maximum is
$$t_{max}=\frac{1}{\alpha}\tag{4}$$
The figure below shows an example of an impulse response $(1)$ with $\alpha_1=0.4$ and $\alpha_2=0.6$, and with the location of the maximum predicted by $(2)$.


How to choose basis functions that contribute most efficiently per term to an approximation of an image f(x,y)?

I do not  understand why you need an infinite number of  terms for a Fourier approxiimation. Yet, if you are looking for tools on non-discretized 2D data,  I would suggest you to start with the concept of Matching pursuit of images. You take a dictionary of 2D atoms that is complete (could be complex exponentials, gaussian bumps). Then, you find the  atom that best matches you  data, subtract it with the necessary scale factor, and repeat.

You could also resort to anisotropic  triangulations, but  this would not suit your "need to be combined in an inner product fashion".

The  quality of approximation could be evaluated over some classes for regular functions ($C^p$ smooth inside $C^q$ curves), and are related to properties of  your  dictionary. You might find some references in section 4.1.1. Matching pursuits from A panorama on multiscale geometric representations, Signal Processing, 2011, L. Jacques et al.

I do suspect that the inner-product  constraint will limit the efficiency: optimal approximations in 1D do not extend generically  in 2D.


Created Matlab echo signal: Is it an actual "echo?" Or am I just manipulating the signal to make it seem like an echo?

To generate a repetitive echo all you need to do is pass your signal through an appropriate AR filter with a delay corresponding to the echo period. For instance the following code generates (for a damped 1 Hz sinusoid sampled at 100 Hz) an echo which repeats every 1000 samples. The strength of the echo can be adjusted by the AR filter coefficient (0.5 here)

Example Matlab code

% generate 10000 samples of a damped sinusoid at 100 Hz
sig = sin(2*pi*(0:9999)/100).*exp(-(0:9999)/200);

% AR filter to generate a repetitive echo, period 1000 samples (10 seconds)
b = 1;
a= [1 zeros(1,999), 0.5];

echo = filter(b,a,sig);

You may want to experiment with your AR filter to match your echo decay. Also look into 'gapped deconvolution'


segmenting region of interest (ROI) from background

Some very cool although a bit more time consuming options:

Active contours will try to find the best "edge" for the cells while maintaining smoothness.
and is implemented in matlab2013a ( active contour)
You might try using mser which is an approach which detects regions which are similiar and in color and stable over the image. which means they stay about the same size over a wide array of thresholds.
I used it for a similar problem with detecting ultrasound tumors.
If you have a good database of images ( good = covering most cases )  I would try some machine learning algorithms. for example, you could learn on the mean and variance of the objects in the image.
If you have trouble separating the cells after segmentation ( i.e. you have two connected cells ) I would try looking at the watershed function.


How to detect all the connected neighboring pixels which have same color value in opencv

I think you're looking for floodFill. From the documentation:

The functions floodFill fill a connected component starting from the seed point with the specified color. The connectivity is determined by the color/brightness closeness of the neighbor pixels. The pixel at (x,y) is considered to belong to the repainted domain if: [...]
Use these functions to either mark a connected component with the specified color in-place, or build a mask and then extract the contour, or copy the region to another image, and so on.

You're probably interested in this flag parameter:

FLOODFILL_MASK_ONLY If set, the function does not change the image ( newVal is ignored), and only fills the mask with the value specified in bits 8-16 of flags as described above. This option only make sense in function variants that have the mask parameter.


Given a spectrogram, how do I calculate loudness over time?

Yes, in fact it's not a bad idea at all to use some kind of spectrogram representation to calculate a loudness envelope. Depending on what exactly you need you can either just sum the spectral magnitude over all frequencies at a given time, or first apply some kind of perception based weighting functions
to make your loudness measure more related to actual auditory perception.

So if you give some more detail about the application you have in mind you will likely receive a more accurate answer.


Difference between jpg compression and mp3 compression

Actually, it's kind of the other way around. If you reuse the same JPEG encoder at the same quality level (without any smoothing steps as built-in prepcosessing) and a decoder which faithfully decompresses the images, I expect the image quality not to degrade from generation to generation. This is because quantization (the lossy part) is done the same way during the 2nd compression. You can imagine scalar quantization like this: given some value v and a quantizer step s, the quantized value is round(v/s)*s. Example:

v = 6.3
s = 2
q = round(v/s)*s = 6

So, the first "compression" includes an error of 0.3 in this case. But if you feed this into the encoder a 2nd time

v = 6
s = 2
q = round(v/s)*s = 6

it does not change anymore. It stays 6. If you change the quality level from one compression to the next, you also change the quantizer step s. This would add a new rounding error.

Keep in mind that this is just a simple example. JPEG Quantization may include a "dead zone quantizer" which changes the threshold a bit near 0. For example, you might want to round any value v to zero in case |v|<0.7. Something like this might give you a slightly higher quality-per-bit ratio because zeros can be much more compactly represented. But this should not change what I described above.

With MP3 the "quantization steps" are computed dynamically and depend on the actual signal. This is what the "psychoacoustic model" is for. It analyzes the audio data and tells us how strongly we can quantize certain portions of the signal without us noticing a difference. But this quantization adds "quantization noise" which could make the psychoacoustic model compute other quanzitation steps during the 2nd compression even in the case where you told the encoder to use the same "quality level". In the case the 2nd compression uses different quantization steps, you'll introduce new rounding errors.

Actually, something like this has been done before. Some artist did such an experiment where he compressed some audio file over and over again (like 1000 times) and the result had little to do with the original. Unfortunately, I don't remember any names or references.

There are other issues which make lossy audio compression more complicated. For example, a decoder might add or remove some leading samples (encoder+decoder delay mismatch) which would make the encoder work on different blocks for the 2nd pass. This would also add new rounding errors. But if I remember correctly, such a "block position mismatch" was shown to be beneficial w.r.t. generation loss by some experiments done by the Hydrogenaudio crowd a couple of years ago. So, while there was still a generation loss, the loss was less annoying than with matching block positions.


How to merge 2 WAV files together? Can I take the average of the 2 audio file float values?

If your goal is to simply make it sound like they were playing at the same time, you can simply just add them together element-wise. There is need to "average". The resulting waveform is the same, but scaled by a factor of 2 in your case. Technically, that number could really be anything.
If you want different properties, where maybe one sounds louder than the other, then you'll have to do something else.


Extracting DCT coefficients from H.264/AVC bitstream

It is possible to extract the transform coefficients without fully decoding the entire picture, but you would still have to find your way through the bit stream, parsing whatever is necessary in order to arrive at every transform unit, which might also not be present for every coding unit (if some macroblock matches exactly one in a reference picture then you don't need to encode any residual, for example).
How you would use only the transform information to do your computer vision tasks, I have no idea, but it is entirely possible to only extract the transform coefficients. What i suggest you do is grab whatever standard you want (H.264 or H.265 or some other) and go through all of the tables with the syntax elements, and the parsing process for arriving at the transform coefficients would be to ignore every fixed width syntax element not useful for decoding the transform (anything related to motion vectors or reference pictures) and parse whatever is necesary to obtain all the transform units in the coded picture. I don't know if there is any tool out there that does this, perhaps you can go into ffmpeg or one of the libraries like x264 or x265 and modify them to extract what you want. It is a difficult task, but definitely possible.

Hope it helps


What's the point of MATLAB's scaling array for secord-order-section filters?

Of course an SOS matrix is capable of completely describing the filter including the gain. However, in the Matlab SOS matrix the $b_0$ coefficients are all equal to $1$ by default (if no extra scaling is used), which is a kind of "neutral" scaling (and gets rid of one multiplication per second-order section). The gain factor $G$ is the overall gain necessary to describe the filter in addition to that scaled SOS matrix.

Note that for fixed-point implementations each second-order section should be scaled individually. This can be done with different goals in mind, e.g., to minimize the probability of overflow, or to minimize the quantization noise.


I am wondering whether it is possible to express $X_{2F_s}$ in terms of $X_{F_s}$.

In general, no.

This is the sampling theorem: unless you've bandlimited $X_{2F_s}$ first to be alias-free sampled by $F_s$, then the digital signal obtained by sampling twice as fast will contain information that the digital signal sampled at $F_s$ does not contain. And information that is not contained cannot be recovered.

If you, however, did not band limit $X$ to the Nyquist bandwidth dictated by the sampling rate, your $X_{F_s}$ will already contain aliases. Then you can simply look at the "upper" half of $X_{2F_s}$, shift it down to overlay exactly the lower half, which is exactly what's called aliasing, and reconstruct $Y_{F_s}$ from $Y_{2F_s}$.


Turn acausal Butterworth output into output of a causal Butterworth without input signal

This is possible, as long as you know the filter which was applied to the data in the first place.

Assuming that the original acausal filter is $H_a(\omega)$ and your filtered data is $Y_a(\omega)$, then the original data $X(\omega)$ can be retrieved by inversely applying the filter to the data:
\begin{equation}
X(\omega) = \frac{Y_a(\omega)}{H_a(\omega)}
\end{equation}

This is simply a re-arrangement of the equation which can be used to describe the initial  filtering operation:
\begin{equation}
Y_a(\omega) = X(\omega) H_a(\omega)
\end{equation}

I've illustrated the concept using frequency domain notation, but the concept is the same for the time domain.  Depending upon your acausal filter, you may encounter some problems inverting it, I don't think this will be an issue in this case.


Compute the Second Order Derivative of Digital Image with Finite Differences

Pay attention that convolution mean flipping the kernel both on the x and y axis.
Hence the first element is a multiplication of $ -1 \cdot -1 $ which yields $ 1 $ as in the answer in the original post.
Pay attention that this is a discrete approximation of the gradient based on Finite Differences. This specific one is based on the Sobel 1st derivative filter.
In Finite Differences Coefficients you may find more variations.
As usual, the longer the filter the better the accuracy of the approximation with longer tarnsients.


What does a "vector" in a hidden Markov model mean?

The way speech recognition is carried out with HTK (or any other tool) is sort of similar to the way speech recognition is carried out in the brain. When you hear a word, you instantly break it down into its constituent phones and then compare the phones with an internal mental "model" of the phones. These "models" are constructed over years of listening to speech and gives you the ability to distinguish between similar-sounding sentences like "How to wreck a nice beach" and "How to recognize speech". Speech recognition with HTK or any other model-based scheme works in a similar manner. Here, in a few steps, is how you do it:

You take the input speech signal and convert it into a feature vector representation.
Take a large number of sentences and perform step 1 on each of them.
Use the feature vectors in step 2 to build a statistical model for each of the phones/words in the sentences (there are a limited number of phones/words as against an infinite number of ways of saying them - so you reduce the unknowns by modelling).
When a new word comes in, break it into phones and compare with each of the known models. The sequence of phones with the highest probability wins!

All the above steps are critical to the successful completion of any speech recognition task. By decomposing a sound into its feature vector, you are taking it into a model-space, giving it a representation that makes it more suited to making-a-model-out-of than other representations (say the time-amplitude representation). Most such representations lie in the frequency, or the time-frequency domain. One of the most popular such representations is the MFCC (Mel Frequency Cepstral Coefficient). In a way, this technique mimics the human hearing response with a set of filters. An input signal is decomposed with this set of filters that have a logarithmic spacing of their center frequencies. The MFCC coefficients of any one sentence (say) are then used to model each of the phones that the sentence is made of. As an example, consider,

Sentence: HI.   Phonetic description: hh aa ey
When you feed the MFCC coefficients into HTK, it will associate the MFCC coefficients of a portion of the sentence with hh, another with aa and so on. When this is repeated many times over, the models for the phones begin to form.

HTK uses the tool HCopy to convert an input sentence into its feature-vector representation. There are many "flavours" to MFCC's as well (E_D_A or E_D_A_Z representations). It would be a good idea to read up on the documentation for HCopy within the htkbook.

The MFCC coefficients are written to a file with extension .mfc by HTK. It is not possible to read that file using any one of the text editors because (I think) the coefficients are written in binary. You can try to read the files with C though.

HTH.