Can I determine a system's $z$-domain transfer function from its pole-zero plot?

Without knowing specifics (ignoring a proportionality constant we'll call K as inidicated in the more compete answer above).  Note that K does not affect the dynamic behavior of the system, so it may not be of particular interest.

it goes like this...

H(z) is a rational function of the form
$$H(z) = \frac{(z-Z1)(z-Z2)...(z-Zm)}{(z-P1)(z-P2)...(z-Pn)}$$

where Z1,Z2...Zm are your zero values and P1,P2...Pn are your pole values.

I assume you know the exact values of you're poles and zeros. If you don't, you'll have to determine approximations from the PZP.


How to solve this equation (about complex envelope)?

If $e(t)$ is band-limited and if the frequency $\omega_0$ is chosen such that the complex signal

$$y(t)=z(t)e^{j\omega_0 t}\tag{1}$$

has no negative frequencies, then $y(t)$ is called the analytic signal, and it can be determined from $e(t)$ using the Hilbert transform:

$$y(t)=e(t)+j\mathcal{H}\{e(t)\}\tag{2}$$

This analytic signal $y(t)$ is unique. However, the complex envelope $z(t)$ is not, unless we define the frequency $\omega_0$. Even though $\omega_0$ can be called center frequency, there is no unique definition of such a center frequency. That's why $\omega_0$ is often called reference frequency, and the complex envelope $z(t)$ depends on the choice of that reference frequency.

If $\omega_0$ is not given, you cannot uniquely solve for $z(t)$. It's up to you to define a convenient value for $\omega_0$, from which the complex envelope is easily obtained:

$$z(t)=y(t)e^{-j\omega_0 t}\tag{3}$$

Usually you would choose $\omega_0$ such that $z(t)$ is a low pass signal.

Note that since

$$e(t)=\Re\{z(t)e^{j\omega_0 t}\}=\frac12\left\{z(t)e^{j\omega_0 t}+z^*(t)e^{-j\omega_0 t}\right\}\tag{4}$$

(where $*$ denotes complex conjugation), the complex envelope can also be determined form $e(t)$ by down conversion and low pass filtering:

$$2e(t)e^{-j\omega_0 t}=z(t)+z^*(t)e^{-2j\omega_0 t}\tag{5}$$

Assuming that $z(t)$ is a low pass signal and that $\omega_0$ is chosen such that it is greater than the bandwidth of $z(t)$, low pass filtering $(5)$ will yield $z(t)$.


What are some of the differences between DFT and FFT that make FFT so fast?

The naive implementation of an $N$-point DFT is basically a multiplication by a $N \times N$ matrix. This results in a complexity of $\mathcal{O}(N^2)$.

One of the most common Fast Fourier Transform (FFT) algorithm is the radix-2 Cooley-Tukey Decimation-in-Time FFT algorithm. This is a basic divide and conquer approach.

First define the "twiddle factor" as:
$$W_N \triangleq e^{-j\frac{2\pi}{N}}$$
where $j \triangleq \sqrt{-1}$ is the imaginary unit, then the DFT $X[k]$ of $x[n]$ is given by
$$X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{kn} \, .$$
If $N$ is even (and $\tfrac{N}{2}$ is an integer), the sum can then be divided in two sums as follows
$$X[k] = \sum_{n=0}^{N/2-1} x[2n]W_N^{2kn} + \sum_{n=0}^{N/2-1} x[2n+1]W_N^{k(2n+1)}$$
where the first summation deals with the even samples of $x[n]$ and the second with the odd samples of $x[n]$. Defining $x_e[n] \triangleq x[2n]$ and $x_o[n] \triangleq x[2n+1]$ and using the fact that

$W_N^{k(2n+1)} = W_N^{2kn}W_N^k$, and
$W_N^{2kn} = W_{N/2}^{kn}$,

this can be re-written as
$$
\begin{align}
X[k] &= \sum_{n=0}^{N/2-1} x_e[n] W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1} x_o[n]W_{N/2}^{kn} \\
& = X_e[k] + W_N^k X_o[k]
\end{align}
$$
where $X_e[k]$ and $X_o[k]$ are the $\tfrac{N}{2}$-point DFT transforms of the even and odd samples of $x[n]$ respectively. So we just transformed a single $N$-point DFT into two smaller $\tfrac{N}{2}$-point DFTs.  This reduces computational cost because
$$ 2 \left( \frac{N}{2} \right)^2 + N < N^2 $$
when $N > 2$.

We can then re-iterate the same process on those two smaller DFTs. This divide-and-conquer approach allows to reach complexity of $\mathcal{O}(N\log N)$, which is way better than the $\mathcal{O}(N^2)$ we had with the naive DFT implementation (as is greatly illustrated by leftaroundabout's answer).


What is the difference between cubic interpolation and cubic "Spline" interpolation?. How to use it for upsampling purpose?

The difference between cubic interpolation as described in your question and cubic spline interpolation is that in cubic interpolation you use 4 data points to compute the polynomial. There are no constraints on the derivatives. Cubic spline interpolation computes a third order polynomial only from two data points with the additional constraint that the first and second derivative at the interpolation points are continuous. So if you have 4 points, then you compute 3 different polynomials (between points 1-2, 2-3, and 3-4), and these polynomials are smoothly connected in the sense that their first and second derivatives are equal at the given data points. This is not the case with standard polynomial interpolation.

But if you want an efficient implementation I wouldn't necessarily use polynomial or spline interpolation but I would go for a standard polyphase implementation of a low-pass interpolation filter. If $L$ is the upsampling factor you get $L$ parallel polyphase filters operating at the lower sampling rate, followed by $L$ upsamplers, the results of which are added. The low-pass filter can be designed using the window method (e.g. a Kaiser window). There is a lot of material available on the polyphase implementation of interpolators. This Master's thesis is very accessible. Just have a look at single-stage polyphase implementations using linear-phase FIR filters and ignore the rest.

EDIT:

If you want or need to use cubic interpolation, here's what you need to do. First of all, some coefficients in your equations are wrong, that's why also the result is not correct. These would be the correct equations:

    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 27d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 25c + 125d

However, you don't want to do it this way because if you always use the sample index as the $x$-coordinate you'll end up with a very large range for the coefficients of the system of linear equations (because you need to take the third power), and your coefficient matrix will sooner or later become ill-conditioned (asking for numerical problems). Second, and even more importantly, you would get a different system of equations for every new segment of the input signal, which would be very inefficient. You can do the interpolation by solving a $4\times 4$ system of linear equations with the same coefficient matrix each time. Just the left-hand side vector (containing the values of the given data points, V in your example) will be different for each segment.

You compute a cubic polynomial for each segment (i.e. for each range between two given data points) by considering the data points defining the segment and the two adjacent data points, just as in your example. If you define the cubic polynomial as

$$P(x)=a_0+a_1x+a_2x^2+a_3x^3$$

and if you define the $x$-axis values of the 4 current input data points as $-1$, $0$, $1$, and $2$, you get the following system of equations:

$$\begin{align}y_0&=a_0-a_1+a_2-a_3\\
y_1&=a_0\\
y_2&=a_0+a_1+a_2+a_3\\
y_3&=a_0+2a_1+4a_2+8a_3\end{align}$$

where $y_i$, $k=0,1,\ldots, 3$ are the current data points, and the current segment is between the points $y_1$ and $y_2$. The system above is readily solved as

$$\begin{align}
a_0&=y_1\\
a_1&=-\frac13 y_0 - \frac12 y_1 + y_2 - \frac16 y_3\\
a_2&=\frac12  (y_0 + y_2) - y_1\\
a_3&=\frac12  \left[(y_1 - y_2) + \frac13 (y_3 - y_0) \right]\\
\end{align}$$

So you just need to evaluate the above 4 expressions for the polynomial coefficients for each set of 4 data values. As soon as you have the coefficients, you can compute the interpolated output values in the current segment by evaluating the polynomial at $x_i=i/L$, $i=1,2,\ldots,L-1$, where $L$ is the upsampling factor ($L=10$ in your case). Then you forget the oldest input value, read the new input value, and you have a new set of 4 values from which you compute a new polynomial for the current segment, etc.

Here is a little Matlab/Octave script implementing this method:

x = randn(20,1);        % some input data
L = 10;                 % upsampling factor

% append points at both ends for first and last segment
N = length(x);
x2 = [x(1);x(:);x(N)];
Nseg = N-1;         % # segments

% matrix to compute polynomial coefficients from the current 4 data points: a = A*y
A = [0 1 0 0;-1/3 -.5 1 -1/6;.5 -1 .5 0;-1/6 .5 -.5 1/6];

% x-axis values of interpolated values in the current segment
mu = (0:L-1)/L;

% allocate output vector
Ny = Nseg*L + 1;
y = zeros(Ny,1);

% loop over all segments
for i = 1:Nseg,
    xt = x2(i:i+3);
    a = A*xt;
    yt = a(1) + a(2)*mu + a(3)*mu.^2 + a(4)*mu.^3;
    Iy = (i-1)*L + (1:L);
    y(Iy) = yt;
end

y(Ny) = x(N);

% plot result
tx = 0:N-1;
ty = (0:Ny-1)*(N-1)/(Ny-1);
plot(ty,y,'r-'), hold on
stem(ty,y,'r'), hold on
stem(tx,x), hold off


Are circular buffers relevant to transposed direct form IIR filters?

Circular buffers are only useful if the contents of the delay elements are not changed when shifted. So you're right that circular buffers cannot be used for the transposed structures. However, this doesn't necessarily mean that the transposed structures are less efficient than the direct forms, it just means that the data aren't just shifted through a delay line. The complexity for computing the new state vector in the transposed structure is comparable to the complexity for computing the output from the states in the direct form.


What is the type number of a discrete time system given $H(z)$?

I couldn't find a definition of system type for discrete-time control systems, but in analogy with continuous-time systems you should count the number of poles at $z=1$, not at $z=0$, because in discrete time, $z=1$ is the DC point, just like $s=0$ in continuous time.

A discrete-time integrator (an accumulator), has a pole at $z=1$. Its transfer function is given by

$$H(z)=\frac{z}{z-1}\tag{1}$$


Are there any efficient algorithm to detect circular shape of arbitrary size in an image?

Usually this problem is tackled using the Hough transform. Check it here:

OpenCV
Scikit

After you have some hints for the possible location of the circles and their diameter then you can calculate the distance (e.g. Euclidean) from the pixels you have to an ideal model (circle parametric equations: x = r.cos a; y = r.cos a)


How to calculate PSD if i have power vs freq spectrum?

i used this formula saw on wiki, https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density
it works okay for me,
thanks for the comments, it helped me to understand the data what I was receiving.


Fail to separate sound signals by FastICA on real-world recording

ICA in raw form is only suitable for use with phase synchronised observation mixtures. Using microphones as you have described will introduce a phase delay as pointed out by other posters. However this phase delay can be used to great avail. The best known algorithm that deals with stereo separation in the presence of delays is DUET. The links are broken but the references you are looking for are here >http://eleceng.ucd.ie/~srickard/bss.html.

This is the paper you should look for >
A. Jourjine, S. Rickard, and O. Yilmaz, Blind Separation of Disjoint Orthogonal Signals: Demixing N Sources from 2 Mixtures, IEEE Conference on Acoustics, Speech, and Signal Processing (ICASSP2000), Volume 5, Pages 2985-2988, Istanbul, Turkey, June 2000


assuming it's a low-pass filter and assuming that the filter order equals the number of poles and there are no zeros in the s-domain transfer function, then the rolloff at the edge of the passband through the stopband is -20 dB/decade (or -6 dB/octave) multiplied by the order of the filter.
now Butterworth and Tchebychev Type 1 filters have no zeros.
Tchebychev Type 2 and Elliptical filters have zeros, so the ultimate rolloff is not there, the gain remains at the stopband level.


The mechanical stresses get amplified by the vibration motor, giving me erroneous readings.

Nothing to do with amplification. Basic high school-level physics!

"Weight" is just the force that an object exerts due to gravity acting on it. Your load cells hence simply measure the force pressing down on them, and, that's the whole idea of vibration motors, that is the sum of the gravity pull and the acceleration applied to the vibrating mass (see: Newton's laws: when your vibration motor moves some mass upwards, it's getting pushed down by a force identical to the lifting force, and vice versa).

  Are there any filtering techniques I could use in order to be able to get better real time data from the load cells?

Any kind of sufficiently low-cutoff low pass filter will do, because no matter how much you vibrate, the mass of the system doesn't change, and what goes up needs to come down.

You'd typically either simply do this with a mechanically dampening system, or with a simple RC low pass operating on your analog electrical signal, which you need for anti-aliasing anyway, because otherwise you're exposing your digital signal to aliases.


narrow band frequency response creation and general FR question

Let's assume that speaker, free space and microphone can be modeled as linear systems and that the overall system transfer function is given by $H(\omega)$. Then the received signal is given by
$$
B(\omega) = A(\omega)H(\omega),
$$
where $A(\omega)$ is the signal sent to the speaker. So, yes, if you know the transfer function $H(\omega)$ you can calculate $B(\omega)$ and it should align with what you record with the microphone.

To determine $H(\omega)$ in practice you will have use discrete-time signals and play a unit impulse $a(n) = \delta(n)$, record the received signal $b(n)$ and calculate the transfer function with help of the discrete Fourier transform (DFT):
$$
B(k) = DFT[b(n)]
$$
Now the maximum frequency $\omega_m$ you are able to analyze depends on the sampling rate $\omega_s$ with which the samples of $b(t)$ are acquired after the microphone to obtain $b(n)$: $\omega_m = \omega_s/2$. Frequencies greater than $\omega_m$ can neither be measured nor calculated and will influence your result because of the Aliasing effect unless they're filtered with a low pass filter before sampling. An alternative approach is to transmit a signal $A(k)$ that only contains frequencies in the range you're interested. The transfer function for this limited frequency range can then be calculated by
$$
H(k) = \frac{B(k)}{A(k)}
$$
In this case frequencies lying outside this range won't affect your result as long as it's a linear system.

A comment on processing power: the sampling rate influences the clock rate in a real-time system and doesn't affect complexity otherwise. The cost of the DFT operation (implemented by the FFT algorithm) increases with the length $N$ of the recorded signal $b(n)$. The greater $N$ the greater the frequency resolution of $B(k)$.


Implementing a Band-Pass Filter from Scratch

A simple way to make a real bandpass filter (assuming you know how to make a low pass filter) is multiply the coefficients of your prototype low pass filter by $cos(2\pi f_c nT)$ where $f_c$ is the desired center frequency of your passband and T is the sampling interval for your coefficients.


The problem is in the input signals to your polyphase filters. They should be delayed and subsampled versions of the input signal, like in the following code fragment:

M=8;    % downsampling factor
L=256;  % length of input signal, integer multiple of M
x=randn(L,1);  % input signal
h = fir1(127, 1/M);
yp = zeros(L/M, 1);

for i = 1:M,
    xtmp = [zeros(i-1,1); x];  % delayed input signal
    tmp = filter( h(i:M:end), 1, xtmp(1:M:L) );  % polyphase filtering
    yp = yp + tmp;  % accumulate outputs of polyphase filters
end

% compare with standard filtering
y = filter ( h, 1, x);
yM = downsample( y, M );

n=1:L/M;
plot(n,(yM-yp)/max(yM))


Tutorial on 2nd generation wavelets (with lifting)?

I've bought finally a copy of [Ripples in Mathematics
The Discrete Wavelet Transform][1], and I'm very pleased by this book. The authors explain the DWT with alternating point of views (lifting schemes, filter banks approach, multi-resolution analysis), where each of these viewpoints has its own advantages. Furthermore, the book is implementation oriented, with chapters about boundaries handling and matlab/C implementations.

