For a discrete LTI system, does "bounded memory" imply "rational transfer function?"

Non-trivially, the answer is no, a Bounded In Memory (BIM) LTI system do not imply the system has a rational Transfer Function (TF).
Proof (TL;DR):
Every discrete LTI system $H$ has an Impulse Response (by definition), and hence have a TF in the z-domain (by construction), which could be infinite (IIR) or finite (FIR).
A FIR LTI is trivially BIM, requiring a bounded number of operation for implement every convolution.
A IIR LTI could have a Rational TF (RTF) -such as most ARMA systems- or could have an not rational -an infinite series (IS)- TF -such as the counter-example 1 below.
If the system has an IS TF which can be decomposed as a rational TF -such as the Counter-Example 2-, then the system is bounded in memory, because the rational TF can be decomposed into a pair of finite series, which can be implemented and evaluated with finite memory and finite CPU.
Similarly, if the system has an IS TF with an Known algebraic expansion (KIS), -such as the Counter-Example 1-, then again the system is BIM, because, the general term is computed in each step, consuming an increasing amount of CPU -this is, an unbounded CPU system- but still BIM.
If this is not the case, that is, if the system has an IS TF without an unknown algebraic expansion (UIS) -which is indeed real for most LTIs- then is not BIM in practice.
The condition here for an IS to be not rational is having no poles and no zeros, that is, being bounded:
$$0<H(z)=\sum_0^\infty c_nz^{-1}<\infty$$ $|z|>z_0$ for some $z_0$.
If the series crosses by zero through $z_i$, or diverges for some $p_i$, then it has a rational set of component $(1-z_iz^{-1})$ or $(1-p_iz^{-1})$, which is nothing more than a polynomial that can be factorized, being this condition applicable in a recurrent way.
Hence our definitions:
ALL $\leftrightarrow$ FIR $\lor$ IIR
IIR $\leftrightarrow$ TF $\lor$ IS
IS $\leftrightarrow$ KIS $\lor$ UIS
And our known relations:
FIR $\rightarrow$ BIM (By construction)
TF $\rightarrow$ BIM (By construction)
KIS $\rightarrow$ BIM (Counter-Example 1)
UIS $\rightarrow$ not BIM (Infinite Terms)
And finally, a Bounded in Memory System is a FIR, a TF or a Known Expansion Infinite Series.
BIM $\leftrightarrow$ FIR $\land$ TF $\land$ KIS
Counter Example 1:
The LTI system, with IIR defined by its Z-transform, is BIM, but has not rational TF:
$$H(z)=log(1+az^{-1}), \\ h(t)=\frac{(-1)^{t+1}a^t}{t}\Gamma(t), t>0$$
We are using the exploit that:
$$log(1+x)=\sum_{n=1}^\infty \frac{(-1)^{n+1}x^n}{n}, x<1$$
which in general is not valid, thus this procedure will produce unbounded LTI systems except for the trivial power series functions $sin(x)$, $cos(x)$,$e^x$...
We are also assuming that our TF must be finite in terms, which is the purpose of having a TF, provided that all Impulse Responses can be converted into series in Z space.
http://www.comm.utoronto.ca/~dkundur/course_info/362/6_Kundur_Rationalz-Transf_handouts.pdf
Counter Example 2:
The IIR AR System has both a rational and an infinite series transfer function:
$$H(z)=\frac{1}{1-az^{-1}}=\sum_{n=0}^{\infty}a^nz^{-n}, |a|<|z|$$


Why is discrete cosine preferred to FFT in neuroimaging GLM

Say we have a signal,

The DFT assumes the signal is periodic, e.g.,

While the typical DCT assumes the signal is even at the boundaries, e.g.,

In the DFT example you can see there is a big jump from one side of the signal to the other. Since it assumes a periodic signal, this big jump introduces a lot of high frequency components. Thus the energy of the signal is spread out over more frequencies. For example, the log of the absolute values of the DFT is

Because the DCT assumes the signal is even at the boundaries, generally there is no big jump and thus no extra high frequency components. If this assumption is correct, more of the energy is concentrated in the lower frequencies. For example, the log of the absolute values of the DCT is

Actually, point at which the DCT assumes the signal is even can be moved. It doesn't have to be at the boundaries. However, if your signal generally has some large ramp component, then it makes sense to use the boundaries.

For the edge condition the DCT is equivalent (with multiplication by a scalar) to the DFT of the mirrored signal,

The log of the absolute values of the DFT of this signal is

which is the same as the DCT result above.

Note

The DCT doesn't automatically reduce boundary effects. Consider the extreme example of a sine signal, its DFT and DCT:

Because the even assumption didn't hold, the DCT actually resulted in many more high frequency components.

Why use it

I'm guessing because of drift in the nueroimaging signal, it is not periodic, and thus there would be a big jump from start to end. This would introduce lots of high frequency components with the DFT. Using the DCT assumes the signal is even and reduces the boundary effects.


Transfer function of a frequency shifting system

The system you're looking for cannot be described by a transfer function because it is a time-varying system. Only linear time-invariant (LTI) system can be fully characterized by a transfer function. However, there is no LTI system that can shift frequencies. The output of a (stable) LTI system can only have frequency components that are already present in the input signal.

Frequency shifts are usually achieved by modulation, i.e., by multiplying the input signal with a sinusoid or with a complex exponential. Such a system is linear but not time-invariant.


What is the relationship between the discrete time and continuous time variables?

It's about the relationship between the discrete time signal $x_d[n]$ and the continuous-time signal $x_c(t)$:

$$x_d[n]=x_c(nT)\tag{1}$$

So you formally replace the variable $t$ by $nT$ but this just means that you sample the continuous-time signal at sample instants $t_n=nT$. So $t=nT$ is only true for the values of $t$ that we're interested in, and these are the discrete values $t_n=nT$.


Minimizing profile curvature of 2D image

I am in search for (fast) algorithm to minimize profile curvature of this image on selected scales (in scale-space sense).

  ...
  ... to minimize curvature at specified scale, [...], trying to apply this accross all scales (I am trying to utilize Gaussian / Laplacian image pyramids). It works by adjusting heights by adding height offset according to profile curvature at given point (height offset is negative in convex areas and positive in concave areas).

Modifications of curvature are basically modifications of a "signal's" dynamics.

In effect, what you are trying to say is "If the signal is rising (or falling) faster (or slower) than some threshold then...do something different". This "do something different" can be anything from "don't respond" to "respond slower" or "respond faster" and so on.

In one dimension, the Douglas-Peucker algorithm does more or less what you are after.

The algorithm takes only one parameter which you can relate to the scale of observation but essentially modifies the local slope so that it remains within certain limits.

Variants that do surface simplification, or mesh simplification exist too (or this, in the case of contours) which use some expression of "cost" or "stress" to decide which element to remove during the simplification.

If you are "tied" to a height map (i.e. implicit 3D data), then you are looking at removing data (i.e. decimation) and interpolating / filtering. For example, you can use your $P_c$ to "mark" the areas of high or low curving slope in one pass and compute the "fill-in" values of each patch on a second pass. For examples of these approaches, you might want to take a look at this or this or more generally this.

Hope this helps.


Binary images are clearly smallest in terms of size, followed by indexed (paletted) images. RGB and greyscale images can take different sizes.

There are also (A)RGB images with sizes up to 48 bits but they usually include an extra alpha channel (for transparency purposes) so I haven't described them below.

RGB  (red, green, blue ) - 24 bits (each colour channel takes up 8 bits of space )

RGB (red, green, blue ) - 16 bits ( the green channel usually takes up 6 bits )

RGB (red, green, blue ) - 15 bits ( 5 bits per channel)

RGB (red, green, blue ) - 8 bits ( usually encoded as rrrgggbb )

Greyscale image - 8 - 24 bits ( each colour channel in a greyscale image has the same value)

Indexed image - usually 8 bits  ( the pixels don't carry colour information at all--rather they are indices into a colour table (i.e a palette). An 8-bit palette can hold 256 colours. While larger palettes are possible, they're rarely used because the most important reason to use a palette is to conserve space )

Binary image - 1 bit ( in a binary image, each pixel is stored as 1 bit. Since 1 bit can encode only 2 values (0 or 1), this kind of image can display only 2 colours (e.g black and white etc. )


How to treat noise (in the acceleration) in a Kalman filter when tracking position?

Uh, this is the whole point of a Kálmán filter.  It has a stochastic measurement model (noisy measurements) and a stochastic process state update model and uses the measurements for estimating the current state of the stochastic process.

If the process were deterministic, you would approach perfect information over time.

So you just plug the variation of your acceleration into the covariance matrix for the state update.

The exact way to do that depends on the conventions your course material employs (lots to choose from).


See "ePeriodicity: Mining Event Periodicity from Incomplete Observations" by Zhenhui Li, Jingjing Wang, and Jiawei Han (2013); preprint at https://faculty.ist.psu.edu/jessieli/Publications/tkde14.pdf


I am not really sure which part you didn't understand. I'm assuming it's the step that goes from the rightmost expression at the top to the one at the bottom.
Assuming $\mathbf{q} = \mathbf{x} - \hat{\mathbf{x}}$ is i.i.d.
Let me give it a try...
$
\begin{aligned}
\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) ]
& = \textrm{tr} \{ \textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) ] \} \\
& =  \textrm{E}[ \textrm{tr} \{(\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) \} ]  \\
& =  \textrm{E}[ \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}}  \} ]  \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}} ] \}   \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\textrm{E}[ (\mathbf{x} - \hat{\mathbf{x}}) (\mathbf{x} - \hat{\mathbf{x}})^{\mathrm{H}} ] \}   \\
& =   \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H} ~\sigma_{q}\mathbf{I}] \} \\
& =   \sigma_{q} \textrm{tr} \{\mathbf{H}^{\mathrm{H}}\mathbf{H}] \}.
\end{aligned}$
where I am assuming $\sigma_{q} = \textrm{E}[  \|(\mathbf{x} - \hat{\mathbf{x}}) \|^2 ] = \frac{\textrm{E}[  \|\mathbf{q} \|^2 ]}{N_t} $.
As you can see, I am making a lot of assumptions here. I don't know if they hold since I couldn't access the paper. But I think it can shed some light on your problem (if your question is strictly about the trace operations).


Block by block CCDF and PAPR analysis in MATLAB

PAPR of an OFDM-like system is typically calculated analytically: What's the the output signal with the highest possible peak? How high is that? Divide by average output power, done.

CCDF is just 1-CDF. Taking a wild guess, your OFDM-like system does a linear combination of independent, identical random variables (like in the DFT formula). The central limit theorem applies, and you get a Gaussian CDF for the amplitude with zero mean and variance equal to the average output power. Also not hard to calculate!

I'd argue that CCDF can relatively reliably be approximated by a monte carlo simulation (i.e. you just modulate a lot of random data, and observe a histogram of that, approximating the CDF). PAPR, not so much, because it is a peak statement, and with large OFDM-like symbols, the search space might simply get too large, and unless you do more math than you'd do to find the PAPR directly analytically, you can't make a statement about how close to the actual PAPR your monte-carlo based PAPR is.


How to Generate White Gaussian Noise with Known PSD in MATLAB

The specific PSD is basically white noise which was filtered with a filter which its Magnitude is the same as the PSD.

This is a result of the Wiener Khinchin Theorem.