I'm still looking on a proper way to handle odd-sized signals, but Ripples gave me a good start.

[1]: http://www.control.auc.dk/~alc/ripples.html "Ripples in Mathematics The Discrete Wavelet Transform", by Arne Jensen and Anders la Cour-Harbo


How can I call the area under an impulse/spike?

What you're describing "literally" is just the area under a curve, the integral. I don't think this helps you.
What you want, functionally, is not quite clear, but I'd bet it's either

the energy of the pulse or
the energy-bandwidth product.

The energy of the signal is the integral over the (magnitude)
squared signal in the continuous time case, or the sum over the
(magnitude) squared signal in continuous time.
The bandwidth is a frequency domain concept, and will tend to be larger the "shorter and sharper" the pulse is. The product of bandwidth and energy has important estimation quality effects for example in radar applications.


You may just implements a VAD (voice activity detector) whose parameters are high enough to filter out voiced consonants. One VAD which always amazed me is the zero crossing rate, because it is so simple to implements yet so difficult to model mathematicaly.


Why does the addition of noise improve coherence between two signals?

To elaborate on @user28715's answer:
Even with a minuscule amount of noise say An = 0.0000001 you will get a clean graph like your first graph. Python's implementation of the Welch method certainly uses finite-precision floating point numbers and quantization errors in the input and from the calculation steps such as windowing and Fast Fourier Transform (FFT) propagate to the results. From SciPy docs:

Cxy = abs(Pxy)**2/(Pxx*Pyy), where Pxx and Pyy are power spectral
density estimates of X and Y, and Pxy is the cross spectral density
estimate of X and Y.

The normalizing division in the equation means that the values of Pxx and Pyy can be arbitrarily close to 0 while still giving large Cxy.


How do you find signal similarity when one signal has zero values? Correlation alternative?

The Correlator is really just a Maximum Likelihood detector. If your noise is AWGN, (Additive White Gaussian Noise), then the correlator is optimal in the Least-squares sense. This is because a correlator is doing a match filtering with your incoming signal, and in this sense, will only consider DFT bands where SNR is the highest, or rather, only consider DFT bands where your signal energy exists.

  "This made me wonder if there's an alternative to cross-correlation
  which would work properly for a zero valued signal."

Regarding the second paragraph, the 'signal' you are trying to 'correlate' against, [0 0 0 0 0], has zero energy, so the above would not apply. It does not have to do with its mean (DC band of the DFT is null), it has to do with the fact that it has no energy. (All DFT bands are null). In a special case such as this, then you could perform:

$$
e[n] = \sum_{n=0}^{N-1} |x[n] - y[n]|
$$

where $x[n]$ is your known signal, and $y[n]$ is your given signal. In effect, this is nothing but the L1-norm of the difference of vectors. You would then pick the $y[n]$ vector that minimizes this error.

(Note that if your 'known signal' is always the zero vector, then this collapses to simply the L1-norm of your given vector. The one with the lowest L1-norm is the one closest to the zero vector).

In contrast with the matched filter, this metric will not suffer if your 'signal' has no energy.


Given your seed locations and the threshold of 1, that means that not every pixel can be labelled. See below.
$$
\begin{array}{ccccccc}
5 & 6 & 1 & 1 & 1 & 0 & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}}\\
6 & 7 & 7 & 1 & 1 & 1 & 6\\
3 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & 7 & 6 & 5\\
4 & 3 & 7 & 7& 7 & 5 & 3\\
3 & 3 & 2 & 1 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 3 & 4\\
2 & 2 & 1 & 2 & 1 & 2 & 3\\
1 & 1 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 1 & 1 & 2 & 2\\
\end{array}\\
\Huge \updownarrow\\
\begin{array}{ccccccc}
5 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} &\enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}} & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}}\\
\enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="green"]{\color{black}{6}}\\
\enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & 7 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 5\\
4 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 7 & 7& 7 & 5 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}}\\
\enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 4\\
2 & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}}\\
\enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & 2\\
\end{array}\\
$$
It looks like there should be 4 labels, not three... and, yes, it looks like the 0 and 1 seeded regions have been merged.

In addition to the lack of labels for all pixels (3 vs 4), I wonder if it's because the calculations are being done in (double precision) floating point so that $|1.00 - 0.00| < 1.00$


Your definition looks fine to me, if your goal is to treat the signal at every fixed $(x,y)$ location as a function of time and then compute 1D total variation along time. Is that what you intend to do?

Here are a couple of other ways you can define it (again, I have no clue what your application is, and why you want to compute the TV norm in the first place, so I don't know which of these make sense in your application).

To fix notation, let's denote your signal of interest as $z(n_x,n_y,n_t)$ where $z$ is the quantity measured at discrete locations indexed by $(n_x,n_y)$ in space and at a time sample number $n_t$.

TV norm for each fixed spatial location:
Fix a location $(n_x^*,n_y^*)$.
$$
TV = \sum_{n_t} |z(n_x^*,n_y^*,n_t+1) - z(n_x^*,n_y^*,n_t)|
$$

TV norm for each 2D time snapshot:
Fix a time instant $n^*_t$.
$$
TV = \sum_{n_x,n_y} |z(n_x+1,n_y,n^*_t) - z(n_x,n_y,n^*_t)| + | z(n_x,n_y+1,n^*_t) - z(n_x,n_y,n^*_t)|
$$

TV norm in 3D space-time:

$$
TV = \sum_{n_x,n_y,n_t} |z(n_x+1,n_y,n_t) - z(n_x,n_y,n_t)| + | z(n_x,n_y+1,n_t) - z(n_x,n_y,n_t)| + |z(n_x,n_y,n_t+1) - z(n_x,n_y,n_t)|
$$

See also: https://en.wikipedia.org/wiki/Total_variation_denoising#2D_signal_images


It can have a positive gain, it all depends on your definition of 0dB. Sometimes, the first path is assumed to be 0dB, then it can happen that later paths have a higher gain than the first one, leading to positive path gain.

However, in reality, there is always the path loss between the transmitter and the receiver, which is normally added to all paths similarly.

So, for modeling and simulation, it is possible to have positive path gain. However, in reality the signal at the receiver will always have less power than the transmitted signal.


Creating a 4th order Butterworth filter R

do I need to change the number to “2” so that it is a second-order

Yes

and should I also double the cut-off frequency i.e change “W” to 0.3 from 0.15.

No
Strictly speaking, doubling up the filter will make the gain at the cutoff frequency -6dB instead of -3dB. So you have created a Linkwitz-Riley and not a Butterworth filter.
If you want a gain of exactly -3dB at 150 Hz you would need to increase the frequency slightly. You can work out the math by calculating the -1.5dB point of a Butterworth filter, but trial and error is probably quicker. 185Hz should get you fairly close.
In any case, the shape will still NOT be that of a Butterworth. Cascading two Butterworth filters does not result in another Butterworth filter. The zeros are in the right location but the poles are simply duplicated and are not Butterworth-distributed any more.


Why are Gaussian filters used as low pass filters in image processing?

Image processing applications are different from say audio processing applications, because many of them are tuned for the eye. Gaussian masks nearly perfectly simulate optical blur (see also point spread functions). In any image processing application oriented at artistic production, Gaussian filters are used for blurring by default.

Another important quantitative property of Gaussian filters is that they're everywhere non-negative. This is important because most 1D signals vary about 0 ($x \in  \mathbb{R}$) and can have either positive or negative values. Images are different in the sense that all values of an image are non-negative ($x \in  \mathbb{R}^+$). Convolution with a Gaussian kernel (filter) guarantees a non-negative result, so such function maps non-negative values to other non-negative values ($f: \mathbb{R}^+ \rightarrow \mathbb{R}^+$). The result is therefore always another valid image.

In general, frequency rejection in Image processing in not as crucial as in 1D signals. For example, in modulation schemes your filters need to be very precise to reject other channels transmitted on different carrier frequencies, and so on. I can't think of anything just as constraining for image processing problems.


Improving the communication speed on a 10km cable (with high capacitance)

It seems your hardware is limited. If you can't detect more finite frequencies and you can't implement another method of modulation you might have reached your maximum in physical bits with your hardware.

May I suggest a bit of coding theory? Implement LZW compression. The data rate will vary and it might not give you double the speed but will definitely increase the data rate.

Its already produced in c, try it out.

http://marknelson.us/attachments/1989/lzw-data-compression/lzw.c


SNR stands for Signal to Noise Ratio. It is a ratio and as such does not have any units, it describes the proportion of signal to undesired noise. There is no single correct measure of SNR, it differs depending on the application.

In the equation you have given, the SNR is broken down in the following way:

1) Calculate the power ratio

$$\frac{\sigma_s^2}{\sigma_n^2}$$

Your equation requires us to know the variance of the signal and the noise, dividing one by the other gives the power ratio. To answer one of your questions, this is the ratio of power, not amplitude, we could achieve the same result by taking a ratio of the root mean square of the signal and noise - another measure of power.

2) Express the ratio in decibels

$$10\log_{10}$$

This is simply a convenience, the SNR may be a very large number so expressing the value on a logarithmic scale can make the information more manageable. Converting a unit to decibels simply involves taking the log10 of the value and multiplying by 10.


Well, if the goal is to learn that some recording was similar to some reference, e.g. that the recording sound like the word "waffle", then yes, there's a way to reconstruct things, and that would be speech synthesis of "waffle", for example.

That's not a great example, but it's meant to illustrate one thing:

If your reduced version of the original signal, i.e. your feature vector (MFCCs in your case) is good enough to tell different things apart, and you can restrict the content of the signal to things that are actually well-represented by that feature set (e.g. you say "only voice, not music"), then yes, you can basically use closest-matches from a reference signal database.


What is the official definition of a feature in signal classification?

Features are properties of the signal used to describe it in some way. Sometimes the term "descriptors" is used, for example MPEG-7 audio-visual description standard:

Time domain features are computed from the audio waveform, while frequency-domain or spectral features are derived after performing an FFT.

You should look through a feature-extraction software package like Yaafe, to get some practical intuition. There are also many tutorials and papers online for further learning.


What Does the Total Variation Norm Mean in the Context of Image Processing

The Total Variation of an image $ I $ can be calculated in one of 3 methods (See The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework):

Anisotropic TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_1 $.
Isotropic TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_2 $.
Isotropic Squared TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_2^2 $.

In practice, all yield almost the same result.
The Total Variation is a measure of how close the image is to a Piece Wise Constant function. Namely, How sparse are the changes in the gradient of the image.
You may read more at Total Variation Denoising.


Given a set of lines, find only those who are parallel (perspective)

Parallel lines in the image do intersect at a vanishing point. Therefore simply hypothesizing lines (a gradient direction at a point suffices to describe it) and voting (see Hough voting) would suffice to identify this point. One could then record all the lines that casted votes to this very point and identify them. Care must be taken as it is difficult to vote directly on a vanishing point, i.e. the space is not bounded and explodes to infinity. Instead, it is possible to use a Gaussian Sphere as Barnard suggested.

Vanishing point extraction is a well studied topic and many algorithms exist, some of which are:

  A Novel Line Space Voting Method for Vanishing-Point Detection of
  General Road Images, Zongsheng Wu, Weiping Fu, Ru Xue and Wen Wang
  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4970002/

  Vanishing point detection for road detection, Hui Kong, Jean-Yves
  Audibert, Jean Ponce http://www.di.ens.fr/sierra/pdfs/cvpr09c.pdf

  Vanishing point detection in corridors: using hough transform and
  K-means clustering, R. Ebrahimpour,  R. Rasoolinezhad, Z.
  Hajiabolhasani,  M. Ebrahimi
  http://ieeexplore.ieee.org/document/6135447/

  Hough Transform for Directional Orientation, Ivan Papusha, Matthew
  Ho
  https://stacks.stanford.edu/file/druid:rz261ds9725/Ho_Papusha_RealTimeHoughTransform.pdf

Finally, when multiple lines cause a problem, one could always cluster them using the line parameters $[a,b,c]$. When number of lines are not known, one can use a hierarchical clustering, such as single-linkage or agglomerative clustering.


Real and imaginary part of signal after fft

Here is one way to do it properly:

import numpy as np

#=======================================================================
def main():

        t = np.arange( 0, 200 )
        N = len( t )

        signal1 = 0.1 * np.sin( .1 * t )
        signal2 = 0.4 * np.sin( .1 * t + .35 )

        rms_time1 = np.sqrt( np.sum( signal1 * signal1 ) / N )
        rms_time2 = np.sqrt( np.sum( signal2 * signal2 ) / N )
        q_time = rms_time1 / rms_time2

        print q_time, rms_time1, rms_time2

        dft1 = np.fft.fft( signal1 )
        dft2 = np.fft.fft( signal2 )

        dft1_real = np.real( dft1 )
        dft1_imag = np.imag( dft1 )
        dft2_real = np.real( dft2 )
        dft2_imag = np.imag( dft2 )

        rms_freq_real1 = np.sqrt( np.sum( dft1_real * dft1_real ) / N )
        rms_freq_imag1 = np.sqrt( np.sum( dft1_imag * dft1_imag ) / N )
        rms_freq_real2 = np.sqrt( np.sum( dft2_real * dft2_real ) / N )
        rms_freq_imag2 = np.sqrt( np.sum( dft2_imag * dft2_imag ) / N )

        rms_freq1 = np.sqrt( rms_freq_real1 * rms_freq_real1 \
                           + rms_freq_imag1 * rms_freq_imag1 )

        rms_freq2 = np.sqrt( rms_freq_real2 * rms_freq_real2 \
                           + rms_freq_imag2 * rms_freq_imag2 )

        q_freq = rms_freq1 / rms_freq2

        print q_freq, rms_freq1, rms_freq2

        print rms_freq1 / rms_time1
        print rms_freq2 / rms_time2
        print np.sqrt( 200 )

#=======================================================================
main()

The output is:

0.246138603222 0.0699023232816 0.283995774603
0.246138603222 0.988568136262 4.016306761
14.1421356237
14.1421356237
14.1421356237

Followup:

This is a more straightforward way to calculate the RMS of the DFT bins.  I coded the sample above to more closely align with how I thought the OP was doing it.

        sumsquares1 = np.real( np.sum( dft1 * dft1.conjugate() ) )
        sumsquares2 = np.real( np.sum( dft2 * dft2.conjugate() ) )

        rms_freq1 = np.sqrt( sumsquares1 / N )
        rms_freq2 = np.sqrt( sumsquares2 / N )

        q_freq = rms_freq1 / rms_freq2

        print q_freq, rms_freq1, rms_freq2


How is voice sound on TV or phone distorted through a face mask?

There are many researches on the acoustic effect of masks since the pandemic. You can find them on the website of The Journal of the Acoustical Society of America and most of them are free to access.
Masks have selectivity on both frequency and spatial.
In frequency domain, masks react as low-pass filters that attenuate high frequencies significantly. Interestingly, different research gives different results. They are saying 1 kHz$^{[1]}$, 2 kHz$^{[2]}$ and 3 kHz$^{[3]}$ maybe due to different understanding about "significantly".
In spatial domain, masks lead to reduced frontal sound radiation compared to the lateral or rear radiation at high frequencies$^{[2]}$.
Researchers also tested the effect of microphone placement. They placed microphones on the lapel, cheek, forehead, and next to the mouth$^{[1]}$. These differences are mostly caused by the spatial radiation pattern with a mask.
And of course, different masks have different acoustic effect according to their materials, weaves and types.
References

[1] Acoustic effects of medical, cloth, and transparent face masks on speech signals
[2] Impact of face masks on voice radiation
[3] Effects of face masks on acoustic analysis and speech perception: Implications for peri-pandemic protocols
[4] Acoustic effects of non-transparent and transparent face coverings

[5] Effect of masks on speech intelligibility in auralized classrooms


Pitch Detection / Harmonic Product Spectrum - What's wrong?

The Harmonic Product Spectrum pitch estimation method, using pure multiplication weighting, only works well for signals that have a full set of harmonics with sufficient magnitudes.  Before using HPS, one should either validate that the signal in question meets that criteria, or add a suitably large non-zero floor to all harmonic magnitudes being multiplied together.  Otherwise a pure sinewave with zero harmonics will simply disappear from the result after being multiplied by a zero magnitude.

Guitar and ukulele sounds usually have a far richer set of harmonics than test sine waves.


Applying filter in scipy.signal: Use lfilter or filtfilt?

filtfilt is zero-phase filtering, which doesn't shift the signal as it filters.  Since the phase is zero at all frequencies, it is also linear-phase.  Filtering backwards in time requires you to predict the future, so it can't be used in "online" real-life applications, only for offline processing of recordings of signals.
lfilter is causal forward-in-time filtering only, similar to a real-life electronic filter.  It can't be zero-phase.  It can be linear-phase (symmetrical FIR), but usually isn't.  Usually it adds different amounts of delay at different frequencies.

An example and image should make it obvious.  Although the magnitude of the frequency response of the filters is identical (top left and top right), the zero-phase lowpass lines up with the original signal, just without high frequency content, while the minimum phase filtering delays the signal in a causal way:

from __future__ import division, print_function
import numpy as np
from numpy.random import randn
from numpy.fft import rfft
from scipy import signal
import matplotlib.pyplot as plt

b, a = signal.butter(4, 0.03, analog=False)

# Show that frequency response is the same
impulse = np.zeros(1000)
impulse[500] = 1

# Applies filter forward and backward in time
imp_ff = signal.filtfilt(b, a, impulse)

# Applies filter forward in time twice (for same frequency response)
imp_lf = signal.lfilter(b, a, signal.lfilter(b, a, impulse))