So if you have the shape of the PSD do the following:

Create a filter which in the Frequency domain has the magnitude (More accurately the squared magnitude $ {\left| H \left( \omega \right) \right|}^{2} $) of the PSD you want (Namely, you don't care about the phase).
Create a White Noise.
Filter the white noise with the filter created in 1.

I would note that in the case of a White Noise the PSD is constant.
Hence in this case all you need is to generate white noise with normalized variance and then to multiply it by a constant (Changing its variance) to the value you're after.

Enjoy...


Negative instantaneous frequency with hilbert transform using scipy hilbert

It looks like artifacts due to the derivative. I used this code in Octave:
fs=300;
t=[0:1/fs:2];
c=chirp(t, 20, 2, 100);
s=c.*(1+0.5*cos(2*pi*10*t));
h=hilbert(s);
m=abs(h);
a=diff(unwrap(arg(h)))/2/pi*fs;
plot(m,"",a)

and this is what comes out:

I also tested the equivalent of this in LTspice and, without any form of unwrapping, this is the result:

where V(mag) is the result of hypot() (m), V(signal) is s, V(freq) is the frequency for the FM, and V(demod) is a. All are delayed (for matching I/O) by 23/2/200 which represents the order of the Hilbert transformer (Hamming window), divided by two (group delay) and by the sampling frequency (200).
My guess: don't take for granted the result of the differentiator and try to take a finer comb to it (is this a pun? not intended).


The signal, whose total energy you want to calculate, is periodic therefore it will have infinite energy...

To see that note, the following Fourier transform pair:
$$
x(t) = \cos(2\pi f_0 t) \longleftrightarrow X(f) = 0.5 \delta(f + f_0) + 0.5 \delta(f - f_0)
$$

And based on Parseval's relation,
$$
E_x = \int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df
$$

you can conlude that total energy is infinite


This is where things get ugly:

Can you define what "one-sideband noise power in 1 Hz bandwidth" is?

In my world, it'd just be the PSD of that signal multiplied by 1 Hz. But then, there'd be no need to call it "one-sideband noise power"; it'd just be the noise power density in 1 Hz.

So, "customary" might very much be a "relative" statement. I'd expect that term to only show up in literature that tries to explain signal properties based on a real-valued, analog technological background (maybe trying to bring the complex signal's view to an audience of people that learned the real-valued signal theory first).
 Frankly, I find it an elegant move that much of the modern literature starts right off with complex signals, and just shows that a simple multiplication with a complex sinusoid can move it to baseband, and then introduces the "normal" (as in: non-quadrature) mixer based on that understanding.

This is probably where that "upper and lower sidebands overlap" comes from:

If you take a noisy bandwidth $B$ of passband around a carrier at $f_0$ (you observe $]f_0-\frac B2; f_0 + \frac BS]$), and mix it down using a quadrature mixer (giving you I and Q signals, which we commonly interpret as real and imaginary part of a complex-valued signal), you get signal in $]-\frac B2; \frac B2]$.
The PSD of noise in that range then is the same as the PSD of noise in the passband. That happens because, mathematically, the quadrature mixer just multiplies with $e^{-j2\pi f_0}$, leading to a "clean" shift of the positive frequencies down to baseband (by $-f_0$).
If you, however, just multiply with a $\cos(2\pi f_0)$, you'll mathematically shift the positive to baseband, but also you shift the negative frequencies, which are a mirror of the positive frequencies, up (by $+f_0$).
Now you have the upper sideband being overlaid with the mirror of the lower sideband.
So, looking at this differently, this is an effect of using real-valued signals to represent something that would become complex baseband when mixed with a quadrature mixer.

Note that the statement is obsolete (or plain wrong), because frankly, unless you really want this for some reason, measurement devices either

mix onto an intermediate frequency (IF), so that none of the signal "crosses" 0 Hz when being downconverted, or
employ quadrature mixers, where that overlap doesn't happen.

Update I'll have to revise my above statement about the obsolence of the statement, because your book,

  Phase Noise and Frequency Stability in Oscillators by Enrico Rubiola

is a book on mixer theory, where you of course analyze the workings of the elementary mixers (which can be combined to form e.g. quadrature mixers), so that the author inherently works with the more "physical" real-valued signals, and then he's obviously right:

if you're characterizing the mixing element itself, you can't tell noise above and below the LO apart.


Problem understanding steps towards defining the convolution sum

The point of the whole exercise it to prove that you can calculate the output of an LTI system by convolving the input with the impulse response of the systems. This is done in the following steps.

Rewrite the input as a sum of delay delta impulses
Put that sum through the systems
Exploit LTI, swap the order of applying the system and the sum. Response of sum = sum of responses !

Let's go through the last step a little slower. We have
$$H\{x[n]\} = H\left \{ \sum_{k=-\infty}^{\infty} x[k]\cdot \delta [n-k]\right \}$$

First let's sawp summation and system, assuming that the response to the sum is the same as the sim of the responses (per LTI definition)

$$H\{x[n]\} =  \sum_{k=-\infty}^{\infty} H\left \{ x[k]\cdot\delta [n-k]\right \}$$

Now comes the tricky part: you have two time-like variables in the thing you need to transform: n and k. k is just the counting index for the different basis functions. n is actually the transform variable. You can think of this as having an infinite number of functions and you need to calculate the response for each one of them (and then sum the results). However, the response of each individual function is independent from the other ones. Hence, you can pull anything that's only dependent on k out from the transform. Different way of of looking at it: if H{} where an integral your integration variable would be $\partial n$ not $\partial k$. So we have

$$H\{x[n]\} =  \sum_{k=-\infty}^{\infty}  x[k]\cdot H\left \{\delta [n-k]\right \} = \sum_{k=-\infty}^{\infty}  x[k]\cdot h[n-k]
$$

where $h[n]$ is the impulse response of the system


I do not think there is a method to get a perfect reconstruction of the original signal by using a trimmed/cropped convolution. However, I am hoping there are approximations/assumptions you can make like the measured signal has finite energy.

If that is the case, what you can do instead is to find the inverse of the kernel, which you can then convolve with the measured signal.

You will, however, have some minor amount of error near the end, and this is unavoidable unless you have an infinitely long signal.

Following is what I did with your example:

inv_kernel = np.fft.ifft(1/np.fft.fft(kernel))

...

sig_deconv_full_2 = scipy.signal.convolve(sig_conv_real, inv_kernel, mode='full')
signal_estimate = sig_deconv_full_2[:NPOINT] #signal trimmed to the required length

This is the script run with your original signal. The squared error term is the squared error of signal_estimate and signal.

The impulse is expected as the measured signal suddenly ends. We can just trim the signal to get rid of the impulse.

I also tested it with a random signal I generated using NumPy's random library.


Correct way of zero padding in time domain

Zero padding changes the bin size of the DFT resulting a finer frequency resolution. It is done by appending zeros to the end of the signal because we are only artificially increasing the length to decrease the bin size, $\frac{2\pi k}{N}$. An $N$ point DFT is just evaluating the DTFT on certain frequencies. If you were to append some zeros to the front of the signal then you are doing a time shift, and a shift in time domain is modulation in the frequency domain: $x[n-k]=X(\omega)e^{-j\omega k}$.


Ripple in a SciPy Butterworth SOS Filter

I can find 2 errors at a quick glance :
1 - You should discard the transient when measuring the steady-state gain at for a given frequency $f$. I haven't checked how long the transient last, but you should factor it in your RMS measurement.
2 - The RMS calculation only works when the number of samples corresponds to a whole number of periods. Otherwise, the calculation is biased.
Edit :
Period in this context means the signal period. As you are in the discrete domain, this is the ratio between the sampling frequency and the signal frequency
$T = \frac{f_s}{f}$
For transient, see this site
https://www.dsprelated.com/freebooks/filters/Transient_Response_Steady_State.html
Edit 2 :
Assuming your signal frequency is 1 kHz
$T = \frac{44.1 kHz}{1 kHz} = 44.1$
Since you have 44.1 samples per period, the RMS calculation with only 44 samples will have some bias. However, if you use 441 samples, you will have exactly 10 signal periods and then the RMS calculation will have no bias.
As for the transient, you should plot the time response. You should be able to evaluate how long the transient lasts.


Approximate IIR filter with FIR filter with energy restrictions

This could be homework, so I'll only give you a few hints to help you solve the problem yourself.

Remember that a ZF equalizer just inverts the channel, so if $D(z)$ is the equalizer's transfer function, what must be the result of the product $C(z)D(z)=?$. From this equation you obtain $D(z)$, which is IIR.
Let the coefficients of $D(z)$ be $d[n]$, i.e., $$D(z)=\sum_{n=0}^{\infty}d[n]z^{-n}$$ The energy of the filter is $$E_d=\sum_{n=0}^{\infty}\big|d[n]\big|^2$$
The FIR filter you're looking for just uses the first $N$ coefficients $d[n]$. If $N$ is large enough, the approximation will be sufficiently good. According to the requirement, you need to choose a number $N$ such that $$\sum_{n=0}^{N-1}\big|d[n]\big|^2=0.8\sum_{n=0}^{\infty}\big|d[n]\big|^2$$


Improvement of SNR by sampling different number of cycles using FFT of a fixed length singal

Ok, so, what you need to understand is that FFT has particularities that if not taken into account produces erroneous conclusions out of its results. Firstly, both of your signals have the same SNR. Why? The way you described the problem led me understand that you tried to put two identical signals into different process paths to understand why they came out differently. Am I right? So, what are you really trying to ask is “Why both paths led to different results?”. If that’s the truth, allow me:
- FFT can generate lobes around the real frequency in many cases: 1. when you are trying to see a frequency that is not part of the indexes when you selected the FFT lenght; 2. When you don’t use a proper window function prior to FFT analysis; 3. When you have a small window relative to the frequency you’re trying to see and you produce a discontinuity. This is because the FFT algorithm considers that everything you have in your analysis window is a periodic pattern. Probably the 3rd reason is letting you believe that your second measurement has a higher SNR when the truth is your window has only 10 periods of your frequency, and it’s much likely that you have a discontinuity between the beginning of it and the end. If you fix this discontinuity you should see a clear FFT in your second case.


Output of a system given it's transfer function and input (beginner)

For an LTI system, output $y(t)$ is given by
$$y(t) = h(t)\otimes x(t)$$
Where $x(t)$ is input and $h(t)$ is impulse response of the system. The operator $\otimes$ represents convolution.

Convolution operation is mapped into multiplication in Laplace domain. ie,
$$Y(s) = H(s)\times X(s)$$
Where, $Y(s)$, $H(s)$ and $X(s)$ are the Laplace transform of $y(t)$, $h(t)$ and $x(t)$ respectively.

You can use any one of the above equation to solve your problem.

Two ways:

Find $h(t)$ from $H(s)$ and use first equation. OR
Find $X(s)$ from $x(t)$, use second equation to find $Y(s)$ and then find $y(t)$ from it.

I prefer second method.


Relation between SER and the channel capacity

Yes, there is a relationship between SER and channel capacity.  The channel capacity equation is-

$$
C = B\log_2(1 + \frac{S}{N})
$$
where $C$ is the channel capacity in bits/s, $B$ is the bandwidth in Hz, $S$ is the average received signal power, and $N$ is the average noise power across the bandwidth $B$.

Looking at the equation, it is clear that increasing the SNR will increase the channel capacity.  Increasing SNR does, of course, also decrease the symbol error rate (SER).  Thus, there is an inverse relationship between SER and channel capacity.

  What is {Modulation Order x (1-SER) / symbol duration [bps]} called?

As far as I am aware it does not have a name.  Having a bit stream where some of the bits are good and some of the bits are bad is generally not very useful.  Shannon's theorem assumes that error correction codes are used to reduce the number of errors to an arbitrarily low rate.  This reduces the number of generated bits/s, but increases the reliability of those bits.  The "channel capacity" of a channel is the number of bits it can produce with an arbitrarily low error rate.


Inverse Fourier transform of complex exponential with frequency dependent shift

In the general case: no.
For the (inverse) Fourier transform of a function to be composed of a sum of countable diracs (i.e. to be discrete), the function needs to be periodic.
Your $e^{-j\omega\tau(f)}$ is not periodic (it can be periodic, if $\omega\tau(f)$ happen to be periodic with period rationally related to $2\pi$, but that's a pretty special case). Therefore, its inverse Fourier transform is continuous, and can't be represented by a sum of diracs.


How to get started with fingerprint recognition?

This  may help you.
Also, everything related to fingerprinting in here.
I think is a good documentation to start with and view some source code. Nevertheless, you should learn some other algorithms before you start coding.


Two steps,

do a simple histogram equalization to make the brightness a little
bit more even.
Then use canny edge detector (as suggested by Marcus Muller).

matlab code:

I=imread('Your Image');
G1=rgb2gray(I);
G2=histeq(G1);
E2=edge(G2,'canny');
figure imshow(E2,[]);

Result (I hope it is what your are looking for):


Let us define the terms first:

Downsampling means reducing the sampling rate of a signal. If there is energy outside the new Nyquist frequency, there will be aliasing.
Filtering means eliminating all energy in some frequency bands.
Decimation means anti-alias filtering a signal and then downsampling it.

In a typical SDR application, your captured bandwidth is larger than the actual band of interest. There are two things you want to do:

Eliminate all signals outside the band of interest (by filtering).
Downsample to the minimum required rate, in order to reduce the computational complexity of the baseband processing algorithms.

One common way to accomplish this is the way you describe: downconvert the band of interest to DC, or to a low intermediate frequency if DC offset is a concern, and then decimate.

There could be a few reasons why you're seeing that filtering is slower than decimation. One is that maybe your filter order is larger than necessary. Another is that there are algorithms that perform low-pass filtering and downsampling simultaneously and with very low complexity. See for example half-band filters.


After advice about detecting focus quality of objects in a photo detected using YoloV3

Here is what I would try on the source image:

Split your image into 5x5 pixel blocks (maybe 3 maybe 7, who knows?)

Create output image one fifth (third, seventh, ??) size

For each block
  For each color channel

    Find best fit plane
    Measure RMS of (pixel value-plane value)

  Next

  Set output pixel to  RMS(R,G,B)

Next

In blurry/plain areas the RMS should be low.  With sharp lines it should be high.

Ced

So, I thought I'd give this a go with 5x5 block.

Here are some results.

Clearly (pun intended), not only is this algorithm an edge detector, but it also a focus metric at those locations.

Those are kind of dim, so I multiplied by 5.

If nothing else, they are kind of cool looking, but I think you can tell where the well focused areas are.

The 5x5 takes a while to process, so I decided to try a 3x3 fit on a simple 4 point gradient estimate (Down from a 16 point gradient estimate).

Still (pun intended again), the blurry lady pic proves this is a focus metric.

Here is the relevant code called for each color channel:

'==================================================================
Public Sub FindFocus(ArgV As Float[]) As Float[]

        Dim w, h, x, y As Integer

        w = ArgV.Bounds[0]
        h = ArgV.Bounds[1]

        Dim theFocus As New Float[w, h]

        For x = 1 To w - 2
          For y = 1 To h - 2
            GoSub CalculateFocusAtPoint
          Next
        Next

        Return theFocus

'-------------------------------------------------------------------
CalculateFocusAtPoint:

        Dim dx, dy As Integer
        Dim a, b, c As Float

'---- Plane Estimate:  z = ax + by + c

        a = (ArgV[x + 1, y] - ArgV[x - 1, y]) * 0.5
        b = (ArgV[x, y + 1] - ArgV[x, y - 1]) * 0.5
        c = ArgV[x, y]

'---- Calculate the RMS of the NonPlanar

        Dim v, e, s2 As Float

        s2 = 0

        For dx = -1 To 1
          For dy = -1 To 1
            v = a * dx + b * dy + c
            e = ArgV[x, y] - v
            s2 += e * e
          Next
        Next

'---- Set the Value

        theFocus[x, y] = Sqr(s2 / 9)

        Return
End
'==================================================================

Finally, just because I could, I turned your pic into a "drawing".


I'll add another answer, even though MBaz's answer is correct, because I think that it doesn't actually address the problem you have with arriving at the final solution. Summarizing, you have to solve

$$\frac{1}{2\pi}\int_{-W}^W|P(\omega)|^2d\omega=0.9\int_{-\infty}^{\infty}\text{rect}^2(t/T)dt=0.9\cdot T\tag{1}$$

for $W$ (which is the bandwidth in radians). If you like, rewrite (1) using $f$ instead of $\omega$. $P(\omega)$ is the Fourier transform of the rectangular pulse:

$$P(\omega)=2\frac{\sin(\omega T/2)}{\omega}\tag{2}$$

If I understood your question and your comments correctly, then you understand all of this, but your problem is solving (1) for $W$. Note that the integral has no closed-form solution, unless you consider an expression including the sine integral $\text{Si}(x)=\int_0^x\sin(t)/t\;dt$ closed form. With a little help from WolframAlpha, Equation (1) can be written as

$$\frac{2}{\pi W}\left[TW\cdot \text{Si}(TW)+\cos(TW)-1\right]=0.9\cdot T\tag{3}$$

Equation (3) has no closed-form solution for $W$, so it is not surprising that you got confused. You could of course try to solve (3) numerically but I do not think that this is the idea of the exercise. What you need to know is that obviously the largest portion of the energy is contained in the main lobe of the spectrum, i.e. between $-W_0$ and $W_0$, where $W_0$ is the first zero crossing of $P(\omega)$. From (2) it is clear that

$$W_0T/2=\pi\quad\Longrightarrow W_0=\frac{2\pi}{T}\tag{4}$$

So this is the essential bandwidth in radians.
Checking the exact percentage using the left-hand side of (3) we get

$$\frac{2}{\pi W_0}\left[TW_0\cdot \text{Si}(TW_0)+\cos(TW_0)-1\right]=
\frac{2\;\text{Si}(2\pi)}{\pi}\cdot T=0.9028\cdot T$$

which seems close enough.

In sum, the essential bandwidth of a rectangular pulse is given by the width of  the mainlobe of its spectrum, so you only need to be able to calculate the first zero of the spectrum and you're done.


Differential Equation of a Band Pass Filter

I'm afraid your results are not correct. Let's have a look at the low-pass filter. The current $i(t)$ through $R$ and $C$ must be

$$i(t)=\frac{v_i(t)-v_o(t)}{R}=C\frac{dv_o(t)}{dt}$$

So we get

$$v_i(t)=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{1}$$

For the high-pass filter we have

$$i(t)=\frac{v_o(t)}{R}=C\frac{d(v_i(t)-v_o(t))}{dt}$$
which gives

$$RC\frac{dv_i(t)}{dt}=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{2}$$

Let's call the time constants of the low-pass and high-pass filters $\tau_L=R_LC_L$ and $\tau_H=R_HC_H$, respectively. For the band-pass filter we need a relation between $v_i(t)$ of the low-pass filter and $v_o(t)$ of the high-pass filter. If we use $v_{oL}(t)$ to denote the output of the low-pass filter, which equals the input of the high-pass filter, we get from (1) by taking the derivative

$$\frac{dv_i(t)}{dt}=\frac{dv_{oL}(t)}{dt}+\tau_L\frac{d^2v_{oL}(t)}{dt^2}\tag{3}$$

From (2) we have

$$\tau_H\frac{dv_{oL}(t)}{dt}=v_o(t)+\tau_H\frac{dv_o(t)}{dt}\tag{4}$$

and (by taking the derivative)

$$\tau_H\frac{d^2v_{oL}(t)}{dt^2}=\frac{dv_o(t)}{dt}+\tau_H\frac{d^2v_o(t)}{dt^2}\tag{5}$$

Plugging (4) and (5) into (3) we finally get for the band-pass filter

$$\tau_H\frac{dv_i(t)}{dt}=v_o(t)+(\tau_L+\tau_H)\frac{dv_o(t)}{dt}+\tau_L\tau_H\frac{d^2v_o(t)}{dt^2}$$

By the way, if you were allowed to use the Laplace or Fourier transform, everything would be a lot easier.


Compute the time derivative of a noisy digital signal?

If computational efficiency is important to you, and your signal's bandwidth is not to large relative to the $f_s$ sample rate, three simple causal differentiators that tend to attenuate high-frequency noise are described at DSPRelated.

A plot of the results of that article is below.

Applying these three filters to your example signal (or an attempt to mimic it) yields the following plots.

R Code Below

# Q27420
h_ref <- c(-1/16, 0, 1, 0, -1, 0, 1/16)
d2 <- c(0.5,0,-0.5)
d3 <- c(-3/16, 31/32,0, -31/32, 3/16)

H_ref <- abs(fft(c(h_ref,rep(0,1000))))
f2 <- abs(fft(c(d2,rep(0,1000))))
f3 <- abs(fft(c(d3,rep(0,1000))))

w1 <- seq(0,length(H_ref)-1,1)/(length(H_ref)-1)
w2 <- seq(0,length(f2)-1,1)/(length(f2)-1)
w3 <- seq(0,length(f3)-1,1)/(length(f3)-1)

n1 <- length(w1)/2
n2 <- length(w2)/2
n3 <- length(w3)/2

plot(w3[1:n3],f3[1:n3]/1.2, col="blue", lwd=5,type="l", lty=2)
lines(w1[1:n1],H_ref[1:n1]/1.63,col="black", lwd=5)
lines(w2[1:n2],f2[1:n2], col="green", lwd=5)

legend(0.0, 1.5, c("Proposed","Reference", "Difference"),
       col=c("blue","black","green"), lty=c(2,1,1), lwd=c(5,5,5));

x <- c(seq(-1.6,-0.4,0.01), seq(-0.4,-0.2,0.001))
xn <- x + rnorm(length(x))/10

f <- seq(0,1,0.1)
m <- f
d4 <- signal::fir2(100,f,m)

y1 <- filter(xn,h_ref, circular=TRUE)
y2 <- filter(xn,d2, circular=TRUE)
y3 <- filter(xn,d3, circular=TRUE)
y4 <- filter(xn,d4, circular=TRUE)

par(mfrow=c(4,1))
plot(xn)
title("Original signal")
plot(y4,lwd=4,col="grey")
lines(y1,col="blue")
title("Deriviative calcualted using the reference technique")
plot(y4,lwd=4,col="grey")
lines(y2,col="black", lty=2, lwd=3)
title("Deriviative calcualted using the proposed technique")
plot(y4,lwd=4,col="grey")
lines(y3,col="green", lwd=5)
title("Deriviative calcualted using simple difference")


I must admit I haven't worked on anything like this before, but looking at your data it seems to me that computing the median of e.g. the lower left corner point of all rectangles (i.e. taking x and y dimensions into account) should get rid of outliers. After having removed outliers from your data, you could simply try to average the remaining rectangles. Obviously, such an approach requires some heuristic tweaking, but it may be worth a try.

I just tried Matt's approach (just taking the 1D medians of the top/left/bottom/right X and Y values) and the results seem pretty good. See plot below of 20 realizations of various random rectangles.

scilab CODE ONLY BELOW

function r = getRect(t,l,b,r)
    large =  [1 1 0 0]*rand(1,1,'normal')*abs(t-b) + [0 0 1 1]*rand(1,1,'normal')*abs(l-r);
    small =  [1 1 0 0]*rand(1,1,'normal')*5 + [0 0 1 1]*rand(1,1,'normal')*5;

    noise = small;
    if (rand(1,1,'uniform') < 0.2) then
        noise = large;
    end
    r = [t l b r] + noise;
endfunction

for p = 1:20;

N = 20;

x_points = [2 2 4 4 2];
y_points = [1 3 3 1 1];

subplot(5,4,p)
RR= [];
for n=1:N
    R = getRect(10,10,100,50);
    RR = [RR; R];
    plot(R(x_points),R(y_points),'r');
end

RRm = median(RR,'r');

xset("thickness",4)
plot(RRm(x_points),RRm(y_points),'g');
xset("thickness",1)
end;


Steve Eddins has a nice page that explains it.

It comes down to being able to do a separate 1-D filters in each direction (requiring $MN(P + Q)$ multiplies and adds), rather than a single 2-D convolution ($MNPQ$ multiplies and adds).  Here the image is assumed to be $M \times N$ and the convolution kernel $P \times Q$.

Of course, this assumes you're doing everything in the time domain, and not using FFTs.


Is it possible to use MRC with DS-CDMA-MIMO

OK, let me explain that for you. First, if you want to understand anything in engineering, it's recommended to write it in mathematical form, then try to solve it theoretically and build its code accordingly. Second, regarding your question, is it possible to use MRC for MIMO - CDMA,   Yes, that's possible, why not?

Now, The code you provided is not right. Why ? all simply, first   check when using simple MIMO without spreading spectrum. You will see that first element in x, which is 0.7+.7i, is multiplied by the first column of channel H and so on. In case if you are using spreading spectrum, it means that instead of having one sample 0.7 + 0.7i, you will have N samples, where N is the code length. it means you will have N samples in the same duration of symbol before spreading.  So, all samples should be multiplied by the first column of H.

I think it's clear now, right?  so you should just now check how to write the above answer in matlab. . .  for example check this command :

r = [];
for kk = 1 : size(h,1)
                      r_1 = kron(x_s(:,kk),h(:,kk));
                      r = [rec r_1];
                  end

r = sum(r,2);

Good luck.


Direct form 1 filter equation for a biquad filter

This may be a misunderstanding based on a fixed point implementation. Fixed point implementations can benefit from "fractional multiplies" which requires a number format where all absolute values are  1 or less. In Q notation that would be Q15.0 for a signed 16-bit integer. Many processors have a dedicated hardware instruction for this operation, which makes it fast and efficient.
Assuming that all poles and zeros are inside the unit circle, we have $|b_0|, |b_2|, |a_2| < 1$ but $|b_1|, |a_1| < 2$. In order to avoid shifting associated with mixed Q notation, the $b_1$ and $a_1$ coefficients are split in have, so numbers can be kept in the same Q. Hence it would be implemented as
$$y[n] = b_0 \cdot x[n] + b_1/2 \cdot x[n-1] + b_1/2 \cdot x[n-1] + b_2 \cdot x[n-1] ...$$
or
$$y[n] = b_0 \cdot x[n] + (b_1/2 \cdot x[n-1]) \ll 1 + b_2 \cdot x[n-1] ...$$
where $\ll$ is the left-shift operator.
This may looks like the coefficient is multiplied by two but in reality it's first split in half to make it's magnitude smaller than $1$, so the same type of fixed point multiplication can be used for all operations.


Tricky problem. Any speaker to mic path in a room is actually a linear time invariant system so they can be fully described by impulse response or transfer function. However the transfer functions are very complicated (1000s of degrees of freedom) and vary tremendously with position in the room. The room is also a very effective de-correlator. Signals from the same speaker recorded at two different location will seem mostly uncorrelated.

You could try to use some statistical room acoustical properties: diffuse field power, reverb time, early decay time, reflection density, mean free path lengths etc. However, these are difficult to determine and not that much different between most rooms in the same house.

So unless you can deploy some significant constraints (control the excitation, control location of source and microphone), this will be very hard to do.


Carrier frequency offset estimation for QPSK and asymmetrical spectrum

A likely source of the amplitude imbalance given it is on the actual signal is filtering distortion along the signal chain such that the gain is not flat across the bandwidth of the spectrum. This would not be uncommon with analog filtering and can be compensated for with equalization (if needed).

For an alternate approach to carrier tracking applicable to QPSK, consider estimating the carrier frequency offset using the baseband samples after timing recovery (timing recovery with a Gardner TED converges even under relatively large carrier offset frequencies). With that, the correct sample locations for each symbol (just prior to decision) are used with a cross product frequency discriminator to estimate carrier offset as follows:

$$Err = |I_2| |Q_1|-|I_1| |Q_2|$$

Where the absolute values convert all samples to the first quadrant (so this will work for offset frequencies such that the phase doesn't rotate more than $pi/4$ radians between samples, or up to 1/8th the symbol rate.

$Err$ is proportional to frequency offset as the result above is the imaginary portion of the complex conjugate product of the two samples, which is proportional to the phase between the samples. A change in phase versus a change in time (the symbol duration) is proportional to frequency, since $f=d\phi/dt$.

Thus the cross product frequency discriminator is the cross product phase detector for the phase between two samples at two different times.


Autocorrelation for the sum of a deterministic and stochastic signal

The random process $x(t)$ doesn't have a power spectral density (PSD) in the conventional sense because it is non-stationary, due to the presence of the deterministic time-dependent function $d(t)$.

You can compute the auto-correlation function of $x(t)$ as follows:

$$\begin{align}R_x(t_1,t_2)&=E\{x(t_1)x^*(t_2)\}\\&=E\{[d(t_1)+s(t_1)][d^*(t_2)+s^*(t_2)]\}\\
&=d(t_1)d^*(t_2)+d(t_1)E\{s^*(t_2)\}+d^*(t_2)E\{s(t_1)\}+E\{s(t_1)s^*(t_2)\}\\&=d(t_1)d^*(t_2)+d(t_1)\mu_s^*(t_2)+d^*(t_2)\mu_s(t_1)+R_s(t_1,t_2)
\end{align}\tag{1}$$

where $\mu_s(t)$ is the (generally time-dependent) mean of the random process $s(t)$, and $R_s(t_1,t_2)$ is its auto-correlation function.

If you assume that $s(t)$ is wide-sense stationary (WSS), and if you further assume (for ease of notation) that $d(t)$ and $s(t)$ are real-valued, then Eq. $(1)$ simplifies to

$$R_x(t,\tau)=d(t)d(t+\tau)+[d(t)+d(t+\tau)]\mu_s+R_s(\tau)\tag{2}$$

where I've used the usual substitution $t=t_1$ and $\tau=t_2-t_1$.
Note that $(2)$ generally still depends on the absolute time $t$, and, consequently, $x(t)$ is non-stationary and has no PSD.

As already pointed out in a comment, only if $d(t)=d$ is constant (and $s(t)$ is WSS) does the random process $x(t)$ become WSS with auto-correlation function

$$R_x(\tau)=d^2+2d\mu_s+R_s(\tau)\tag{3}$$

and its PSD becomes (by taking the Fourier transform of $(3)$)

$$P_x(\omega)=2\pi d(d+2\mu_s)\delta(\omega)+P_s(\omega)\tag{4}$$

where $\delta(\omega)$ is the Dirac delta impulse.

EDIT: As mentioned in a comment by Dilip Sarwate, $P_s(\omega)$ also generally includes a Dirac impulse if the mean $\mu_s$ is non-zero. If you define a zero-mean random process $\tilde{s}(t)$ by

$$s(t)=\tilde{s}(t)+\mu_s\tag{5}$$

then $P_s(\omega)$ can be written as

$$P_s(\omega)=P_{\tilde{s}}(\omega)+2\pi\mu_s^2\tag{6}$$

Using $(6)$, Eq. $(4)$ becomes

$$P_x(\omega)=2\pi (d+\mu_s)^2\delta(\omega)+P_{\tilde{s}}(\omega)\tag{7}$$

which explicitly shows all Dirac impulses in the PSD of $x(t)$.


Single pole IIR filter, fixed point design

One thing to consider when implementing an IIR filter, whatever the order, is quantization and limit cycles.

Let me show you with a quick example with your original filter

$y[n] = a*x[n]+(1-a)*y[n-1]$

Let a = 0.005 and say that we use 16-bit signed coefficients.

$a_{fixedpoint} = a * 32768 = 164$

Let's assume that the input and output are 16-bit integers. They will go from -32768 to 32767. The previous output $y[n-1]$ is also a 16-bit integer in this scenario. The implementation would look like this, with temp being a 32-bit integer. All other variables are 16-bit integers.

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1]$

$y[n] = round(temp/32768)$

Let's try this implementation with a step response

It looks good, right? Perhaps we should zoom in.

We don't reach 32767 even though the gain of the original filter is 1. This is called a limit cycle. For order-1 IIR filters, limit cycles take the form of a DC offset that you can't get rid of. For IIR filters with complex poles, limit cycles take the form of a small oscillation that you can't get rid of. The problem is caused by the fact that $y[n-1]$ is stored on 16 bits. If we increase the resolution of $y[n-1]$ to 32 bits, it would solve the problem. Another solution would be to use fraction-saving.

If the input would go from 32767 to 0, we would have the same problem, the output would get close to 0 but would not actually reach 0.

Increased resolution

If instead of storing the output on 16 bits, we store it on 24 bits and name this variable acc (increased resolution of y[n]).

$ temp = a_{fixedpoint}*x[n]*256 + (32768 - a_{fixedpoint})*acc$

$ acc = temp/32768$

$y[n] = round(acc/(256))$

As we can see, with the increased resolution we can reach 32767 even though the coefficients are still on 16 bits. However the price to pay is having intermediate variables with greater resolution, the calculations might be slower.

fraction-saving

Finally there's a clever approach to minimize the added resolution needed for an IIR, especially when the poles are close to the unit circle. It is called fraction-saving. The idea is to memorize the round-off error of the previous sample and apply it to the next calculation to reduce the quantization effects. An additional variable is needed called error, but you don't need to increase the resolution of the previous output. So if x is represented on 16 bits, then y[n] can be stored on 16 bits too. The error variable would be represented with 16 bits too in this example. As Robert Bristow-Johnson pointed out, fraction-saving is a form of noise-shaping

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1] - error[n-1]$