plt.subplot(2, 2, 1)
plt.semilogx(20*np.log10(np.abs(rfft(imp_lf))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('lfilter')

plt.subplot(2, 2, 2)
plt.semilogx(20*np.log10(np.abs(rfft(imp_ff))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('filtfilt')

sig = np.cumsum(randn(800))  # Brownian noise
sig_ff = signal.filtfilt(b, a, sig)
sig_lf = signal.lfilter(b, a, signal.lfilter(b, a, sig))
plt.subplot(2, 1, 2)
plt.plot(sig, color='silver', label='Original')
plt.plot(sig_ff, color='#3465a4', label='filtfilt')
plt.plot(sig_lf, color='#cc0000', label='lfilter')
plt.grid(True, which='both')
plt.legend(loc="best")


What is the low-pass filter that gives minimal transient response?

As mentioned in the comments there is no unique answer to your question, because you're asking for a contradiction. A good low pass filter (i.e. one with a fast transition between pass band and stop band, and with high stopband attenuation) will exhibit a considerable transient response. On the other hand, a filter with a very smooth step response is worse as a frequency selective filter than other filters of the same order. However, a Bessel filter might be a good compromise between these two conflicting goals. It performs reasonably well as a low pass filter, and it has very little overshoot in its step response.


Python audio analysis: which spectrogram should I use and why?

I recommend a Synchrosqueezed Continuous Wavelet Transform representation, available in ssqueezepy. Synchrosqueezing arose in context of audio processing (namely speaker identification), and there's much literature on applying CWT for audio tasks.
Advantage over STFT is the inherently logarithmic nature of the feature extractor, matching audio structuring. That is, we perceive sound differences in relative terms, so 6kHz is to 3kHz what 24kHz is to 12kHz. STFT is unable to map such features without being very wasteful, as it lacks an adaptive window (your Spectrogram A is not logarithmic; it's a logarithmic display of linear data). Analytic CWT (CWT with analytic wavelet) further enjoys superior instantaneous frequency and amplitude representations (see below references).
Synchrosqueezing further enhances these representations, and can be thought of as an attention mechanism. Experiments have shown it to enhance EEG classifier performance, for example. Currently only CWT is supported, but I'll be adding STFT within a week or so. Also see:

Excellent CWT tutorial (also comparing with STFT)
Synchrosqueezing friendly overview
CNN + CWT vs STFT comparison
Synchrosqueezing paper

Examples on EEG data:


Is there a widely available implementation of an adaptive, recursive, numerically stable IIR filter?

While I may not be the final authority, I have looked into this. In short, no.

From my own research, the problem is stability. When the filter is extended to have a wider bandwidth, which of course means that the y[n] values and lags 'move quickly', and then is quenched to a lower bandwidth, where the y[n] values and lags 'move slowly', the y[n] buffer may be loaded with fast oscillated entries that appear as initial conditions to the low bandwidth filter. The filter can easily go unstable. If the bandwidth is altered in an adiabatic manner this instability may be avoided (I'm supposing) but then there's an additional constraint on your adaptivity.

Now, you don't have to take my word for it. In "Adaptive Filters: Theory and Applications" (2nd) page 5 (Amazon), Boroujeny writes, "However, as we shall see in the later chapters, because of the many difficulties involved in adaptation of IIR filters, their application in the area of adaptive filters is rather limited." The author states that the adaptation process can place (digital) poles outside of the unit circle, even if they started within the unit circle. The rest of his book is mostly dedicated to FIR filters.

I hope this information helps.


Examples of mixture of 'almost periodic' signals/data in nature or practice

The sound produced using the lowest pitched strings of large stringed music instruments often includes slightly inharmonic overtones (due to finite string thickness and stiffness), which renders the waveform of the note pitch not quite, but "almost periodic".  You can record these sounds from the bottom octaves of most (physical, stringed) pianos, and bass guitars, etc.


Basic Understanding of Convolution Operation

Since $y[n-n_0] = y[n] \ast \delta [n - n_0]$ and convolution is both associative and distributive, is it not the case that

$$
y[n-n_0] = \left(x[n] \ast r[n] \right) \ast \delta[n-n_0] = x[n] \ast r[n-n_0] = x[n-n_0] \ast r[n]?
$$


Can I calculate a half-zero data's FFT from its non-zero front part's FFT?

The computationally efficient way is to do an IFFT/zero-pad/FFT.  It's usually less efficient to do a circular Sinc interpolation (to create the effect of a rectangular window in the other domain) for any large number of extra points.  Any other interpolation is less accurate.


Why does the Hamming window attenuate the FFT?

You can figure the change in energy applied to your signal by analyzing the energy in your window.  The Hamming window's energy (or scaling) is well known and Google will tell you what it is.  You are doing something different. Can you post your code for further help?

From your code, ans your image, it really appears that you applied the window to your transformed data.  The window (in this case) is applied to the time data before the transform.  The window will reduce sidelobes from strong frequency components.  Take note that taper windows cause the so-called "spectral leakage" where energy leaks into other frequency bins


What is the role of complex exponential?

Euler's relationship says that $e^{j\Theta}$ is equal to $cos(\Theta) + j*sin(\Theta)$.  The Fourier Transform can then be seen as correlating the signal with sinusoids at various frequencies.  The continuous Fourier Transform correlates with an infinite number of sinusoids, while the discrete transform uses $N$ sinusoids, where $N$ is the length of the transform.


Physically, what does the usage of two variables mean for convolution

My intuition of convolution is that it is just a way to depict multiplication of two signals where each signal is made up of various frequencies and phases.

Not really, no.  I can't think of a case where you'd actually receive two signals and convolve them (correlate them -- yes, but not convolve).  This may be why you're not understanding the rest.
Convolution is used to find the output of a linear system* when you know the system's impulse response.

And intuitively, convolutions are seen as sliding windows.  ... $\tau$ is very clearly some dummy variable but beyond that I feel a little confused.

$\tau$ is the dummy variable that implements the sliding window in the convolution.  For each time $t$, you multiply the signal by the time-reversed impulse response, and you integrate that product.  $\tau$ is just the variable over which you do the integration.
* I purposely left the "time-invariant" part of that out -- convolution applies to time-varying systems, too, it's just more complicated.


There is no problem in your plot if by fft you mean DTFT (or Discrete Time Fourier transform) and it is true to say "the points on this dash line are the DTFT result of corresponding frequency":
If your function is a cosine:
$$
x[n] = \cos(\omega_0n) \Longrightarrow X(e^{j\omega}) = \pi\delta(\omega-\omega_0)+\pi\delta(\omega+\omega_0)\ ; \qquad|\omega|<\pi
$$
and your window function is defined as follow:
$$
w[n] = \left\{
        \begin{array}{ll}
            1 & \quad |n| \leq N \\
            0 & \quad |n| > N
        \end{array}
    \right.
\Longrightarrow W(e^{j\omega}) = \frac{\sin\left(\omega(\frac{N}{2}+1)\right)}{\sin\left(\frac{\omega}{2}\right)}
$$
as you mentioned if we multiply two functions in time domain, it's equal to periodic convolution in frequency domain:
$$
y[n] = x[n]w[n]\\ \Longrightarrow Y(e^{j\omega})=X(e^{j\omega})\circledast W(e^{j\omega})=\frac{1}{2\pi}\int_{2\pi}^{}X(e^{j\omega})W(e^{j(\omega-\theta)})  d\theta
$$
in this case (in one period):
$$
Y(e^{j\omega})=\frac{1}{2}\frac{\sin\left((\omega-\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}+\frac{1}{2}\frac{\sin\left((\omega+\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)}
$$
if you plot this function in one period for $\omega_0=\frac{\pi}{5}$ and $N=2$ it look like this:

but if you change the window length to $N=4$ the result changes to:

and the cause for this, is that two main lobes overlap with one another.


Fourier transform diagonalizes time-invariant convolution operators

TL;DR: Fourier is the logarithm of the convolution.
Doing a convolution (linear time-invariant operator) is a quite involved processing in the original (primal) domain: it involves intricate sums and products of an operator $h$ and signals $s$ at different indices:
$$ \sum_k h[k]\times x[n-k]\,.$$
With a Fourier transform, you transform the whole information for the operator $h$ and data $x$ in the primal domain onto a novel dual (frequency) domain, in $\hat{h}$ and $\hat{x}$. Then the result of convolution becomes:
$$  \hat{h}[m] \times \hat{x}[m]$$
which is a simple index-by-index product. Hence the term "diagonalization": from products of two indices, $[k]$ and $[n-k]$, Fourier turns them into a product of a single index $[m]$ times  $[m]$.
To make the above a bit simpler, Fourier is to convolution what logarithms are to multiplication. Let us forget the carry for a moment. Doing a multiplication is complicated, especially for large numbers: one has to multiply thousands to tens and add them to products of hundreds.  Logarithms turn multiplication into adds: you (almost) have to add figures in the same position (lest the carry), instead of combining figures in different powers.
For ages, people have used logarithm tables or logarithmic sliding rulers as analog computers. Logarithmic tables and rulers have been very important in days before computers, to allow computations in engineering, astronomy, etc.
Fourier principles, tools and fast librairies are the nowadays analogs of logarithms, that helped shaped our digital world. The FFT algorithm is often considered one of the most important algorithms. Indeed, there exists FFT based multiplication of large numbers.


Can I decompose a sine signal into its component sine waves?

Curve fitting is another option after Fourier transform, you can fit a sum of sine function to your signal and the fitting coefficients show the amplitude and frequency of the signal.
Check this: https://stackoverflow.com/questions/38832474/curve-fitting-in-matlab-for-a-sinusoidal-function-with-more-than-8-terms


One reason is that higher frequencies are envisioned. With higher frequencies, the path loss grows (cf. Friis equation). Also, the wavelength is reduced and thus, $\lambda/2$ radiators start becoming quite small. The power you can radiate from a small aperture is limited as well so overall our power budget suffers. The only way out is directivity (cf. Friis equation again): we need to make sure our energy is directed towards the users of interest. This is possible by employing large antenna arrays (i.e., arrays comprised of a large number of elements). If each element radiates a certain power, $N$ elements can add up coherently to the $N$-fold power, provided we send it in the right direction. This is what beamforming does.
Another reason is that beamforming allows us to spatially separate users so that communication links can share the same time and frequency resources and be separated only spatially. In theory, an $N$-element antenna array can serve $N$ data streams at the same time and frequency, either to one user (that receives $N$ spatially multiplexed data streams, as it is common in WiFi) or to $N$ users (that each receive one). This increases network capacity even more.
Both factors play well together: higher frequencies and smaller wavelengths lead to a more rapid decorrelation of the channel in space so that even relatively closely spaced users can still be spatially multiplexed. Think of very narrow "pencil" beams, each pointed to a different receiver (if we are serving their lines of sight) or path (if we are using multipath).
While this sounds good, it also poses significant challenges. If we only slightly depoint the beam, the signal may be lost entirely. For moving users, this can be quite a challenge. Initial acquisition becomes harder too because signals transmitted/received omnidirectionally might be too faint to detect.


Extra frequencies after applying low pass filter

The issue is that you're violating Nyquist which requires that the sampling rate be at least twice the rate of the highest frequency in the signal. In your case all of the sines that you're generating above 5000 Hz are folding back. For example, a 5100 Hz tone (5000+100) would fold back to 4900 Hz (5000-100) and so on.

The 1 kHz tone you highlighted is the 9 kHz tone you generated, but folded back. If you look at your graph more closely you'll notice that this tone is in the unfiltered signal and is in no way related to the filtering.


Decision directed TED (timing error detector). What does it mean?

"Decision directed" means that the inputs to the timing error determination use decisions (best guesses) of what the ideal symbol is based on the sample read.
This is in contrast to "symbol-aided" where a known symbol is transmitted and the replica symbol can be used in the receiver. You can visualize how this would allow for much further error determination due to a much wider range before the results alias.
I have diagrams specific to phase detection, but well demonstrates the concept and utility of decision-directed approaches applicable to timing recovery as well. As an example of this for the purpose of phase detection (used in carrier recovery) below shows a Decision Directed Phase Detector for QPSK and showing at the ideal sampling location in the symbol the actual sample measured at $V_2$ from which we assume (best guess) that the symbol transmitted must have been $V_1$. This is just an assumption based on nearest neighbor and only serves to drive the correction loop to that state, even though the whole symbol may end up being rotated once that correction converges (and then other techniques such as pilot symbols or higher level inspection or the data itself can be used to determine if such rotation occurred).

In contrast below shows a QAM Constellation where a symbol-aided approach is used for phase detection. The lock range as limited by the nearest-neighbor feature of a decision directed approach is greatly extended as depicted in the figure where this shows an actual symbol received as the dark circle and based on the known location (the outer symbol at +45 degrees) a much greater error can be determined. In a Symbol-Aided approach the known symbols are transmitted over a portion of the symbol such as in a pilot which is then used for carrier and timing recovery, and can operate in much lower SNR conditions. Decision-directed approaches out-perform  other approaches at higher SNR levels where all samples are based on the received signal itself (such as the Gardner TED which typically does not use a pilot or decisions) since a decision when good removes all noise for that input to the timing detector.


These fall under signal processing and mostly you will find electrical engineers who are working in domain of speech recognition. The thing with computer vision is that it overlaps between both field.But as you mentioned you only have option of Electronics/Communications, Electronics is not much related to these field ,well in communications, we all use statistical signal processing which is also useful in speech processing and computer vision.

If you can get into Electrical Engineering(signal processing) that would be much better but out of these two I would suggest go for Computer Science if you plan a future in these field.


How does voice affect a sound wave or its corresponding waveform?

I would view it as the following:
The letter A is the input to a system, one is your vocal chords the other one mine. Since the system response of your and my vocal chords are different, we would get different outputs.
Now, what would characterize the system response, it would be factors like, rise time/bandwidth, decay response etc. All of this would impact the tone and pitch terminologies usually used in audio/speech. These variables will vary for each user. With different input the response will then also be different . So "b" will not be same as "a"


MATLAB Image-Processing Toolbox: asking for a roadmap!

If the distinct contents/ingredients really have different colors than you just have to take a picture, always in the same location/distance, and then count the amount of pixels for each color.

You can transform your RGB image to the HSV color space to make it easy to identify specific colors.

Some pseudocode:

image_hsv = rgb2hsv(image_rgb)
greenmax = 120+10
greenmin = 120-10
num_green_pixels = sum( and( image_hsv[H]>greenmin, image_hsv[H]<greenmax ) )


How to filter data to remove noise in MATLAB?

You can try to remove noise (especially occasional spikes) by non-linear filter. I suggest to use median filter (http://en.wikipedia.org/wiki/Median_filter), as your spikes has length only few points (as I see). There is matlab function - medfilt1. It is probably that such filtration will be enough.


Confused about DFT symmetric behavior for real-valued signals?

$N=1000$ and your signal frequency is $10000$. So there should be a peak at bin $k=10$. And another one at $N-k=1000-10=990$. However since MATLAB starts indexing at $1$ instead of zero, you see peaks at $11$ and $991$.


Nop. For an infinite sequence, the convergence of Z transform depends on the value of $z$ also. So as far as $n\rightarrow\pm\infty$, equation (1) will not be satisfied for all values of $z$.

Consider a simple case, $x(n) = u(n)$
$$X(z) = \sum_{n=0}^{\infty}z^{-n}$$
$$X(z) = 1 + z^{-1}+z^{-2}+z^{-3}+ \cdots$$
for $z=0.5$,
$$X(z) = 1 + 2+4+8+ \cdots = \infty$$
for this sequence to converge, $|z|>1$.


Look at the forms of the functions, and remember the duality property.

I invite you the read the Wikipedia article: http://en.wikipedia.org/wiki/Fourier_transform


How Do I Plot Color Distribution of HSV Image in MATLAB

Have a look at File Exchange Submission - RGB/HSV Distribution on  visualization of RGB/HSV distribution in a given image:

hsv_distribution(inputImage, 5)

It should do what you're after as the result is like:


Convolution by differentiation property of convolution

First of all let me state that, the property of differentiation under (convolution) integral operator may not be holding true for arbitrary signals but only for some subclass, such as absolutely integrable ones. That's stated in the wikipedia link OP provides about convolution. All I can make sense from that statement is that interchange of differentiation and integration is only possible if the integral is convergent. That's formal mathematics and I'm not very sure whether that applies to generalized functions, such as $u(t)$ and $\delta(t)$ as well. So I'll leave the discussion to mathematicians for a moment and proceed with the solution, treating the functions as if they're well behaving.

So far you've correctly obtained the signal $y(t)'$ from the differentiation of convolution property as: $$y(t)' = h(t)' \star x(t)  = -3 u(t) + 2u(t-1) + u(t-5)$$

Now it comes to proper integration of $y(t)'$ to obtain $y(t)$ which was the most misleading step; $$ \int_{-\infty}^{t} y(\tau)' d\tau = \int_{-\infty}^{t} -3 u(t) + 2u(t-1) + u(t-5) d\tau$$

the left side evaluates as;

$$ \int_{-\infty}^{t} y(\tau)' d\tau = \int_{-\infty}^{t} dy d\tau = y(\tau)+C |_{-\infty}^{t} = (y(t)+C)-(y(-\infty)+C) = y(t)-y(-\infty) $$

Now here I shall make a few simple comments. As below we argue with @msm whether this is a definite integral or an indefinite one? He seems to claim it as an idefinite one with constant of integration at the output. (I don't actually know it because he doesn't put his equation there inspite of all my attempts, also he seems(?) to delete that comment as well? Rather he says it's simple to see. I don't know. It's always better to make it clear when requested to avoid any confusion. Of course this is a definite integral with clear limits on it, hence at the output you won't have a constant of integration whatsoever. The integrand is antiderivatived, and evaluated at the limits and subtracted so that the constant of integration is cancelled. What remains there as $y(-\infty)$ the integrand evaluated at the limits. For causal systems and signals applied, $y(-\infty)$ is by defintion zero, hence it's ignored, but when you have a noncausal system (like in the example system with $h(t) = u(-t-1)$ you may not ignore it and hence the term $y(-\infty)$ remains.

The right side of the equation is integrated as follows:
$$  = \int_{-\infty}^{t} -3 u(\tau) + 2u(\tau-1) + u(\tau-5) d\tau  = \int_{0}^{t} -3 u(\tau) + 2u(\tau-1) + u(\tau-5) d\tau $$

$$  = (-3 r(t) + 2r(t-1) + r(t-5) + C)-(-3 r(0) + 2r(0-1) + r(0-5) + C)$$
$$  = -3 r(t) + 2r(t-1) + r(t-5) $$

Where we have indeed antiderivated each integrand and subtracted their evaluation at $\tau=t$ and $\tau=0$, the latter one all providing $0$

Hence we have;
$$ y(t) = -3 r(t) + 2r(t-1) + r(t-5) + y(-\infty) $$

At this point it's worth noting that, the output $y(t)$ depends on itself due to the $y(-\infty)$ term. Where you can fortunately (easily) evaluate it from the convolution integral as $$y(t) = \int_{-\infty}^{\infty} h(t-\tau)x(\tau)d\tau $$ as $t \rightarrow -\infty$ resulting in;
$$ y(t) = \int_{-\infty}^{\infty} h(-\infty-\tau)x(\tau)d\tau = y(t) = \int_{-\infty}^{\infty} 1 x(\tau)d\tau = 7 $$

Hence the solution becomes, $$y(t) = 7 -3 r(t) + 2r(t-1) + r(t-5) $$
and finally $$y(4) = 7 - 3 r(4) + 2 r(3) + r(-1) = 7 - 3 \times 4 + 2 \times 3 + 1 \times 0 = 1$$

Below you can see the verification by direct evaluation of the colvolution integral.

(EDIT: see the first paragraph to get the point)
As warned by the OP, that the differentiation property under the convolution is not valid for the signal $h(t) = u(-t-1)$, so we then resort to direct evaluation to check the answer.

$$y(t) = x(t) \star h(t) = x(t) \star u(-t-1)$$

The graphical method provides the easiest evaluation of this integral for particular output value of $t=4$ as:

$$y(4) = \int_{\tau= 5}^{\tau=6} x(\tau)h(4-\tau) d\tau = \int_{\tau= 5}^{\tau=6} 1 \cdot 1 d\tau = 1 $$

So the requested value is: $y(4) = 1$ , hence you are right the answer is $y(4) = 1$ and (EDIT: see the first paragraph to get the point) [it seems the differentiation property is not applicable for arbitrary signals but those which are absolutely integrable.]

An example to see the diferentiation property is; $x(t) = u(t) - u(t-1)$ and $h(t) = u(t) - u(t-1)$  where we have $y(t)' = u(t) - 2u(t-1) + u(t-2)$ and the differentiation property yields: $$y(t)' = h(t)' \star x(t) = (\delta(t) - \delta(t-1)) \star x(t) = u(t) - 2u(t-2) + u(t-2)$$

Sorry for the confusion...


Why negative harmonics copies in sampling

First, about "why we have negative frequencies and which one is copied"...
In a real signal, there is no difference between a negative frequency and a positive frequency in this sampling context. Consider that you have captured 1000 samples that you know contains a single sine wave. You analyze it and find it contains exactly ten cycles of a sine wave. Knowing the sample rate, you determine that it's a 100 Hz sine wave. Now analyze the capture from the other direction—end to start. You would also find that it's 100 Hz, but since we moved backwards through time, it would be -100 Hz.
That is, if you choose to plot negative frequencies, they will always be a mirror of the positive frequencies. You will never find -80 Hz without 80 Hz, for instance. There is never a choice of "which on is copied". If you find somehow that there is a -80 Hz component, you can be sure there is an 80 Hz component.
Sampling is a type of amplitude modulation (AM). Pulse Amplitude Modulation (PAM) in the analog domain, Pulse Code Modulation (PCM) in the digital domain. In the analog domain, it's a multiplication of the signal with an impulse train at the sampling frequency. Looking at just the first harmonic of the impulse train, it's a sinusoid (cosine in particular) at the sampling frequency. In your case this means you are multiplying your signal by 800 MHz (as well as its pulse harmonics, but we need only consider the closest one here).
Amplitude modulation produces sum and difference frequencies, but for this problem we only care about the difference frequencies at the moment. That means if your signal has a 120 MHz first harmonic, sampling will also produce 800-120 HMz = 780 Mhz. And the fourth harmonic of your signal, 480 MHzHz, would produce a new component of 800-480 = 320 MHz. Note that this is exactly the same answer as your "-4*120M+800M". Except that I came up with my answer without considering negative frequencies (800M minus 4x120M).
To summarize, negative frequencies always have positive counterparts. And there is no particular need to consider negative frequencies with what you are looking at—it seems it just adds confusion, but may make more sense in context of the document you refer to.


The difference about WPCN and TDMA or FDMA

TDMA and FDMA are paradigms; they provide different methods for sharing communications resources between multiple users.  WPCN is a technology for using wireless communications resources, and can be used with TDMA or FDMA depending on the particular application, though neither is "best suited" to the technology.


How can I create a histogram showing the RMS amplitude for each frequency bin of an audio file

If what you want is to determine the amount of energy your signal has in particular frequency bands you can perform a FFT of your signal. For that you need to make sure that you have amplitude data and not energy data. So you would want to have sound pressure Pa and not Pa^2. After having calculated a FFT with that the values have to be multiplied by $1/\sqrt{2}$ to get to the RMS-format.
Depending on your time data length and sampling time you get a specific frequency resolution of your signal in the frequency domain. With $N$ time samples and a sampling time of $\Delta t$ you get a frequency resolution of $1/(\Delta t \cdot N)$ Herz. Now if you want to know the energy contribution of a specific frequency range you can calculate the RMS-value of that frequency range using the earlier calculated linear spectra in RMS-format. According to Parseval's theorem one can use the following equation to calculate the RMS-value of a specific frequency range:

Therein the $A_k$s are the amplitudes of different frequencies in the spectra:

I used that approach to calculate the Energy of 10 Hz bins of a signal containing two sines with amplitude $2$ and Frequencies $5$ and $45$ Hz. Sampling frequency is 0.0098.
Here is the time Signal:

The resulting linear Spectra with Peak to Peak values is:

From that I calculated the RMS-Values of $\Delta 10$ Hz bins. So the first bar indicates the RMS-Value of the signal from $0 - 10$ Hz. The second from $10 - 20$ Hz and so on.

If some noise is added the spectra changes to:

The RMS-values of the $\Delta 10$ Hz bins become then:

I further suggest looking at Root Mean Square and Overall Level
Here is the Matlab Code I used to do the above (probably still containing errors though):
clear all
close all

t=linspace(0,10,1024); % time vector
y=2*sin(5*2*pi*t)+2*sin(45*2*pi*t)+randn(1,length(t)); % time signal
Fs=1/(t(2)-t(1));
Fn=Fs/2;
N=length(y);
yf=fft(y,N);
df = Fs/N; % Frequency Resolution

plot(t,y)
grid on
title('Time Signal')
ylabel('Amplitude')
xlabel(['Time Resolution: ',num2str(1/Fs),' s'])

amplH = abs(yf);
amplitudengang = fftshift(amplH/N);

x_fn = 0 : df : Fn-df;
x_fa = 0 : df : Fs-df;

figure
stem(x_fa-Fn, amplitudengang, 'b.-')
axis([-Fn Fn 0 max(amplitudengang)])
title('Two sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),' Hz'])
grid

amplitudengang=[amplitudengang(513) amplitudengang(514:end).*2];

figure
stem([0:df:(Fn-df)], amplitudengang, 'b.-')
axis([0 Fn 0 max(amplitudengang)])
title('Single Sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),'Hz'])
grid

amp10=amplitudengang(1:round((10-df)/df));
amp20=amplitudengang(round((10-df)/df):round((20-df)/df));
amp30=amplitudengang(round((20-df)/df):round((30-df)/df));
amp40=amplitudengang(round((30-df)/df):round((40-df)/df));
amp50=amplitudengang(round((40-df)/df):round((50-df)/df));

rms10=amp10./sqrt(2);
rms20=amp20./sqrt(2);
rms30=amp30./sqrt(2);
rms40=amp40./sqrt(2);
rms50=amp50./sqrt(2);

rms10=sqrt((rms10(1)/2)^2+sum(rms10(2:(end-1)).^2)+(rms10(end)/2)^2);
rms20=sqrt((rms20(1)/2)^2+sum(rms20(2:(end-1)).^2)+(rms20(end)/2)^2);
rms30=sqrt((rms30(1)/2)^2+sum(rms30(2:(end-1)).^2)+(rms30(end)/2)^2);
rms40=sqrt((rms40(1)/2)^2+sum(rms40(2:(end-1)).^2)+(rms40(end)/2)^2);
rms50=sqrt((rms50(1)/2)^2+sum(rms50(2:(end-1)).^2)+(rms50(end)/2)^2);

figure
stem([5 15 25 35 45],[rms10,rms20,rms30,rms40,rms50])
axis([0 Fn 0 max([rms10,rms20,rms30,rms40,rms50])])
title('RMS values per 10Hz')
ylabel('Amplitude')
xlabel(['Frequency in Hz'])
grid on

Edit
To reduce leakage the time data needs to be windowed. After that the resulting spectra needs to be multiplied by the matching correction factor:


Effects of linear interpolation of a time series on its frequency spectrum

Duane Wise and i wrote a paper back in the 90s that we presented to an AES convention that spelled out how to model time-domain polynomial interpolation (of which linear interpolation is an example) in the frequency domain.

i think you can get a copy here: Performance of Low-Order Polynomial Interpolators in the Presence of Oversampled Input


Continuous-time frequency response question

Is it possible to express $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ in terms of $X_c\left(j\Omega \right)$ using some continuous-time filter?

an LTI system cannot do that.  continuous-time "filters" need not be LTI but usually we mean LTI when we consider "filters" in analog or digital filters.

you can change the spectrum of $X_c\left(j\Omega \right)$ to $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ by speeding up (or slowing down) the time index of $x_c(t)$ to $x_c\left(t \frac{T_1}{T_2}\right)$.

it appears you want a pitch shifter or similar.  is that it?


The states depend on the form of an IIR filter. Take second order filters for simplicity. The following picture shows the structure of a transposed direct form II filter.

The differential equation is
$y(n) = b_0  x(n) + s_1(n)$
and the filter states can be calculated by recurrence
$s_1(n) = b_1 x(n-1) + s_2(n-1) - a_1 y(n-1)$
$s_2(n) = b_2 x(n-1) - a_2 y(n-1)$
As for (transposed) direct form I, it should have a filter states with 4 elements. The following is the transposed direct form I implementation of a second order filter.


Should we use all $N$ correspondence points to compute Homography or just $4$ points?

It's generally better to use all the points you have, since often, errors are independent per point.

given my correspondence points are good, shouldn't there be an exact
solution

Sure, but in realistic situations, it's not usually the case and/or you don't know that the points are good / which points are good.
A typical approach is to use all the points you have, but reject "bad" points using outlier-rejectiong  (e.g. RANSAC ) These slides have an example of such a procedure.


Applying Circular Cross Correlation in MATLAB

Cross correlation can be implemented in the frequency domain using FFT by multiplying signals.
For a cyclic cross correlation you can use:
cc = ifft(fft(x).*conj(fft(y)));
You can also calculate the linear cross correlation using FFT by zero-padding the signals before the FFT.
Related: Cross correlation with FFT and fftshift
Edit:
In order to calculate the cross correlation of two images you have 2 options:

x and y should be the grayscale representation of the images
Repeat the process 3 times. Each time x and y are a single channel (r,g or b) from each image


You need to be careful with the definition of the DC value of a signal. The actual time average, which is often called DC value is given by

$$\overline{s(t)}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}s(t)dt\tag{1}$$

whereas the value of the signal's Fourier transform at $\omega=0$ is given by

$$S(0)=\int_{-\infty}^{\infty}s(t)dt\tag{2}$$

If $(2)$ is finite, $(1)$ equals zero, and if $(1)$ is finite but non-zero, $(2)$ doesn't exist, and the Fourier transform $S(\omega)$ has a Dirac impulse at $\omega=0$. See also this related answer.

For the given signal $s(t)=1/t$ both integrals in $(1)$ and $(2)$ do not converge in the conventional sense. However, their Cauchy principal value exists and equals zero, as mentioned in a comment by Robert Bristow-Johnson. Hence, the DC-value of $s(t)=1/t$ equals zero, regardless whether you define it by $(1)$ or by $(2)$.

Note that the Fourier transform of $s(t)=1/t$ equals

$$S(\omega)=-j\pi\,\textrm{sgn}(\omega)\tag{3}$$

which is just the frequency response of a Hilbert transformer (scaled by $\pi$).

From $(3)$, $S(0)$ doesn't exist, but the average of the left-sided and right-sided limits $\frac12 (S(0^-)+S(0^+))$ equals zero, and - using the Cauchy principal value according to the definition of the Hilbert transform - the Hilbert transform of a constant is indeed zero.


Help in sketching a CT signal with discontinuities

You are right, the solution manual is wrong and here is a piecewise linear definition of the signal $[x(t)+x(-t)]u(t)$ based on the acceptance of $\text{u}(0) = 0.5$ (i.e. the value at a discontinuity is the half way between right and left ends)
$$[x(t)+x(-t)]u(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t < 0}}\\&[1.5] &\scriptstyle{\text{, t = 0}}\\&3 &\scriptstyle{\text{, 0<t<1}}\\&[2] &\scriptstyle{\text{, t = 1}}\\-&2t+3 &\scriptstyle{\text{,  1<t<2}}\\&[-0.5] &\scriptstyle{\text{, t=2}}\\&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Where the values within $[.]$ are points of discontinuities.

Let me put the steps to reach this conclusion in more details:
Looking at the provided graph of $x(t)$ we can deduce a piecewise definition of the signal. However we must be causious about the possible discontinuities, for the sake of mathematical (continuous) function theory, those isolated points of discontinuities can be omitted entirely or as an engineering convention the value of the signal can be defined to be the half of the sum of its values at the right and left limits of the function at the discontinuity point. I'll accept this second definition, as DSP is concerned more with engineering than with pure mathematics.

Therefore, we have the definition of $x(t)$ as follows:
$$x(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t < -2}}\\
&[-0.5] &\scriptstyle{\text{, t = -2}}\\
&t+1 &\scriptstyle{\text{, -2<t<-1}}\\
&[0.5] &\scriptstyle{\text{, t = -1}}\\
&1 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&2 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5] &\scriptstyle{\text{,  t=1}}\\
&-t+2 &\scriptstyle{\text{,  1<t<2}}\\
&0 &\scriptstyle{\text{,  t>2}}\\
\end{align}
\end{cases}$$

From which we can also deduce $x(-t)$ as follows:
$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  -t < -2}}\\
&[-0.5] &\scriptstyle{\text{, -t = -2}}\\
&-t+1 &\scriptstyle{\text{, -2<-t<-1}}\\
&[0.5] &\scriptstyle{\text{, -t = -1}}\\
&1 &\scriptstyle{\text{,  -1<-t<0}}\\
&[1.5] &\scriptstyle{\text{, -t=0}}\\
&2 &\scriptstyle{\text{,  0<-t<1}}\\
&[1.5] &\scriptstyle{\text{,  -t=1}}\\
&t+2 &\scriptstyle{\text{,  1<-t<2}}\\
&0 &\scriptstyle{\text{,  -t>2}}\\
\end{align}
\end{cases}$$

Observe carefully how every occurence of $t$ is replaced with $-t$ both in the defining formula of the function and also in the definition of its intervals. We shall reformat intervals of $x(-t)$ with positive $t$ as follows:

$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t > 2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&-t+1 &\scriptstyle{\text{, 1<t<2}}\\
&[0.5] &\scriptstyle{\text{, t = 1}}\\
&1 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&2 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{,  t=-1}}\\
&t+2 &\scriptstyle{\text{,  -2<t<-1}}\\
&0 &\scriptstyle{\text{,  t<-2}}\\
\end{align}
\end{cases}$$

Now adjust the order of intervals of $x(-t)$ so that they are in left to right order:

$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t<-2}}\\
&t+2 &\scriptstyle{\text{,  -2<t<-1}}\\
&[1.5] &\scriptstyle{\text{,  t=-1}}\\
&2 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&1 &\scriptstyle{\text{,  0<t<1}}\\
&[0.5] &\scriptstyle{\text{, t = 1}}\\
&-t+1 &\scriptstyle{\text{, 1<t<2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Now we can add these two signals for each respective interval. In fact we are here a bit lucky that the signals do have matching intervals, which would not be the case if the original signal $x(t)$ involved nonsymetric pieces in the first place. The addition produces the following signal:

$$x(t)+x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t<-2}}\\
&[-0.5] &\scriptstyle{\text{,  t=-2}}\\
&(t+1)+(t+2)=2t+3 &\scriptstyle{\text{,  -2<t<-1}}\\
&[0.5]+[1.5]=[2] &\scriptstyle{\text{,  t=-1}}\\
&1+2=3 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5]+[1.5]=[3] &\scriptstyle{\text{, t=0}}\\
&2+1=3 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5]+[0.5]=[2] &\scriptstyle{\text{, t = 1}}\\
&(-t+1)+(-t+1)=-2t+3 &\scriptstyle{\text{, 1<t<2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Finally we shall multiply this with the funciton $\text{u}(t)$ whose piecewise definition is:
$$\text{u}(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{, t<0}}\\
&[0.5] &\scriptstyle{\text{, t=0}}\\
&1 &\scriptstyle{\text{, t>0}}\\
\end{align}
\end{cases}$$

And the resulting function $[x(t)+x(-t)]\text{u}(t)$ was already provided in the first part above...


Digital implementation of first order analog filter using bilinear transformation

The expected 3dB frequency is wrong because of radian conversion. With $s=j\omega$,
$$H(j\omega) = \frac{1}{1+\tau j\omega} $$
and $20\log_{10}|H(j\omega)|\approx-3$
when $\omega=\omega_{3dB}=\frac{1}{\tau}$. However $\omega=2\pi f$, so
$$f_{3dB}=\frac{\omega_{3dB}}{2\pi}=\frac{1}{2\pi\tau}.$$
With $\tau=0.1\text{s}$, $\ f_{3dB}=1.59\text{Hz}$ as on the original plots.


Comparing FFT for increasing time intervals in a signal

I'm afraid the premise of your approach is faulty as can be shown with a single counterexample.  Suppose your signal is a simple real pure tone sinusoidal.  When the interval length reaches 3 cycles, bin index 3 (zero based) of the DFT will be non-zero and all the rest of the DFT bins will be zero.  (Considering only the lower half of the DFT values because a real valued signal as the upper half is a conjugate mirror of the lower half.)  Then as your interval increases the DFT starts to smear with leakage.  When you reach 3 and 1/2 cycles, the smearing will be the greatest and the magnitudes of bins 3 and 4 will be nearly equal and their phases nearly 180 degrees apart.  When you get up to 4 cycles in your sampling interval, once again all the bin values are zero except for bin 4.  That behavior will repeat for each whole number of cycles your interval reaches and therefore there will never be any kind of convergence happening.

It is not clear to me from your question, but are you assuming your signals are mixes of steady pure tones?


The arma approach is to model the current output of the system as the sum of past outputs and past inputs explicitly. The assumption of gaussian model for the noise statistics still can be used for the unpredictable signal which cannot be modeled as arma.

From system's frequency response spectrum modeling aspect: an ar model is able to model only the spectrum peaks as it fits the peaks better (using poles) and the ma model is able to model the valleys better (using zeros) in the spectrum. Hence, by using an arma model we will be better enable both peaks and valleys using a reasonable order for both ar and ma.

However, it is always possible to model a zero (in z-domain) using a infinite number of poles so an approximation to an arma model (whose parameter estimation is relatively involved) is to use a high order ar model (whose paramter estimation is relatively easy).


Why are magnitudes normalised during synthesis (IDFT), not analysis (DFT)?

Whether you scale the output of your DFT, forward or inverse, has nothing to do with convention or what is mathematically convenient. It has everything to do with the input to the DFT. Allow me to show some examples where scaling is either required or not required for both the forward and inverse transform.

Must scale a forward transform by 1/N.

To start with, it ought to be clear that to analyze a simple sine wave, the length of the transform should be irrelevant, mathematically speaking. Suppose N=1024, Freq=100 and your signal is:

f(n) = cos(Freq * 2*Pi * n/N)

If you take a 1024 point DFT of f(n), you will find that bin[100] = 512. But this isn't a meaningful value until you scale it by N. 512/1024 = 1/2 and of course, the other 1/2 is in the negative spectrum in bin[924].

If you double the length of the DFT, N=2048, the output values would be twice those of the 1024 point DFT, which again, makes the results meaningless unless we scale by 1/N. The length of the DFT should not be a factor in these sorts of analysis. So in this example, you must scale scale the DFT by 1/N.

Must not scale a forward transform.

Now suppose you have the impulse response of a 32 tap FIR filter and want to know its frequency response. For convenience, we will assume a low pass filter with a gain of 1. We know that for this filter, the DC component of the DFT must be 1. And it should be clear that this will be the case no matter the size of the DFT because the DC component is simply the sum of the input values (i.e. the sum of the FIR coefficients).

Thus, for this input, the DFT is not scaled by 1/N to get a meaningful answer. This is why you can zero pad an impulse response as much as you want without affecting the outcome of the transform.

What is the fundamental difference between these two examples?

The answer is simple. In the first case, we supplied energy for every input sample. In other words, the sine wave was present for all 1024 samples, so we needed to scale the DFT's output by 1/1024.

In the second example, by definition, we only supplied energy for 1 sample (the impulse at n=0). It took 32 samples for the impulse to work its way through the 32 tap filter, but this delay is irrelevant. Since we supplied energy for 1 sample, we scale the DFT's output by 1. If an impulse were defined with 2 units of energy instead of 1, we would scale the output by 1/2.

Must not scale an inverse transform.

Now let's consider an inverse DFT. As with the forward DFT, we must consider the number of samples we are supplying energy to. Of course, we have to be a bit more careful here because we must fill both the positive and negative frequency bins appropriately. However, if we place an impulse (i.e. a 1) in two appropriate bins, then the resulting output of the inverse DFT will be a cosine wave with an amplitude of 2 no matter how many points we use in the inverse DFT.

Thus, as with the forward DFT, we don't scale the inverse DFT's output if the input is an impulse.

Must scale an inverse transform.

Now consider the case where you know the frequency response of a low pass filter and want to do an inverse DFT to get its impulse response. In this case, since we are supplying energy at all points, we must scale the DFT's output by 1/N to get a meaningful answer. This isn't quite as obvious because the input values will be complex, but if you work through an example, you will see that this is true. If you don't scale by 1/N you will have peak impulse response values on the order of N which can't be the case if the gain is 1.

The four situations I've just detailed are end point examples where it is clear how to scale the DFT's output. However, there is a lot of gray area between the end points. So let’s consider another simple example.

Suppose we have the following signal, with N=1024, Freq=100:

f(n) = 6 * cos(1*Freq * 2*Pi * n/N)  n = 0 - 127
f(n) = 1 * cos(2*Freq * 2*Pi * n/N)  n = 128 - 895
f(n) = 6 * cos(4*Freq * 2*Pi * n/N)  n = 896 - 1023

Notice the amplitude, frequency, and duration differences for the three components. Unfortunately, the DFT of this signal will show all three components at the same power level, even though the 2nd component has 1/36 the power level of the other two.

The fact that all three components are supplying the same amount of energy is obvious, which explains the DFT results, but there is an important point to be made here.

If we know the duration for the various frequency components, then we can scale the various frequency bins accordingly. In this case, we would do this to accurately scale the DFT's output:
bin[100] /= 128; bin[200] /= 768; bin[400] /= 128;

Which brings me to my final point; in general, we have no idea how long a particular frequency component is present at the input to our DFT, so we can't possibly do this sort of scaling. In general however, we do supply energy for every sample point, which is why we should scale the forward DFT by 1/N when analyzing a signal.

To complicate matters, we would almost certainly apply a window to this signal to improve the DFT's spectral resolution. Since the first and third frequency components are at the beginning and end of the signal, they get attenuated by 27 dB while the center component gets attenuated by only 4 dB (Hanning window).

To be clear, the output of the DFT can be a pretty poor representation of the input, scaled or not.

In the case of the inverse DFT, which is usually a pure mathematics problem, as opposed to the analysis of an unknown signal, the input to the DFT is clearly defined, so you know how to scale the output.

When analyzing a signal with a spectrum analyzer, analog or FFT, the problems are similar. You don't know the power of the signal displayed unless you also know its duty cycle. But even then, the windowing, span, sweep rates, filtering, the detector type, and other factors all work to goof the result.

Ultimately, you have to be quite careful when moving between the time and frequency domains. The question you asked regarding scaling is important, so I hope I have made it clear that you must understand the input to the DFT to know how to scale the output. If the input isn't clearly defined, the DFT's output has to be regarded with a great deal of skepticism, whether you scale it or not.


It's often said that pulse compression gives you a gain proportional to the time-bandwidth product (otherwise known as the pulse compression ratio, or $PCR$). This is a really misleading statement, and it had me confused enough to sit down and think about it for awhile. I thought I'd share some of my findings that I pieced together from both reading the internet and from my own DSP simulations.

The claim

Before looking into this, I had assumed that, if I take two pulses of equal amplitude $A_0 = \sqrt{P_0}$ and equal pulsewidth $\tau$, one of them a simple sine (no frequency modulation) and the other an LFM chirp, then the matched filter output $SNR$ of the chirp would be larger than that of the sine, precisely by a factor of $PCR$. I had reasoned that they both have equal energy $(E = P_0 \tau)$, but the chirp has a much narrower peak width than the sine $(\delta t = \tau/PCR)$, which would result in a larger peak height, inversely proportional to keep the total energy fixed: $$E = (P_0) (\tau) = (P_0 PCR)(\frac{\tau}{PCR}) $$

This is wrong, yet for some reason it's stated exactly this way on Wikipedia's article on the subject. Someone should fix this?

What actually happens

It turns out that both signals have the same peak power, despite the fact that the chirp's energy is squeezed into a smaller width. Therefore, we conclude that the total energy of the matched filter output of an LFM is actually smaller!

At first this looks like a problem. Moreover, you're not going to do any better with any other filter, because matched filters are optimal. But as has been noted already, the same peak power will be observed regardless of the waveform. Improved $SNR$ isn't the advantage of an LFM. The real benefit is the resolution. With the same $\tau$, I can resolve echoes that are less than $\tau$ apart, without losing $SNR$. That's pretty awesome, it turns out.

How to actually observe "Pulse Compression Gain"

If you really want to see pulse compression gain, then you've gotta compare something different. Indeed, the subject of pulse compression is usually introduced by explaining the tension between time resolution and $SNR$, before pulse compression. With a simple sine, you need to have a longer $\tau$ to improve $SNR$ (more total energy), but then you have the problem that you can't resolve echoes as well as you did before. If you want to resolve echoes better, you have to either crank up your peak power or sacrifice $SNR$.

As I mentioned above, pulse compression allows you to resolve shorter time-scales while keeping $SNR$ fixed. Alternatively, you could try comparing two signals with equal bandwidth: a simple sine wave with pulsewidth $\tau$ and an LFM chirp with pulsewidth $\alpha \tau$. Both pulses have the same resolving capability when $\alpha = PCR$ of the chirp. But the chirp will a larger peak power, and thus larger $SNR$, by a factor of precisely the $PCR$. This is simply because you have a longer pulse, not because more energy has been "squeezed" into a smaller time bin.


sparse representation for image denoising

Here are two interpretations of why sparse representations are useful for denoising:

1) Your clean image is well represented by your dictionary (usually a wavelet or 2D-DCT dictionary), but the noise isn't. This means that when calculating your sparse representation from the noisy signal, you will only capture the clean signal, and not the noise.

2) Enforcing the representation vector to be sparse (i.e. with only a few non-zero elements), means that you will eventually threshold out all the other low-energy noisy components.

As you said sparse representations can be pretty useful for compression. Compression is merely about keeping the relevant information, while discarding the irrelevant one. Well denoising is about keeping the relevant part of your signal (the one that fits your model, or your dictionary), while discarding the irrelevant one (the noise).


sound source localization - finding source vector

In general, 2 microphones are not enough.  In 3 dimensions, you actually need 4 (non co-linear) for the 3 polar coordinates.

There are cases where if you can bring more information to your problem, like knowing the signal of the source and the transmit and receive times of the signal, you can estimate the range for simple known geometries.

A good reference is:

Strang, Gilbert, and Kai Borre. Linear algebra, geodesy, and GPS. Siam, 1997


What is the difference between the Point Spread Function and Sampling Aperture?

Sounds correct to me. The point-spread function indicates how the intensity of a point source is spread among multiple pixels in an image. In an imaging sensor, the sampling aperture of a particular pixel indicates how it responds to intensity from points in space as a function of the spatial offset from the pixel's center.


Create low-frequency triangular or square waveform. Start from constant value, written into the wave. Generate one in audacity and save. Compare. Do not generate sine or anithing more compled right away. Start building complexity from simple components. Be smart, use simplicity and evolutionary approach.


What methods can be used to identify and remove echo from an audio system?

You are correct. Many methods of echo cancellation exist, but none of them are exactly trivial. The most generic and popular method is echo cancellation via an adaptive filter. In one sentence, the adaptive filter's job is to alter the signal that it's playing back by minimizing the amount of information coming from the input.

Adaptive filters

An adaptive (digital) filter is a filter that changes its coefficients and eventually converges to some optimal configuration. The mechanism for this adaptation works by comparing the output of the filter to some desired output. Below is a diagram of a generic adaptive filter:

As you can see from the diagram, the signal $x[n]$ is filtered by (convolved with) $\vec{w}_n$ to produce output signal $\hat{d}[n]$. We then subtract $\hat{d}[n]$ from the desired signal $d[n]$ to produce the error signal $e[n]$. Note that $\vec{w}_n$ is a vector of coefficients, not a number (hence we don't write $w[n]$). Because it changes every iteration (every sample), we subscript the current collection of these coefficients with $n$. Once $e[n]$ is obtained we use it to update $\vec{w}_n$ by an update algorithm of choice (more on that later). If input and output satisfy a linear relationship that does not change over time and given a well-designed update algorithm, $\vec{w}_n$ will eventually converge to the optimal filter and $\hat{d}[n]$ will be closely following $d[n]$.

Echo cancellation

The problem of echo cancellation can be presented in terms of a adaptive filter problem where we're trying to produce some known ideal output given an input by finding the optimal filter satisfying the input-output relationship. In particular, when you grab your headset and say "hello", it's received on the other end of the network, altered by acoustic response of a room (if it's being played back out loud), and fed back into the network to go back to you as an echo. However, because the system knows what the initial "hello" sounded like and now it knows what the reverberated and delayed "hello" sounds like, we can try and guess what that room response is using an adaptive filter. Then we can use that estimate, convolve all incoming signals with that impulse response (which would give us the estimate of the echo signal) and subtract it from what goes into the microphone of the person you called. The diagram below shows an adaptive echo canceller.

In this diagram, your “hello” signal is $x[n]$. After being played out of a loudspeaker, bouncing off the walls and getting picked up by the device’s microphone it becomes an echoed signal $d[n]$. The adaptive filter $\vec{w}_n$ takes in $x[n]$ and produces output $y[n]$ which after convergence should be ideally tracking echoed signal $d[n]$. Therefore $e[n]=d[n]-y[n]$ should eventually go to zero, given that nobody is talking on the other end of the line, which is usually the case when you’ve just picked up the headset and said “hello”. This is not always true, and some non-ideal case consideration will be discussed later.

Mathematically, the NLMS (normalized least mean square) adaptive filter is implemented as follows. We update $\vec{w}_n$ every step using the error signal of the previous step. Namely, let

$$\vec{x}_n = \left ( x[n], x[n-1], \ldots , x[n-N+1] \right)^T$$

where $N$ is the number of taps (samples) in $\vec{w}_n$. Notice what samples of $x$ are in reverse order. And let

$$\vec{w}_n = \left ( w[0], w[1], \ldots , x[N-1] \right )^T$$

Then we calculate $y[n]$ via (by convolution) finding the inner product (dot product if both signals are real) of $= \vec{x}_n$ and $= \vec{w}_n$:

$$y[n] = \vec{x}_n^T \vec{w}_n = \vec{x}_n \cdot \vec{w}_n $$

Now that we can calculate the error, we’re using a normalized gradient descent method for minimizing it. We get the following update rule for $\vec{w}$:

$$\vec{w}_{n+1} = \vec{w}_n + \mu \vec{x}_n \frac{e[n]}{ \vec{x}_n^T \vec{x}_n}= \vec{w}_n + \mu \vec{x}_n \frac{\vec{x}_n^T \vec{w}_n - d[n]}{ \vec{x}_n^T \vec{x}_n}$$

where $\mu$ is the adaptation step size such that $0 \leq \mu \leq 2$.

Real life applications and challenges

Several things can present difficulty with this method of echo cancellation. First of all, like mentioned before, it is not always true that the other person is silent whilst they receive your “hello” signal. It can be shown (but is beyond the scope of this reply) that in some cases it can still be useful to estimate the impulse response while there is a significant amount of input present on the other end of the line because input signal and echo are assumed to be statistically independent; therefore, minimizing the error will still be a valid procedure. In general, a more sophisticated system is needed to detect good time intervals for echo estimation.

On the other hand, think of what happens when you’re trying to estimate echo when the received signal is approximately silence (noise, actually). In absence of a meaningful input signal, the adaptive algorithm will diverge and quickly start producing meaningless results, culminating eventually in a random echo patter. This means that we also need to take into consideration speech detection. Modern echo cancellers look more like the figure below, but above description is the jist of it.

There are plenty of literature on both adaptive filters and echo cancellation out there, as well as some open source libraries you can tap into.


Why sampling rate is not an option with normalized frequency in Matlab Filter Designer

...the normalized frequency mapped $[0,\pi]$ to $[0,1]$...

The normalized frequency $F$ is defined as $f/f_s$ where $f$ is the "real" frequency (i.e. not the normalized one) and $f_s$ is the sampling frequency. The map is thus from $[0, f_s]$ to $[0, 1]$ (or to $[0, 2\pi]$ in normalized pulsation $\Omega = 2\pi F$). When using normalized frequency $F$, you thus specify frequencies as a percentage of the sampling frequency. Hence, the sampling frequency itself is not directly needed.


Difference between discrete time fourier transform and discrete fourier transform

alright, i'm gonna answer this with an argument that "opponents" to my rigid nazi-like position regarding the DFT have.
first of all, my rigid, nazi-like position:  the DFT and Discrete Fourier Series is one-and-the-same.  the DFT maps one infinite and periodic sequence, $x[n]$ with period $N$ in the "time" domain to another infinite and periodic sequence, $X[k]$, again with period $N$, in the "frequency" domain.  and the iDFT maps it back.  and they're "bijective" or "invertible" or "one-to-one".
DFT:
$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi nk/N} $$
iDFT:
$$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] e^{j 2 \pi nk/N} $$
that is most fundamentally what the DFT is.  it is inherently a periodic or circular thing.
but the periodicity deniers like to say this about the DFT.  it is true, it just doesn't change any of the above.
so, suppose you had a finite-length sequence $x[n]$ of length $N$ and, instead of periodically extending it (which is what the DFT inherently does), you append this finite-length sequence with zeros infinitely on both left and right.  so
$$ \hat{x}[n] \triangleq \begin{cases}
x[n] \qquad & \text{for } 0 \le n \le N-1 \\
\\
0 & \text{otherwise}
\end{cases} $$
now, this non-repeating infinite sequence does have a DTFT:
DTFT:
$$ \hat{X}\left(e^{j\omega}\right) = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} $$
$\hat{X}\left(e^{j\omega}\right)$ is the Z-transform of $\hat{x}[n]$ evaluated on the unit circle $z=e^{j\omega}$ for infinitely many real values of $\omega$.
now, if you were to sample that DTFT $\hat{X}\left(e^{j\omega}\right)$ at $N$ equally spaced points on the unit circle, with one point at $z=e^{j\omega}=1$, you would get
$$ \begin{align}
\hat{X}\left(e^{j\omega}\right)\Bigg|_{\omega = 2 \pi\frac{k}{N}} & = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} \Bigg|_{\omega = 2 \pi\frac{k}{N}} \\
& = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi k n/N}  \\
& = X[k]  \\
 \end{align} $$
that is precisely how the DFT and DTFT are related.  sampling the DTFT at uniform intervals in the "frequency" domain causes, in the "time" domain, the original sequence $\hat{x}[n]$ to be repeated and shifted by all multiples of $N$ and overlap-added.  that's what uniform sampling in one domain causes in the other domain.  but, since $\hat{x}[n]$ is hypothesized to be $0$ outside of the interval $0 \le n \le N-1$, that overlap-adding does nothing.  it just periodically extends the non-zero part of $\hat{x}[n]$, our original finite-length sequence, $x[n]$.