$y[n] = round(temp/32768)$

$ error[n] = y[n]*32768 - temp$


Can Aliasing Happen in Laplacian Pyramids?

The reconstruction will remove aliasing, no matter which filter kernel is used.

The lower layer of the 2-layer Laplacian pyramid on page 8 of the Toronto lecture notes, adapted here:

Figure 1. A 2-layer Laplacian pyramid.

can be redrawn as:

Figure 2. An alternative presentation of the lower layer of the Laplacian pyramid of fig. 1.

This layer is an identity operation: what is subtracted is also added. This allows to draw the next (and in this case the only) layer above it the same way and it also can be identified to be an identity operation. And so on until the whole process becomes an identity operation.


How to apply "Substitution Method" in Difference Equations?

Substitution is not a method to obtain the particular solution $y_p(n)$, but the author simply claims that $y_p(n)=e^{-\beta n}u(n)$ is a particular solution of the recursion (1.70), and encourages you to check his claim by substituting $y_p(n)$ for $y(n)$ in Eq. (1.70) and see if it's correct or not:

$$e^{-\beta n}u(n)=e^{-\beta}e^{-\beta (n-1)}u(n-1)+\delta(n)=e^{-\beta n}u(n-1)+\delta(n)\tag{1}$$

Now check if (1) is satisfied. Distinguish 3 cases: $n<0$, $n=0$, and $n>0$. As a result you should be able to see that (1) (and, therefore, Eq. (1.70)) is satisfied for $y(n)=y_p(n)$, showing that the suggested sequence $y_p(n)$ indeed is a particular solution of (1.70).


How to warp a pixel patch to another camera

$x  K_{inv}$, with $x(3) = 1$, produces the pixels as 3D points "stuck on the screen" of the first camera. After transforming you get the same 3D points, but in the coordinate frame of the second camera, and then you project in that camera. So, overall, what you are doing is looking at the first camera from the point of view of the second one.

I suspect that is not what you want....


information filter instead of kalman filter approach

Information filters are an application of a learning algorithm. One example of an information filter would be a recommendation system. On the other hand, a kalman filter is restricted to linear tracking. Its not learning per se.

If you are looking for machine learning/AI based approaches to SLAM, then you should probably look at Particle Filters.


In the case $a$ is a random variable (r.v.) we distinguish between two types of unbiased: (hereinafter I'm switching to $\theta$ and $\hat{\theta}$ instead of $a$ and $\hat{a}$ respectively):

Unbiased in the strong sense: given by $E\left \{ \left (  \hat{\theta} - \theta \right )| \theta\right \} = 0, \quad \forall \theta$. Most often tha n not it is hard to satisfy this definition.
Unbiased  in the weak sense: given by $E\left \{  \hat{\theta} - \theta \right \} = 0, $ where the expectation is executed on all random terms inside. The meaning of this criterion is that you want the estimator to be unbiased "in the mean" and not for all values of $\theta$. This criterion can be calculated using the law of total expectation (smoothing theorem).


What is the value of last sample of a discrete sinewave

No.  Assuming the period of the signal is N samples since your indexing is zero based the $Nth$ sample will be indexed as $N-1$.

Your signal would be modelled by:

$$ x[n] = A \sin \left( \frac{2\pi}{N}n \right) $$

Thus $x[N]=0$ and $x[N-1] = x[-1] = -x[1]$.


Determining optimal binary decision rule threshold from observations with unknown priors?

My intuition is that it would be difficult to derive the right decision threshold you expect to find:

$$\tau = \frac{1}{2}\left(\mu_0 + \mu_1\right) - \frac{\sigma^2}{\lVert\mu_0 - \mu_1\rVert^2} \log \frac{\pi}{1 - \pi}\left(\mu_0 - \mu_1\right)$$

From the global statistics you are considering (sample mean: $\pi \mu_0 + (1 - \pi) \mu_1$ ; standard deviation: more complex expression but I doubt it would involve a log).

I would approach the problem this way:

If the assumption that $\sigma$ is small can be made

I'm mentioning that, because keep in mind that the decision threshold is affected by $\pi$ only if $\sigma$ is sufficiently high to allow both classes to overlap. If the $\mu$s are distant by more than a few $\sigma$, class prior probabilities have nothing to say in the decision process!

Run k-means on your observations ($\sigma$ is small and is shared by both classes, so k-means is in this case EM for the mixture model). If you just want to binarize these observations and no other data, you can stop here.
If you have new observations to binarize, and you know they are generated by the same process, you can use the class centroids found by k-means on your training data as estimates of $\mu$, and use the middle as a decision threshold.

If no assumption about $\sigma$ can be made

Run the EM algorithm (with pooled, diagonal covariance) on your training data. Use the inferred "soft class membership" variables to binarize your observations.
Compute the decision threshold $\tau$ from the parameters given by EM to binarize new data generated by the same process.


Why does a higher sampling frequency mess up my bandpass filter?

With an IIR filter, the higher the ratio is of sampling frequency to bandpass filter frequency, the closer the poles are to the complex unit circle in the Z-transform.  Poles close to the unit circle can produce numerical instabilities (ill conditioned computations) when computing the filters response.  This is because any numerical quantization or "rounding errors" can get amplified in the arithmetic recurrences if there is only a tiny difference in location between the pole and the unit circle, even if the pole is "supposed" to be inside.

A better solution is to first downsample the signal before filtering to use an IIR filter that is possibly of a lower order and more stable.  You can upsample back up, if needed, after the filter has been applied.


Difference between Uniform Linear Array (ULA) 3 dB beamwidth and bearing resolution

The answer will be NO, but to be precise especially in case of conventional beam forming.

Theoretical Derivation

The antenna electric field pattern of array antenna can be given by

where N is the number of antenna elements

d is the spacing between antenna elements

In order to find the beamwidth (3 dB), the above equation should be equated to
  and solve for

The solution will  come to be as

where D is the total aperture distance and can be approximated as

For spacing of   the equation simplifies and can be approximated as

Thus beam width of planar antenna can be represented as

Angular Resolution Using FFT Over Antenna Elements

In order the prove the lemma that the angular resolution obtained by performing the FFT across antenna dimension equals to the beamwidth of the antenna array, we have to obtain the angular resolution using FFT.

The below diagram shows the frequency obtained due to path difference between the antenna elements which occurred due to the angle of arrival other than the broadside angle

The frequency resolution using FFT can be represented as  where N is the number of samples but in our case it is equal to number of antenna elements.

The angular resolution   can be found by equating the difference in frequency of different angles
 and

The frequency resolution becomes:

Solving we get

and if we put  and
we get the same resolution as the beamwidth i.e

Summary

Thus for summarizing the above discussion using FFT we can only achieve the max angular resolution equals to the beamwidth of the antenna array which is equal to angular resolution = antenna array beamwidth = 2/N

Kindly note that this equality is true for the case of conventional beam forming only, there exist also some advance beam former where your angular resolution is better than the 3db beam width

This is the same answer which I provided on
ULA beamwidth and angular resolution, might this will be helpful.


If group delay is a constant for all $\omega$, Does this system have a linear phase response?

This is a homework type question, so I will only give hints and no solutions. You should know that the group delay is the negative derivative of the phase (with respect to frequency), so if the phase is a linear function of frequency, the group delay must be a constant. This should answer your first question.

Concerning your second question, why don't you just compute the phase of the given frequency response? It may be helpful to rewrite $H(e^{j\omega})$ as

$$H(e^{j\omega})=2e^{-2j\omega}\left(1-\frac{e^{j\omega}+e^{-j\omega}}{2}\right)\tag{1}$$

and to realize that the term in parentheses is real-valued.

EDIT:

The frequency response $H(e^{j\omega})$ can be written as

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

where $M(\omega)$ is the magnitude, and $\phi(\omega)$ is the phase. Now compare $(1)$ and $(2)$ and figure out what $\phi(\omega)$ must be.


When a system is linear, effects are proportional to causes. You can multiply the input by any constant, you'll get the same output multiplied by the same constant.

If you look at how the rise time is defined, you'll see that it is independent of any amplification factor.

You can also add two inputs, and the output will be the sum of the separate outputs corresponding to the separate inputs.

So it you apply a second, delayed step, there can be two situations

steady state (constant output) has been reached after the first step: then you can consider two distinct rises with equal times;
steady state has not been reached yet: there is a single, combined rise, with a different rise time.


Assume $Y = g(X)$ be the function of RV $X$, then by using the following

$$E\{ g(X) \} = \int g(x) f_X(x) dx $$

variance of $Y$ can be computed without the computation of pdf $f_Y(y)$ as:

$$
\begin{align}
\text{Var(Y)} &= E\{ (Y-\mu_Y)^2 \} = E\{ Y^2 \} - (\mu_Y)^2 \\
& =  E\{ g^2(X) \} - E\{ g(X) \}^2 \\
& = \int g^2(x) f_x(x) dx  - \left(\int g(x) f_X(x) dx \right)^2 \\
\end{align}
$$


Matlab code for Simulating of Phase Noise for OFDM systems

Below is a short code to do this for 1024 samples. "phi" array below is the required phase for each samples

N=1024; %number of samples

phi=zeros(1,N);
phi(1,1) = pi/4;. %initial phase, could be set to zero if the phase is set to zero by the front end processing

Phase_snr=10; %snr in phase estimation,equivalent to 1db

for n=1:N-1
phi(1,n+1)=phi(1,n) + (1/sqrt(phase_snr))*randn(1);
end


I tried http://www.codeproject.com/Articles/196168/Contour-Analysis-for-Image-Recognition-in-C once. It did not work well for me, but that guy seem to have better results.


Personally, I would just use a lowpass filter  try and keep it as small as possible to reduce the effect of asymmetry.

The fact that there is significant noise already implies some uncertainty in the peak position, I expect this will be of similar or greater magnitude to any bias the filter introduces.

If you are really concerned about the bias you could look at using morphological opening and closing filters. This is basically running a ball over/under the surface so either one will over/underestimate the peak height. You can then get 2 lines to give an envelope of where the peak may lie.  The average of the two curves is probably a decent estimate of the "true" curve prodived the noise frequency is suitably high.


You would have to connect both cameras to your computer, and either create two webcam objects or two videoinput objects to read images from them.


what is the type of filter that has values ones?

This, on its own, would perform a summation operation of L elements of your input. If you modified it to p = ones(1,L)./L, then it would become a moving average (or in other words a low-pass) filter.

For more information please see this and this link.

Hope this helps


Gesture recognition with Hidden Markov Models

I have used HMM for gesture recognition (not pose recognition). What I did was: tracking the hand and recognize the gesture the hand was drawing in the air, you can image it as a trail.

You can use HMM as sequence recognizer, so first of all you need to transform your image into a discrete number sequence.

For each gesture you want to recognize, you need to train an HMM for that.

So you have a dictionary with some known word. Each one is a trained HMM. If you have a new word (unknown obdervation) you can compute the probability for each word of the dictionary to be likely the unknown one.

Some pseudo-code:

##################### training phase

examples = [112233, 11233, 123, 1122223333]
word1 = train(example)

examples = [222333, 22222223333, 2222333333]
word2 = train(example)

examples = [124555, 1122445, 1111111222224444555]
word3 = train(example)

dictionary = [word1, word2, word3]

##################### recognition phase

#let's say you have a new unkown word: 12245, you want to know what word of the dictionary it is more likly

unkown = 12245
probabilities = []
for w in word:
    probabilities.append( calculate_prob(unkown) )

Now, see what the is the maximum value in probabilities and you get the most likely word of the dictionary!

See here:

http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf

http://www.codeproject.com/Articles/69647/Hidden-Markov-Models-in-C

http://www.creativedistraction.com/demos/gesture-recognition-kinect-with-hidden-markov-models-hmms/


From your description, here is what is happening in the video:

The short-term Fourier transform (aka spectrogram) of the signals is computed. The output of this operation is a matrix of complex values, which cannot be represented as images. Thus, the magnitude or the square of the magnitude is extracted to yield a single positive real value converted into a pixel intensity.
The magnitudes (pixel intensities) are multiplied.
The inverse short-term Fourier transform of the product is synthesized, presumably with made-up phases (or with the original phase of the carrier).

There are so many reasons why this has nothing to do with a convolution. I'm sure other members will point a few more:

What is multiplied are magnitudes or their squares; but not the actual complex values produced by the Fourier analysis.
Even if the quantities that were multiplied were the actual complex values, keep in mind that when dealing with discrete signals of finite length, the inverse discrete Fourier transform of the product of the discrete Fourier transforms is the Circular Convolution, not the convolution.
And still... If you split two signals $x$ and $y$ into blocks of length N, and compute the pair-wise circular-convolution of length N of these blocks, you'll get something very different from $x * y$.

Convolving two audio-signals is a rather meaningless operations. One usually convolve an audio signal with the impulse response of a system; and when convolution is performed in the frequency-domain, the details are trickier than just pair-wise multiplications of FFTs (there's some blocking issues; + necessary overlap add or save).

What you saw would be more accurately described as spectral cross-synthesis. And there is one major reason why this is a very different beast from a vocoder. The goal of a vocoder is to apply the spectral envelope of the modulator onto the carrier. I emphasize on spectral envelope, because when you use a vocoder to make a stack of sawtooth waves say "hello", you apply the formants and overall loudness envelope of the original speech signal to the saw waves, but the last thing you want is the pitch information of the speech signal to get involved. The analysis filter bank of the vocoder should be designed to make abstraction of the individual spectral peaks in the modular signal - what matters is the rough spectral envelope - the location of the bumps (formants). This is why a vocoder doesn't need more than 12-30 channels - too few channels and it doesn't capture the formants, too many channels and it starts capturing pitch-dependent fine spectral peaks. Just like the features used in speech recognition...

Let me give you an example. Let us say you have a modulator speech signal with a 200 Hz f0 ; and formants at 1kHz and 1.5kHz - its spectrum is a sequence of narrow peaks at 200 Hz, 400 Hz, 600Hz, 800 Hz, 1000 Hz, 1200 Hz, 1400 Hz, 1600 Hz, 1800 Hz ; with the peaks at 1000 Hz and 1400 / 1600 Hz emphasized (formants). Let's say your carrier signal is a 140 Hz sawtooth - the spectrum is made of narrow peaks at 140 Hz, 280 Hz, 420 Hz... with a decreasing $1/n$ amplitude. If you go with what you described in your question (product of STFTs), you wouldn't have much left because the two spectra have very little overlap - the first common frequency in the sequence would be at 1400 Hz! What you want to do is somehow capture that the modulator spectrum has a bump at 1kHz and 1.5kHz, and use this information to boost the carrier spectrum in this frequency area. That's how vocoders work.

To do so, and if you really want to go the STFT route, an option would be to smooth the modulator spectrum with a rather wide kernel before doing the multiplication - so that what you are really doing is applying the spectral envelope of one signal onto the other one. This would be akin to applying a motion-blur on the Y axis in Photoshop using this silly image manipulation presentation...


Spectrogram Time Resolution with Overlap

A sliding FFT (maxed-out overlap of N-1) used for a spectrogram will turn a single sample event into a fat wide hump (shaped like the window applied) of width 2N-1.  You might call the fat wife blurry peak of a single hump “resolution”, but multiple overlapped events of different heights plus a bit of noise might make it impossible to “resolve” (as in separate) adjacent event peaks in time.

To make the blurry humps for short events appear narrower you have to use shorter FFT windows, no matter what the overlap.

A 50% overlap does provide more time resolution, as the typical quantized Von Hann window eats edge events.  The extra windows in a 50% overlap brings that info back up into the picture, as the sum of 50% overlapped Hamming or Von Hann windows has a constant gain.


The algorithm sets the high frequency DFT coefficients to zero hence effectively providing a lowpass filtering on the time domain signal. Your problem was on the line of creating Y1, which I replaced with the corrected form. I've changed the code a little bit also.

% part 1:
[y, f] = audioread('C:\path to your wav file');
y = y(1:1024*64);                % select a portion of it.

Ya = fftshift(fft(y));           % Compute DFT Ya[k] and SHIFT it
N = length(Ya);                  % Length of DFT Ya[k]

figure,plot(-f/2+f/N:f/N:f/2 , abs(Ya))  % plot the DFT magnitude.

% part 2:
Yr = zeros(1,N);                % Yr[k] of length N!
Yr((N/4)+1:(3*N/4)) = Ya((N/4)+1:(3*N/4));  % Assign nonzero coefficients

% part 3:
yr = real(ifft(fftshift(Yr)));  % RECONSTRUCT the filtered signal.
                                % Take REAL part for convenience.
figure,plot(-f/2+f/N:f/N:f/2, abs((Yr)));

% part 4:
sound(yr,f,16);                 % Listen to the reconstrcusted audio


Polyphase Barker sequences have been studied (Google throws up multiple hits), but from a practical viewpoint, it is usually simpler to use binary Barker sequences on the I and Q channels.


Help generating coefficients for basic FIR LPF

These even aren't great coefficients for any software implementation that supports floating point – you could easily omit the last and the first coefficient without changing the filter significantly; the rest of the coefficients seem relatively well-conditioned.

So, basically, what you do is that you pick a fixed point bit width $B$ that suits the signal you're dealing with and the precision you need, and multiply each coefficient with $k=\frac{2^{N}}{\max\limits_{x\in b}|x|}$ (or $2^{N-1}$, as we're probably talking about signed numbers); you then usually just round, and try your filter with pseudorandom white noise, as you'll be inferring some degree of quantization error. Alternatively, you analytically calculate the frequency response in $\mathcal Z$-domain (freqz does that for you) of the quantized filter. Of course, you'll see an additional gain of $k$, but DSP engineers typically just "carry through" such factors (as most operations are linear), and interpret the results accordingly at the very end.

Now, as said, your coefficients (aside from the edge ones) look relatively well-conditioned, and I daresay that if you for example decided that 24 bit are enough precision for your coefficients, then that quantization error might be negligible.

If, however, the resulting frequency response is suboptimal, you can increase bit depth of the coefficients in bit steps¹, or by selectively increasing or decreasing symmetrical pairs of taps. For example, if you find that your frequency response is distorted to have too much gain at low frequencies, it might be a good idea to alternatingly in- and decrease adjacent coefficients by one (be slow with that, and always make sure the errors you're fixing even matter at all – who cares if your passband ripple is 0.2 dB instead of the floating point 0.18 dB, if all you needed was 0.4 dB?).

Without loss of generality, yes, there's also methods to generate integer taps directly (they shouldn't be much better for the window method you're using than rounding afterwards, but might be better for iterative approaches), but I'm not versed in these, or would even know of any toolbox that does them.

¹: as FPGAs really do not care whether you're using 23, 24, or 25 bits – unlike CPUs, they are not inherently limited to specific bitwidths; now, this is a bit of a lie, as for multiplication you might be using "DSP slices", which do have finite bit width and are fixed, not implemented on demand. However, it's thoroughly possible that your synthesizer might decide that for "fixed factor" multiplications, a straightforward combinatorial multiplier is appropriate. Maybe don't worry about these details right now.


Relation between a Fourier series harmonic component and its power

See first of all, you have confused magnitude spectrum to power spectrum. The figure you have shown in green is not power but amplitude spectrum hence y axis is represented as "Magnitude"

what is happening in the figure you have shown is the periodic signal is expressed as a trigonometric series where

$$x(t)= a_0 +\sum_{k=1}^{\infty} [a_k\cos(k\omega_0t) +b_k\sin(k\omega_0t)]$$

where $a_k$ and $b_k$ are trigonometric Fourier series coefficients and $\omega_0$ is fundamental period of signal x(t).

Now to express the amplitude and phase of the periodic signal in terms of Fourier series coefficients, trigonometric series representation is inconvenient as determining phase and magnitude of spectrum cannot be done directly from the above mentioned equation.

So we generally express periodic signal in terms of polar form of Fourier series where the polar form is given by
$$x(t) = M_0 + \sum_{k=1}^{\infty} M_k \cos(k\omega_0t - \phi_k)$$

where from the above expression you can see that $M_k$ denotes amplitude at kth harmonic and $\phi_k$ represents the phase of the kth harmonic.

where
$$M_0 =a_0$$
To find the remaining terms $a_k$ and $b_k$ in terms of $M_k$ and $\phi_k$, let's expand the summation of both the expressions of $x(t)$ and equate them.

i.e.  $$a_k\cos(k\omega_0t) +b_k\sin(k\omega_0t) = M_k\cos(k\omega_0t-\phi_k)$$

now expand the right hand side equation which is of form $\cos(A-B)$

So $$a_k \cos(k\omega_0t)+b_k \sin(k\omega_0t) = M_k\cos(k\omega_0t)\cos(\phi_k) +M_k \sin(k\omega_0t)\sin(\phi_k)$$

so equating coefficients on both sides we get

$$ a_k = M_k \cos(\phi_k)$$

and $$b_k = M_k \sin(\phi_k)$$

now to find magnitude $M_k$, square and add both the above terms, then we get

$$M_k =\sqrt{a_k^2 +b_k^2}$$

and phase is obtained by dividing the expressions for $a_k$ and $b_k$

$$\phi_k = \arctan(b_k/a_k)$$

Now the graph you have shown in green is the plotted $M_k$ values.

Note that you can also find magnitude and phase using exponential Fourier series coefficients, but then you will have both positive and negative harmonics in exponential Fourier series unlike trigonometric Fourier series which has only positive harmonics. So polar form can be seen as a direct equivalence with the trigonometric Fourier series coefficients.


Can I sample at Nyquist rate if I know that my samples are taken only at the signal's maxima or minima?

when you have a quantized signal (like in your case) you also have a noise floor in the frequency spectrum. This noise floor is proportional also to the sampling rate. If you sample a signal at twice its maximum frequency you can distinguish all the frequency components above the noise floor. In your case you are sampling at $F_p = 2$ because there are 2 samples every second therefore you are assuming your signal has a maximum period of 1.
There is one case when you are interested on sampling at the same frequency of the signal. For example if you know that your signal has an unwanted ripple at the frequency $f_{ripple}$ you can sample at this frequency $F_p = f_{ripple}$ preserving the lower frequencies while the ripple is transformed into a constant error which basically depends on the ripple waveform.


Confused about anti-aliasing in analog filters

Aliasing does not happen in the analog world. It only happens when you sample a signal.

If you have an ADC with sampling frequency $f_s$, then you'll have aliasing if the input signal has frequency components larger than $f_N=f_s/2$. To avoid aliasing, it's common to low-pass filter the signal before it is input into the ADC. That way, you make sure that the digital signal out of the ADC has no aliasing.

Such an analog filter is called an "anti-alias filter", not because it will experience anything related to aliasing itself, but because its purpose is to avoid aliasing in the ADC.


Indeed you can do that.
You may look on my answer to How to Prove a 2D Filter Is Separable?

By the SVD for any filter $ A $:

$$ A = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Since we're talking about separable filter then:

$$ A = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

So the columns filter is $ \sqrt{\sigma}_{1} {u}_{1} $ and the rows filter is $ \sqrt{\sigma}_{1} {v}_{1}^{T} $.

Since what's important is that $ A = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $ you can actually choose any 2 vectors which their outer product is $ {\sigma}_{1} {u}_{1} {v}_{1}^{T} $. The choice isn't unique. For example you could take for columns $ {\sigma}_{1}^{2} {u}_{1} $ and $ \frac{1}{{\sigma}_{1}} {v}_{1}^{T} $ for the rows.


Different ways of decomposing signals into sinusoids

Setting aside the issue of a Spectrogram or a sequence of windowed Fourier transforms.

A single windowed Fourier Transform can be shown to correspond to a Holder p norm minimization where p=2.

One can also use any p norm. so yes there are many possible Fourier types of compositions but the p=2 Holder norm can be bound in terms of the p=1 and p=infinity norms, and vice verse. This implies that these alternate spectra will not be radically different from the conventional DFT ( for discrete time).

A lot of sparse techniques depend on the zero norm.

The conventional Fourier Transform is a direct calculation while the others require iterative optimization.

In some cases like a sparse expansion, it is worth it but in general, the conventional Fourier transform is satisfactory. Also through Parsevals theorem, the conventional transform has a straightforward physical interpretation in terms of power.

Given the Wiener Khinchin theorem for random processes, the other norms would have seem to have less justification.

In a nutshell, yes you can have these other possible expansions but not without a lot of justification that the conventional Fourier series has.

There are cases where other expansions have been shown to be useful but not as a general approach


plotting phase transformations of a signal

My attempt to get the phase, see figure below:


Prewarping both resonant frequency $f_0$ and bandwidth (or $Q$) when using bilinear transform

alright, there are a couple of different issues that may (or may not) need to be de-conflated.

i'm gonna try to keep the number of symbols minimized.

$dB_\text{gain}$ is the number of dB gain of the peak (for $dB_\text{gain} > 0$) or cut (for $dB_\text{gain} < 0$).  it appears to be 6 dB in the plots.

$$A^2 \triangleq 10^{dB_\text{gain}/20}$$

is the linear gain for the boost or cut.

$f_\text{s}$ is sample rate.

$f_0$ is frequency of the boost or cut in the same units as the sample rate.

the analog transfer function of a resonant second-order filter (a.k.a. "biquad" or "SOS") is

$$\begin{align}
\\
H(s) &= \frac{b_0 + b_1 s + b_2 s^2 }{1 + a_1 s + a_2 s^2} \\
\\
 &= \frac{b_0 + b_1 s + b_2 s^2}{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\\
\end{align} $$

don't mistake the $a_k, \ b_k$ coefficients of the analog prototype with those in the resulting digital filter from the recipe in the cookbook. this is how the "EE definition" of $Q$ is defined.

an analog BPF with passband gain of 0 dB, has transfer function:

$$ H(s) = \frac{ \frac{1}{Q}\frac{s}{2 \pi f_0} }{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} $$

note when $f = f_0$ then $H(j 2 \pi f_0) = 1$.  note also that there are two bandedges $f_+$ and $f_-$ such that $f_- < f_0 < f_+$ and

$$ f_+ = f_0 \cdot 2^{bw/2} $$
$$ f_- = f_0 \cdot 2^{-bw/2} $$

and

$$ |H(j 2 \pi f_-)|^2 = |H(j 2 \pi f_+)|^2 = \frac{1}{2} \quad . $$

we define those bandedges for the BPF to be the "half-power frequencies" (a.k.a. "-3 dB frequencies").  and the $bw$ parameter is the bandwidth expressed in octaves.  the higher bandedge $f_+ = 2^{bw} f_-$ is $bw$ octaves higher log frequency than the lower bandedge $f_-$.  turns out that this bandwidth $bw$ in octaves is related to $Q$ as follows:

$$\begin{align}
\\
\frac{1}{Q} &= \frac{2^{bw} - 1}{2^{bw/2}} \\
\\
 &= 2 \ \sinh \left( \frac{\ln(2)}{2} \ bw \right) \\
\end{align} $$

keeping the same definition of $Q$, the "bell-shaped" boost/cut parametric EQ takes that BPF, gives it some gain (with sign) and adds it to a wire.  it has transfer function:

$$\begin{align}
H(s) &= (A^2 - 1)\frac{ \frac{1}{Q}\frac{s}{2 \pi f_0} }{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \ + 1 \\
\\
 &= \frac{1 + \frac{A^2}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}\\
\end{align}$$

this is the "traditional" analog parametric EQ.  note that
$$|H(j 2 \pi f_0)| = A^2 = 10^{dB_\text{gain}/20} \quad .$$
the problem is that, leaving $f_0$ and $Q$ constant, the curve for $dB_\text{gain} > 0$ is not a mirror image for $dB_\text{gain} < 0$, given the same number of dB, $|dB_\text{gain}|$, boost or cut.  the cut will be much skinnier than the boost.

some people might want the cut to exactly undo the boost given all other parameters being the same.  so in the cookbook (and in some other papers), a redefinition of $Q$ is made.  this redefinition is the substitution:

$$ Q \ \leftarrow \ A \cdot Q $$

so that makes the transfer function for the parametric EQ (with adjusted $Q$):

$$\begin{align}
 H(s) &= \frac{1 + \frac{A^2}{A \cdot Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{A \cdot Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\\
  &= \frac{1 + \frac{A}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{A \ Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\end{align} $$

note that if $dB_\text{gain}$ is replaced with $-dB_\text{gain}$, then $A$ is replaced with $\frac{1}{A}$ and the numerator and denominator of $H(s)$ are essentially swapped, which causes the frequency response of the cut to mirror that of the boost.

now, using the same relationship between bandwidth $bw$ and $Q$:

$$ \frac{1}{Q} = 2 \ \sinh \left( \frac{\ln(2)}{2} \ bw \right) $$

then the bandedges

$$ f_+ = f_0 \cdot 2^{+bw/2} $$
$$ f_- = f_0 \cdot 2^{-bw/2} $$

satisfy this gain definition:

$$ |H(j 2 \pi f_-)| = |H(j 2 \pi f_+)| =  A = 10^{dB_\text{gain}/40} $$

which are the "mid-gain frequencies" having dB gain of $\frac{dB_\text{gain}}{2}$.

so the definition of bandedge gain is a bit different, but at least for the analog prototype, we're keeping the relationship between $bw$ and $Q$ the same.  higher $Q$ means tighter $bw$.

so, first, before we discuss "warping Q", let's be completely consistent about which Q to compare.  for this, i might recommend leaving the "EE definition" of Q behind to not confuse.

so the bilinear transform which compensates for the frequency warping at the "significant frequency" or the resonant frequency $f_0$ makes this substitution:

$$ \text{normalized }s \triangleq \frac{s}{2 \pi f_0} \ \leftarrow \ \frac{1}{\tan(\pi f_0/f_\text{s})} \ \frac{1 - z^{-1}}{1 + z^{-1}} $$

$H(z)$ is the resulting digital filter transfer function after making that substitution for normalized $s$.

the frequency warping function turns out to be:

$$ f_\text{analog} = \frac{f_0}{\tan\left(\pi \frac{f_0}{f_\text{s}} \right)} \ \tan\left( \pi \frac{f_\text{digital}}{f_\text{s}} \right)  $$
or
$$ f_\text{digital} = \frac{f_\text{s}}{\pi} \ \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \frac{f_\text{analog}}{f_0} \right) $$

(note that DC gets mapped to DC and $f_0$ gets mapped to $f_0$, but when $f_\text{analog} \to \infty$, then $f_\text{digital} \to \frac{f_\text{s}}{2}$.)

fixing (or "prewarping") $f_0$ does not compensate for a cramped bandwidth because, besides the resonant frequency $f_0$, so also are the bandedges $f_+$ and $f_-$ warped by the bilinear transform.    measured in octaves the bandwidth, $BW$, in the digital filter is:

$$\begin{align}
 BW &= \log_2\left( \frac{f_\text{s}}{\pi} \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \right) \frac{f_+}{f_0} \right) - \log_2\left( \frac{f_\text{s}}{\pi} \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \frac{f_-}{f_0} \right) \right)  \\
\\
&= \log_2\left(\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{+bw/2} \right) \right) - \log_2\left(\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{-bw/2} \right) \right)  \\
\\
&= \log_2\left(\frac{\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{+bw/2} \right)}{\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{-bw/2} \right)} \right)  \\
\\
 \end{align}$$

you can see that the mapping of $bw$ to $BW$ is an odd-symmetry function, so it goes through $0$ and has no even-order terms in a Maclaurin (a.k.a. Taylor series) expansion.  now if you fix $f_0$ and plot digital $BW$ vs. analog $bw$, you will see that $BW < bw$ and that is the bandwidth cramping done by frequency warping of the bilinear transform.  to uncramp the bandwidth, you would have to solve for $bw$ in terms of $BW$ and $f_0$ and $f_\text{s}$.  and that is (how shall we say?) a female canine.

if you read the cookbook a little, you will notice a slightly adjusted mapping between $BW$ and $Q$:

$$ \frac{1}{Q} = 2 \ \sinh \left( \frac{\ln(2)}{2} \ BW \frac{2 \pi f_0/f_\text{s}}{\sin(2 \pi f_0/f_\text{s})} \right) $$

so the digital filter bandwidth $BW$ was increased by a factor of $\frac{2 \pi f_0/f_\text{s}}{\sin(2 \pi f_0/f_\text{s})}$ as a first-order attempt to compensate for the bandwidth cramping.  this result comes from assuming a narrow bandwidth in the first place and evaluating the derivative:

$$ \frac{\partial \, BW}{\partial \, bw} \Bigg|_{bw = 0} = \frac{\sin(2 \pi f_0/f_\text{s})}{2 \pi f_0/f_\text{s}} $$

the first term of the Maclaurin series is

$$ BW \approx \frac{\sin(2 \pi f_0/f_\text{s})}{2 \pi f_0/f_\text{s}} bw $$

compensating (or "prewarping") that first-order cramping of the bandwidth (which is what i think you mean by "prewarping Q") is done in the cookbook and always has been.  but it is only a first-order compensation.  if you want to do it better, be my guest and try to compute the next third-order term for the Maclaurin expansion (and flip it around so you have an expression of the approximate analog $bw$ in terms of digital $BW$).  doing this with Maclaurin series expansion is the only way i know because i don't think you will invert the $bw \to BW$ mapping above directly.


Using phase information to extract sample values

If the values are stored in a linear array indexed by $n$, read the values out
sequentially from a different starting point in the array, that is beginning with a nonzero value of $n$, and have the read out "wrap" around the end of the array, that
is, the incrementing of the index is done modulo $N$. This gives sinusoids
with possible initial
phases $0, 2\pi/N, (2)2\pi/N, (3)2\pi/N, \ldots, (N-1)2\pi/N$ depending on the
choice of starting point of the read out.  If you want other values
for the initial phase, you will need to do some computations, and not just
a plain reading out from memory.


Particle filter localization with different types of measures

I implemented many types of particle filter for one sensor. Since your case(i.e sensor fusion) includes multiple sensors, I try to make analogy using KF which I used many times for sensor fusion. In KF you perform time update and measurement update for each upcoming measurement sequentially. And if your measurements are asyncronous you modify your time update matrix each time. Thus, in particle filter it should be the same.

You, first, propagate the particles using time update equation then calculate the measurement likelihood for the corresponding sensor. You need to use timetag of each sensor measurment and likelihood equation should be modified depending on your sensor type.

The rest of the algorithm is the same. You perform resampling if necessary and move on.


Fourier Transforms and Series for the NON mathematically inclined.

The book "Blip, Ping & Buzz", by M. Denny (2007, John Hopkins University Press), explains signal processing, as used in sonar and radar, at the level of a supermarket science magazine.  Some algebra required.  More on the use of Fourier analysis than any theory.


Suggesting method for removing noise and image restoration?

The question is: "Suggest the method for removing noise with transforming $g$ and filter this image in the transformed domain." I believe it is more about "finding a method" than "doing the job practically".

Your first step is quite sound: remove the power complexity. But you are left with a product of a noise and a ($\log$) image. Most well known tools operate in a linear way. So use your hypotheses a little further: since $f>1$, its logarithm is positive. hence $v\log f$ is positive too. So you can turn the noise/signal product into a sum, by applying a logarithm again:

$$ \log \log g = \log v + \log \log f\,.$$

Now  $ \log v$ can be considered as an additive noise to $\log \log f$. With this linearized model, you can try get insights on the nature of the noise, looking at flat regions of the image, checking histograms, etc. Once done,
you can  apply your best filter (adapted to the statistics of $ \log v$, yet unknown), and $\exp \cdot \exp$ the result. The function $v\to \exp  \exp v$ grows very fast (from rechneronline.de):

So it relatively flattens low values with respect to higher values, with can result in  contrast unbalance. This effect could be counterbalanced by intensity bounding, nonlinear mappings, etc., but more importantly  by proper image restoration.

For your information, using linear filters in a $\log$ domain is an instance  of homomorphic filtering, which your solution may belong too:

  Homomorphic filtering is a generalized technique for signal and image
  processing, involving a nonlinear mapping to a different domain in
  which linear filter techniques are applied, followed by mapping back
  to the original domain


Matlab `filtfilt` provides excessive transient

It must have to do with the initial conditions used by the function filtfilt.m. The idea is to match initial conditions in a way such that startup and end transients are minimized. This, however, doesn't always seem to work, and it appears that for your filter specifications it actually does more harm than good. As far as I know there is no way to tell filtfilt.m not to mess with initial conditions. What remains in this case is to simply implement a straightforward version, just as you did.

Also take a look at this related answer.

EDIT:

I tried to reproduce your results with Octave, but I don't get any of the effects that you describe. I use a $10^{th}$ order Chebyshev 2 filter (which is the lowest order satisfying your specs), and I used a modified version of filtfilt allowing me to both switch on and off the use of initial conditions and the extrapolation of the input signal. There are two explanations for this: either your version of Matlab's filtfilt has a bug, or your filter is substantially different from mine and it introduces numerical problems (due to poles very close to the unit circle), which happen to manifest themselves when the input signal is extrapolated, and/or when special initial conditions are used (as is the case with filtfilt). If you provide your filter coefficients we can probably solve this mystery.

EDIT 2:

Tendero kindly reproduced the filter coefficients using the Matlab commands provided in the question:

b = [0.682585947577945 -0.012906667355225 -3.354515884719406 0.025370108447479 6.651816094305592 -0.000000000000000 -6.651816094305592 -0.025370108447479 3.354515884719406 0.012906667355225 -0.682585947577945];
a = [1.000000000000000 -0.194376467633997 -4.144930396202897 0.592472615834755 7.013126532356848 -0.691483043645839 -6.033921428993976 0.364443043103532 2.633924737223088 -0.072944768631229 -0.465923573748158];

I used these coefficients in Octave with the function filtfilt.m with a random input signal x = randn(2^20,1) just like in the question. As mentioned above, I used a modified version of filtfilt.m to be able to switch on and off the different measures for reducing transients:

no extrapolation of the input signal, zero initial conditions
no extrapolation of the input signal, optimal initial conditions (as implemented in the standard version of filtfilt.m)
extrapolation of the input signal, zero initial conditions
extrapolation of the input signal, optimal initial conditions

The last combination is the processing implemented in the standard version of filtfilt.m.

The plot below shows the first and last $50$ samples of the $4$ different output signals. For all other samples the signals are virtually identical.

As expected, one can see slight differences during startup and at the end, but these differences are very small, and the problem of the OP cannot be reproduced using the original filter coefficients and the Octave function filtfilt.m. There must be either a bug in the Matlab version of filtfilt.m used by the OP, or another problem that is neither related to the filter coefficients nor to the function filtfilt.m.


How to find the convolution kernel in frequency domain?

If you have noise present in your signal, the straightforward Fourier domain division will cause plenty of errors in your result. Some ways to avoid it are by using the so-called dual channel FFT (Part 1 and Part 2). I can also suggest deconvolution via adaptive filters, LMS or NLMS ([Normalized] Least Mean Squares) filters in particular are easy to understand are not hugely expensive in terms of CPU cycles in case your signals are long. LMS adaptive filters are very robust to noise.


I guess that your problem may not really be 1/f noise, but rather the effects from using the FFT directly on your time signal (without any windowing). If that is true, then try applying the Hann window to the time signal before the FFT and check the result.

There is plenty of information available about the window functions and how to use them.


Determining type of filter given its equation

The $\mathcal Z$-transform of the given difference equation is :

$$H(z) = 1 + z^{-2}$$

Put $z= \exp(j\omega)$ and remember that in discrete time systems the low frequencies are at $2n\pi$ $(n=0,1,2, \ldots)$ and the high frequencies occur at $(2n+1)\pi$; this is just a consequence of the periodic behavior of the discrete time complex exponentials.

So, $H(z) = 1 + \exp(-2j\omega)$ at $z = \exp(j\omega)$.

When $\omega=0; H(z) = 2$ and $w = \pi$ gives $H(z)= 2$.

Thus, both at high and low frequencies the the system function provides same
gain and hence the filter with the given $H(z)$ is a BAND REJECT/ NOTCH FILTER with $H(z) = 0$ at $\omega = \pi/2$.


Working with a sound's magnitude instead of amplitude

There is no reason why your piezo shouldn't be able to produce a bipolar output, if you use proper biasing and/or preamp. See for example https://www.homemade-circuits.com/diy-contact-mic-circuit/
Working with the magnitude only is probably a non-starter. $y = |x|$ is a highly non-linear operation and will dramatically change your spectrum. For example a sine wave at $f_0$ has just a single spectral line. The spectrum of the magnitude has a strong DC bias, lots of even harmonics $2f_0, 4f_0 ...$ but no content at $f_0$ at all anymore


That is simply converting analog frequency into digital frequency. For example, if you want a cutoff of 1kHz at a sample rate of 44.1kHz, the corresponding digital frequency is 2*pi*1000/44100 = 0.142... radiants.
In terms of DFT the digital frequency is the $ \omega $ in the $e^{^-j\cdot \omega }$ term of the DFT


How do I measure frequency of an audio file into an array?

Speech processing is a very long and detailed topic, but let's make a start.

Correlation

This probably won't work, but the first thing you can try is cross-correlation of the signals.

xc = xcorr(speaker_1, speaker_2);

Note that really all you may get out of this is the delay between each sample (if you're lucky).

You might also get some level of correlation coefficient by taking the peak value of

xc = xcorr(speaker_1/std(speaker_1), speaker_2/std(speaker2));

Where Next?

If that doesn't work (I'm not hopeful), you will need to look at lots of other options:

Blocking the data into frames
Finding features within these frames (LPCs, etc.).
Using some method to time-align these between speech samples (DTW)
Finding some measure of "likeness" between the time-aligned features.


I think wikipedia can more than adequately answer your question.

http://en.wikipedia.org/wiki/Phase_(waves)

However, in brief, the phase term describes the relationship between a waveform and a fixed reference point in time. For example the sinusoid $sin(\omega t)$ is zero at $t = 0$ whereas $sin(\omega t - \phi)$ is zero at $t = \phi$. $\phi$ could be referred to here as the phase offset.


How to Use the Walsh Hadamard Transform to Filter 24 bit RGB Color Values in a Localized Spatial Region?

The JPEG-XR Standard is using the Walsh Hadamard transform (Instead of DCT / Wavelets).

It has pretty good results (Unfortunately, it didn't get enough traction, yet it is supported on Windows).

You can follow the standard and have idea how to use it for image compression.
You'll be able to even get some source code.


It seems to me that we can directly calculate ... which can be done in a few second in matlab.

Who says that Matlab is calculating it directly, or that it isn't using Gershgorin's circle theorem in its algorithm?  Can you look at the code to see what they're doing?
Finding the eigenvalues of an $n \times n$ matrix is tantamount to factoring an $n^{th}$-degree polynomial.  Factoring any polynomial of degree greater than $n = 4$ is an iterative process, which works better if you can start with a good guess at the roots.
So, when Matlab does an eigenvalue decomposition of your $1000 \times 1000$ matrix, it is not "directly calculating" it -- it's guessing at values, then refining guesses, then reducing the matrix, etc., etc.
Pursuant to that, Gershogorin's circle theorem looks like a handy way to not only do get a set of initial guesses, but to bound the set of possible eigenvalues in a manner that's easy to test for if your algorithm starts to stray.

Is there anything wrong with the way I think?

I would say that believing that typing something into Matlab and getting an answer is the same thing as directly calculating it is an extreme error in thinking.  Someone knows how to solve the problem - you hope.  But if all you know how to do is press buttons in Matlab and get an answer, the set of persons who really know how to solve the problem doesn't include you.


When should I calculate PSD instead of plain FFT magnitude spectrum?

Now, I would like to show what frequencies the speech has. However,
  I'm not sure what would be the best way to do that. It seems sometimes
  one calculates the absolute value of a Fourier transform, and
  sometimes power spectral density.

If you want to attach physical meaning to your analysis, then go with the power spectral density, (PSD). This is because this will simply give you the power of your signal, in each frequency band. On the other hand if you do not want/care about a physical meaning, but want to know how the fourier amplitudes of each band vary relative to each other, you can stick to absolute magnitude.

In practice, you can compute the PSD as simply the absolute magnitude of the fourier transform squared. For example, if your signal is $x[n]$, and its DFT is $X(f)$, then the absolute magnitude of the DFT is $|X(f)|$, while the PSD is $|X(f)|^2$.

  If I understand correctly, the latter works so that I divide my signal
  into parts, do FFT part-by-part and somehow sum these. Window
  functions are somehow involved. Can you clarify this a bit for me? I'm
  new to DSP.

No, this is not true. What you are talking about here refers to the Short Time Fourier Transform, (STFT). This is simply chopping up your time domain signal, widowing it, and then taking the fourier trnasform. At the end of the day though, you will still have a complex matrix. If you choose to take its absolute magnitude, you will have an absolute magnitude fourier transform matrix. If you take its absolute magnitude squared, you will have a power spectral density matrix.