Two-View 3D reconstruction using the sparse Levenberg–Marquardt algorithm

First, please consult this tutorial about the LM minimization and this document on the bundle adjustment.

Along with many other derivative based descent algorithms, the Levenberg-Marquardt algorithm relies on the partial derivative matrix, a.k.a. Jacobian Matrix, which is the matrix of all first-order partial derivatives of a vector-valued function:

$$
\begin{align}
J_i &= \frac{\partial f(x_i, \beta)}{\beta} \\
\mathbf{J} &= \begin{bmatrix}
    \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n}\\
    \vdots & \ddots & \vdots\\
    \dfrac{\partial f_m}{\partial x_1} & \cdots & \dfrac{\partial f_m}{\partial x_n} \end{bmatrix}
\end{align}
$$

Now how you compute Jacobian depends on what you actually want from the minimization. You could for example do a full bundle adjustment where Jacobian might have an involved expression. But, let me give an example of the simplest scheme I could think of.
Triangulation from multiple views is generally written in the form of minimization of the reprojection error. Here is the cost function for the two view problem:

$$
C(X) = d(x,\hat{x})^2 + d(x',\hat{x}')^2
$$
where $x$ and $x'$ are points in 2D views and $\hat{x}$ and $\hat{x}'$ are the projections of 3D points. In the matrix notation, the error function reads:

$$
\begin{pmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_m \\
\end{pmatrix} =
\begin{pmatrix}
\|x_1^a-P_a X^1\|+ \|x_1^b-P_b X^1\|\\
\vdots \\
\|x_m^a-P_a X^m\|+\|x_m^b-P_b X^m\| \\
\end{pmatrix}
$$
where $P_a, P_b$ are the projection matrices of cameras $a$ and $b$. Since we are searching for the 3D points only $\{\mathbf{X}_i\}$, the Jacobian is $Mx3$:

$$
\begin{align}
\mathbf{J} &= \begin{bmatrix}
    \dfrac{\partial \epsilon_1}{\partial x} & \dfrac{\partial \epsilon_1}{\partial y} & \dfrac{\partial \epsilon_1}{\partial z}\\
    \vdots & \ddots & \vdots\\
    \dfrac{\partial \epsilon_m}{\partial x} & \dfrac{\partial \epsilon_m}{\partial y} & \dfrac{\partial \epsilon_m}{\partial z} \end{bmatrix}
\end{align}
$$

Of course this problem requires a good initial solution, which could be obtained from a linear closed form solution. The quality of the non-linear minimizer is highly dependent on the initialization and the quality of 2D points.

For the rest, and different parameterizations, check out the links.

Here are some remarks:

If you are not looking for re-inventing the wheel, then I would recommend Ceres Solver, which includes Automatic Differentiation. This saves the headache of Jacobians (and is much more accurate than numerical derivatives - numerical approximations are highly discouraged). Ceres also includes samples for bundle adjustment.
If you are also optimizing for the pose (camera orientation), then I do not recommend operating on the 12 parameter vector of rotation and translations. Rather choose a representation such as quaternions or the exponential map (twist coordinates). As this preserves the manifold-ness, you will get much higher quality solutions, without ambiguities. Some examples are given in the links I provided.


Hardware implementation of Gaussian noise generator

My friend, you dont need to break your head against a paper...

Just with any example of this list, the Blum Blum Shub Method:

$$x'_{{n+1}}=x_{n}^{'2}{\bmod M}$$

with $M=pq$, $p$, $q$ primes, with $x=\frac{1}{M}x' \approx U(0,1)$ gives you a fairly acceptable uniform $U(0,1)$ random number generator. You can use $p$=13331 and $q$=131 or whatever other.

And for converting the $U(0,1)$ Uniform Random Number into a $N(0,1)$ Gaussian Number Generator, as in here, just can use the Central Limit Theorem to generate N numbers each step, and make the sum (Bates Distribution). You can use $N$=15:

$$y'_{n}=\frac{1}{N}\sum_{i=1}^N x_{n}^i$$

with $y=\sqrt{12N}(y'-\frac{1}{2}) \approx N(0,1)$

For randomizing, just use a clock or cpu time seed as $x_0'$.


Your understanding in not correct now. And it can be seen from the equation which states the relation between DTFTs of $h[n]$ and $h_{new}[n] = h[2n]$ as follows :
$$ H_{new}(e^{j\omega}) = \frac{1}{2} \left( H(e^{j\frac{\omega}{2}}) + H(e^{j\frac{\omega-2\pi}{2}}) \right) $$ The general form of this relation for $h_{new} = h[Mn]$ is

$$H_{new}(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} H(e^{j\frac{\omega+2\pi k}{M}}) $$

It's up to you to sketch the resulting expanded and shifted spectrums that add up to create the final spectrum, to see that it's a lowpass filter with cutoff frequency of $\omega_c = \pi/2$

The following MATLAB/Octave code also helps to see it practically:

h = fir1(256, 0.75, 'high');    % Create a HIHPASS filter with cutoff wc = 3*pi/4
h2 = h(1:2:end);                % Create h_new[n] = h[2n] , downsampled filter

figure, plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
hold on
plot(linspace(-1,1,1024), abs(fftshift(fft(h2,1024))),'r');
title('DTFT magnitudes for h[n] and h2[n] = h[2n]');
legend('h[n]','h[2n]')

The resulting plot is:


Metric Spaces: Why $L_\infty$ selects the maximum value

In the general case, let  $x = (x_1,\dots,x_n)$ be a finite-length vector (in a finite dimensional space).
The finite sequence of absolute values $|x_i|$ does attain its maximum (because the sequence is finite), denoted $M = \max_i |x_i|$.

Let $m$ be the (exact) number of coordinates in $x = (x_1,\dots,x_n)$ whose absolute value  is equal to $M$. Thus, $1\le m\le n$. Then, we can lower and upper bound the $\ell_p$ norm of $x$ as follows: $$(mM^p)^{\frac{1}{p}}\le\ell_p(x)\le (nM^p)^{\frac{1}{p}}\,.$$

Both $m^\frac{1}{p}$ and $n^\frac{1}{p}$ tend to $1$ as $p\to\infty$, thus $\ell_p(x) \to M$, by the squeeze theorem.

[EDIT] To provide more concrete substance, let us see what happen with your example: $x=[1,2,3]^T$ (the transposition does not change the result): $$\|x\|_p = d_p(x,0) = (1^p+2^p+3^p)^{1/p} = 3\times\left(\left(\frac{1}{3}\right)^p+\left(\frac{2}{3}\right)^p+1\right)^{1/p}\,.$$

Both $\left(\frac{1}{3}\right)^p$ and $\left(\frac{2}{3}\right)^p$ tend to $0$ as $p \to \infty$. Thus, $\|x\|_p \to 3$.


Does decoding jpeg images uncompresses them?

Image processing is mostly done on frames. The digital image frame is a rectangular raster of pixels, either in an RGB color space or a color space such as YCbCr. So, as you noticed, you can be sure that your compressed images are decoded to rectangular rasters in your application. Each encoded image file stores with its data the source bitmap size (width and height) and the source color space (pixel format and bitness). When decoding, these parameters are restored. So, irrespective of what lossless/lossy algorithms are used, the sizes (in bytes) of recovered raster rectangles are identical, if both lo-res and hi-res compressed image files are received by encoding the same source bitmap data.
EDIT: A "multiresolution" comment fairly corrects my answer, and there is the other explanation why the application may not allocate in CPU memory the entire memory occupied by a source raster: a decoder can use GPU memory. But OP notices that decompressed images take "pretty much the same amount of space in memory" irrespective of encoded image quality. Even with multiresolution, the application can decode the image to source dimensions. Also, code performance dependence on image quality is not part of the question.
Still, the "multiresolution" comment is absolutely necessary for completeness.


Interpolation algorithms in computer graphics

Some ideas:

rotate the image in 2D
compute the Radon transform.  This is used in tomography applications.
Other 2D effects, like swirling or adding ripples
Map the 2D image onto the surface a 3D shape, such as a sphere
Simulate the effect of a lens, or refraction by some other transparent 3D object in front of your bitmap
Simulate gravitational lensing by a black hole

2D image rotation is what comes to mind first as the simplest application of 2D interpolation after rescaling.  The obvious approach is to simply apply the well-known coordinate transformation; an alternative trick is to decompose the rotation into a sequence of three sheers so that only 1D interpolation is needed.

The 3D examples are more involved, since they will probably require implementation of raytracing.


As I've said in the (copious) comments, the signal modeled using

$$
y[t]=A_1(\sin(\omega_1 t+\phi_1)+A_2 \sin(\omega_2 t+\phi_2)+....+A_p \sin(\omega_p t+\phi_p) + z(t)
$$

is stationary.

Your statement i am not interested estimation of parameters,just i am concerned about spectral resolution seems like an oxymoron.  I think what you mean is that you're not concerned with estimating $A_k$ and $\phi_k$, but you are interested in $\omega_k$ for the case where $\omega_k \approx \omega_{k+1}$.

Can you please refine your question, based on this information and the discussion in the comments above?


For large number of samples both will be indistinguishable.

The Biased Version is the Maximum Likelihood Estmator (MLE) of the problem.
It means it has many nice properties for $ N \to \infty $.
The Unbiased version is Minimum Variance Unbiased Estimator (MVUE) of the problem.
Yet in practice when $ N $ is large, as in your case, it won't matter much and there won't be any difference between the two.

Yet, you should remember those are for the case the mean isn't known.
In case it is known, one should always use the Unbiased Form.


Normalizing loudness level of a set of songs

You can spare yourself time by using off-the-shelf software such as http://mp3gain.sourceforge.net/


Separation of sinusoidal oscillations from the frequency spectrum

You don't have to perform a full FFT. Take the definition of DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}$$
What you need is one element of the sequence $X[k]$. For example, suppose the sampling rate is $f_s$, and the signal has a length of $N$, so the frequency resolution is $\Delta f=f_s/N$, and $k$ corresponding to 1000 Hz equals to $1000/\Delta f + 1$. Just do DFT at that discrete frequency $k$.
On the other hand, if you want to convolve your signal with something, use sinusoidal wave of 1000 Hz. We know that a convolution in the time domain equals a multiplication in the frequency domain. The Fourier transform of sine wave is delta function at the specific frequency. So convolving with sine wave can extract the harmonic component.


How to decide what platform is best to implement real-time audio processing on?

If you are keeping to "standard algorithms" like IIR, FIR, radix-2 or 4 FFT (ie stuff that fits DSP architectures well without much control flow), you can try this:

Count up how many "multiply accumulates" you need per second in all your algorithms.

< 10 million you can probably get a fast microcontroller to do the job (or even a slow one if you are <1M)
< 100M is easy DSP territory
< 1G is likely to be doable fast DSP territory
1G-10G is where the cross-over between DSP and FPGA comes
> 10G is multiple DSPs or FPGA


Relating bandwidth and wavelength of a carrier

I will try to derive your formula.  The formula you were given is technically not correct, but on a practical level it is.  I will try to explain below.

The formula for the frequency of any wave (sound, water, light, etc.) is $f = \frac{\nu}{\lambda}$, where $\lambda$ is the wavelength and $\nu$ is the wave's velocity.  This makes intuitive sense if you think about it because you will see more peaks the faster the wave travels, and you will see fewer peaks the longer the wavelength is.  In your case, of course, the wave velocity is $c$.

Bandwidth is the highest frequency minus the lowest frequency.

$$
\begin{align}
B &= f_h - f_l \\
&= \frac{c}{\lambda_l} - \frac{c}{\lambda_h} \\
&= \frac{c}{\lambda - \frac{\Delta\lambda}{2}} - \frac{c}{\lambda + \frac{\Delta\lambda}{2}}
\end{align}
$$
At this point to make the typing a little easier I will call $\frac{\Delta\lambda}{2}$ "$d$".
$$
\begin{align}
B &= \frac{c}{\lambda - d} - \frac{c}{\lambda + d} \\
&= \frac{c(\lambda + d) - c(\lambda - d)}{(\lambda - d)(\lambda + d)} \\
&= \frac{2cd}{\lambda^2 - d^2} \\
&= \frac{c\Delta\lambda}{\lambda^2 - \frac{\Delta\lambda^2}{4}}
\end{align}
$$
The formula you were given is, as a practical matter, correct, because usually $\Delta\lambda$ is much smaller than $\lambda$.  In your case $\Delta\lambda$ is 1/10 of $\lambda$, which is actually much bigger than any normal situation.  Usually $\Delta\lambda$ is orders of magnitude smaller than $\lambda$.  Even at 1/10, though, $\frac{\Delta\lambda^2}{4}$ is $\frac{\lambda^2}{400}$, which is negligible.  Thus, it is dropped.

I hope this helps.


Why does MATLAB's dualtree3 function not return the lowpass subband of the imaginary tree?

Prof. Nick Kingsbury kindly provide an answer to my question!.

In 1-D, the lowpass basis functions (scaling functions) from the two
trees (a) and (b) of the dual tree WT tend to look very similar to
each other, apart from a shift of half the output sample period
between them. Hence it usually makes most sense to regard the two sets
of lowpass samples as interleaved samples from a single 2x
over-sampled real basis function.
On the other hand, the highpass basis functions (wavelets) from (a)
and (b) look approximately like a Hilbert pair of functions. so they
are most logically interpreted as the real and imaginary parts of an
analytic basis function.
Hence the over-sampling by 2 in the lowpass channel replaces the
spitting into real + imaginary channels for the bandpass channels of a
multi-scale wavelet transform.  Similar things apply when 4 trees are
used for the 2-D dual-tree WT, or when 8 trees are used for the 3-D
version.
I hope this helps.  Do refer to our 2005 paper in IEEE Signal
Processing Magazine for more details.

[UPDATE at 9/18/2020]
Since I quoted Prof. Kingsubury's answer above, I should also note one more thing. He mentioned in his email that he is not familiar with dualtree3 function in MATLAB. However, based on my question, he gave the above general answer.


How to include phase in a sinusoidal Kalman Filter

The answer is that the filter I detailed can correct for the phase, the issue was that the values in my Q matrix were too small and therefore it was overly trusting of the model and was correcting for the change in phase very slowly. By increasing the values in the Q matrix this model can quickly correct for the change in the phase of the oscillator!

The output now looks like the following:


Would subtracting the phases of two images be a superior difference metric than subtracting the images directly?

You don't have to work with FFT's and their complex phases, it is far simpler to work directly with the data.  It is particularly simple if the two images are registered.  At the simplest level you find the edges in both of the images and take the difference between the edge maps.  All structure that is common to both images is then zero and all that is left is the pixels that have changed.  This is robust against changes in illumination provided you match the histograms of the edge images prior to subtraction.

Inevitably there are complications, all of which in any case are present when working in the Fourier domain.  I have already touched on the effect of illumination, you can make that slightly more effective by renormalising the histograms of the images before applying the edge filter.  As for which edge filter to use, Roberts, Sobel, or whatever, I would suggest taking the 3x3 Laplacian (this is equivalent to finding the trace of the Hessian of the image density field).

One serious complication is noise: it is essential to remove the noise in the images before starting out.  A simple 3x3 2D smoothing filter passed over the image is fine.  (In the Fourier domain you kill the highest frequencies).  This is important because differentiating random data, which is what edge filters are about, is generally a bad thing to do, and taking second derivatives (the Laplacian) is even worse.  Any noise filter applied in the spatial domain will require a threshold filtering to eliminate noise.  Choosing the threshold visually is effective, but not practical if there is a lot of data.  Then there are auto threshold mechanisms available, and some of those are adaptive to the local noise.  It can get as complicated as you like!

A further complication arises from knowing a priori what is the scale on which you expect to find the difference: how many pixels?  Simple 3x3 filters are inadequate if the features delineated by edges on a scale of >10 pixels.  If you are searching for a particular scale it's worth smoothly shrinking the images so as to make the edges of interest conform to the sought after features.

If you have no idea what the scale is, then my advice would be to use a simple wavelet transform (like the 5-3 wavelet) which automatically locates edges on all scales. I don't have a ready reference to that, but I'll find one if that is of interest.

It can get quite involved depending on how difficult it is to define what you are looking for.  The Multi-scale Morphology Filter was introduced in my field in the paper http://arxiv.org/abs/0705.2072 and refined in a later astronomy paper, http://arxiv.org/abs/1209.2043 .  These dealt with the problem of finding and classifying structure on a point set distribution (galaxies) and, in terms of the question, describe the scale free analysis used to find edges (filaments of galaxies) of arbitrary scale.  The work was also focussed on making the search for structures parameter free, i.e. "objective" and since then has progressed to being able to handle data in N-dimensions of billions of data points (not finished and so not published as yet).

I also use this in industrial process control using thermal imaging cameras, which have weird noise characteristics.

The bottom line is to assess the complexity of your data by doing the very simplest thing contained in the first paragraph.  If that works then don't go any further.


How are signal bandwidth and MSPS related?

The need for modulation:

Your voice signals generally lies in the frequency range 1kHz to 4kHz and music lies in the
range 20Hz to 20Khz. Say your locality has 4-5 AM broadcast stations. Now all the AM broadcast stations cannot transmit their content as is because if all the stations transmit at the same time there will be signal interference. Also, for efficient radiation of electromagnetic energy the radiating antenna should be of the order of the wavelength of signal radiated. If the signal is 20Hz-20kHz the length of antenna has to be of the size of kilometers. So, they place the information carried by 20Hz-20kHz at higher frequencies. Every station places its their content in their own band of frequencies. Like one AM station will operate on 100Mhz with their band of operation being 100Mhz - 20Hz and 100Mhz+20Khz and the second one at 105Mhz with their band of operation being 105Mhz - 20Hz and 105Mhz+20Khz and so and so forth. The process of placing the frequency content from one spectrum on to another is called modulations. Now this information is radiated in to free space by an antenna at the station.

Receiver:

The receiver has a tuning knob, which lets you select a band of frequencies. Say, you have tuned your receiver to receive station 3 which is operating at 110 Mhz. You will select all the frequencies  in the range 110Mhz - 20Hz and 110Mhz+20Khz. The signal being present in the spectrum 110Mhz - 20Hz and 110Mhz+20Khz is brought back to 20 Hz to 20kHz using a process which reverses the operation in encoder side. This is called demodulation. This signal is ready to hear.

For a simple modulation-demodulation read this

Finally, if you want to post-process the signal before listening like reducing the noise, increasing bass/treble, you have to filter the signal.

Software-Defined Radio:

In a software defined radio, there will be a RF receiver which does the job of the receiver mentioned above. The message signal that was recovered is then converted to digital samples using an analog-to-digital converter. The analog to digital converter has 2 steps, 1: Sampling at rate above or equal to Nyquist and 2. Quantization. These digital samples are processed using a software on a pc/embedded system. You can apply filters on these samples. These digital/discrete samples is then converted to analog using a D/A Converter. Main step of a D/A converter is interpolation.

Theoretically, in an ideal system, it is sufficient to sample the signal at Nyquist rate. Before play back, an ideal low pass filter will recover the message/baseband signal with out any information loss. But such a design is not feasible.

At a low sampling rate, the interpolation filter (low-pass filter) will not be able to interpolate intermediate samples effectively. However if  the sampling rate is increased you can have better reconstructed message signal. Hence at higher msps you must receive a better quality signal than at low msps. However, the quality you perceive will become negligible after certain point as your ear might not be able to negotiate the difference.

I believe lowering the msps below Nyquist can introduce cross talk. That doesn't mean that you will receive neighbor radio stations because they got clipped at the RF receiver stage. The cross talk you receive is because of frequency warping.

In the above description i am speculating  that sampling happens on the demodulated signal not on the modulated signal. Because demodulation can be easily achieved via an envelope detector for AM or a phase locked loop for FM. Demodulation in digital domain can be costly as sampling a high frequency signal at nyquist rate and operating at such high bandwidths can be resource consuming.

Also no FM station operates at Ghz. Giga is for satellite. 100 Mhz for FM. and Khz for AM.


Identify if a Specific Edge Pixel Belongs to Horizontal, Vertical or Diagonal Edge Line

Not sure how you do it in python, but the idea can be as simple as the following:

After you select a random point on the edge pixel, calculate the horizontal gradient and the vertical gradient with this pixel at the centre of the gradient filter.

The width of the  kernel for the gradient filters will depend on the maximum width of the edges. For example if the edges are 3 pixel wide have a 5x5 gradient filter.

If the magnitude of vertical gradient is more than horizontal gradient then the pixel belongs to horizontal edge and vice versa for vertical edge.

Also for diagonal edge similarly you can have a kernel which calculated diagonal gradients, if diagonal gradients are not exactly 45 degrees, you might need to to increase the size of the kernel to get more angular resolution of gradients.

Hope that helps.


Can addition of two band limited signals create aliasing?

Ideally no, adding 2 signals should result as a simple sum.

provided you don’t clip

DSP systems have finite dynamic range. A double float has a very large dynamic range so it is not usually a consideration.

Fixed point arithmetic has uniform quantization but the finite dynamic range requires attention.

Floating point is not a uniform quantization, so arithmetic also introduce rounding, which are discontinuous on small scales.


Band edge filter for raised cosine impulse

It looks like the filter you want is indeed $-ith(t)$.  Here is some Octave code to get a visualization in the frequency domain:

t = [-1:0.01:1];
h = 0.5*(1+cos(pi*t));
hd = -i*t*0.5.*(1+cos(pi*t));
H=fftshift(fft(h,512));
HD=fftshift(fft(hd,512));
v = [-256:255];
plot(v, 20*log10(abs(H/512)), v, 20*log10(abs(HD/512)))

BTW, when plotting $h(t)$ in the time domain, your provided $h(t)$ appears to be normalized, spanning from $t = -1$ to $t = 1$.  So, I'm guessing your normalized symbol period is $T = 1$ with some ISI, or $T = 2$ with no ISI.


Why convolution is required, or what is the philosophy behind convolution?

The Idea of Convolution

My favorite exposition of the topic is in one of Brad Osgood's lectures on the Fourier Transform. The discussion of convolution begins around 36:00, but the whole lecture has additional context that's worth watching.

The basic idea is that, when you define something like the Fourier Transform, rather than working directly with the definition all the time, it's useful to derive higher-level properties that simplify calculations. For example, one such property is that the transform of the sum of two functions is equal to the sum of the transforms, i.e.

$$
F\{f + g\} = F\{f\} + F\{g\}.
$$

That means if you have a function with an unknown transform, and it can be decomposed as a sum of functions with known transforms, you basically get the answer for free.

Now, since we have an identity for the sum of two transforms, it's a natural question to ask what the identity is for the product of two transforms, i.e.

$$
F\{f\}F\{g\} = \space ?.
$$

It turns out that when you calculate the answer, convolution is what appears. The entire derivation is given in the video, and since your question is mostly conceptual, I won't recapitulate it here.

The implication of approaching convolution in this way is that it's an intrinsic part of the way the Laplace Transform (of which the Fourier Tranform is a special case) turns linear constant-coefficient ordinary differential equations (LCCODE's) into algebraic equations. That fact that such a transform is available to make LCCODE's analytically tractable is a large part of the reason why they're studied in signal processing. For example, to quote Oppenheim and Schafer:

  Because they are relatively easy to characterize mathematically and
  because they can be designed to perform useful signal processing
  functions, the class of linear shift-invariant systems will be studied
  extensively.

So one answer to the question, is that if you're using transform methods to analyze and/or synthesize LTI systems, sooner or later, convolution will arise (either implicitly or explicitly). Note that this approach to introducing convolution is very standard in the context of differential equations. For example, see this MIT lecture by Arthur Mattuck. Most presentations either present the convolution integral without comment, then derive its properties (i.e. pull it out of a hat), or hem and haw about the strange form of the integral, talk about flipping and dragging, time-reversal, etc., etc.

The reason I like Prof. Osgood's approach is that it avoids all that tsouris, as well as providing, in my opinion, deep insight into how mathematicians probably arrived at the idea in the first place. And I quote:

  I said, "Is there a way of combining F and G in the time domain, so
  that in the frequency domain the spectrums multiply, the Fourier
  transforms multiply?" And the answer is, yes, there is, by this
  complicated integral. It's not so obvious. You wouldn't get out of bed
  in the morning and write this down, and expect that this was going to
  solve that problem. How do we get it? You said, suppose the problem is
  solved, see what has to happen, and then we'd have to recognize when
  it's time to declare victory. And it's time to declare victory.

  Now, being an obnoxious mathematician, you cover your tracks and say,
  "Well, I'm simply going to define the convolution of two functions by
  this formula."

LTI Systems

In most DSP texts, convolution is usually introduced in a different way (that avoids any reference to transform methods). By expressing an arbitrary input signal $x(n)$ as a sum of scaled and shifted unit impulses,

$$
x(n) = \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}, \tag{1}
$$

where

$$
\delta(n) = \begin{cases}
0, & n \ne 0 \\
1, & n = 0,
\end{cases} \tag{2}
$$

the defining properties of linear time-invariant systems lead directly a convolution sum involving the impulse response $h(n) = L[\space \delta(n) \space]$. If the system defined by an LTI operator $L$ is expressed as $y(n) = L[\space x(n)\space]$, then by applying the repective properties, namely linearity

$$
\overbrace{L[\space ax_1(n) + bx_2(n)\space ]}^\text{Transform of the sum of scaled inputs} = \overbrace{aL[\space x_1(n)\space ] + bL[\space x_2(n)\space ]}^\text{Sum of scaled transforms} , \tag{3}
$$

and time/shift invariance

$$
L[\space x(n) \space] = y(n) \space \xrightarrow{\text{implies}} L[\space x(n-k) \space] = y(n-k), \tag{4}
$$

the system can be rewritten as

$$
y(n) = \overbrace{L\left[ \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}\right]}^\text{Tranform of the sum of scaled inputs} = \overbrace{\sum_{k=-\infty}^\infty{x(k)L[\delta(n - k)]}}^\text{Sum of scaled transforms} = \overbrace{\sum_{k=-\infty}^\infty{x(k)h(n-k) }. }^\text{Convolution with the impulse response}
$$

That's a very standard way to present convolution, and it's a perfectly elegant and useful way to go about it. Similar derivations can be found in Oppenheim and Schafer, Proakis and Manolakis, Rabiner and Gold, and I'm sure many others. Some deeper insight [that goes further than the standard introductions] is given by Dilip in his excellent answer here.

Note, however, that this derivation is somewhat of a magic trick. Taking another look at how the signal is decomposed in $\left(1\right)$, we can see that it's already in the form of a convolution. If

$$
\overbrace{(f * g)(n)}^\text{f convolved with g} = \sum_{k=-\infty}^\infty{f(k)g(n-k)},
$$

then $\left(1\right)$ is just $x * \delta$. Because the delta function is the identity element for convolution, saying any signal can be expressed in that form is a lot like saying any number $n$ can be expressed as $n + 0$ or $n \times 1$. Now, choosing to describe signals that way is brilliant because it leads directly to the idea of an impulse response--it's just that the idea of convolution is already "baked in" to the decomposition of the signal.

From this perspective, convolution is intrinsically related to the idea of a delta function (i.e. it's a binary operation that has the delta function as its identity element). Even without considering its relation to convolution, the description of the signal depends crucially on the idea of the delta function. So the question then becomes, where did we get the idea for the delta function in the first place? As far as I can tell, it goes at least as far back as Fourier's paper on the Analytical Theory of Heat, where it appears implicitly. One source for further information is this paper on Origin and History of Convolution by Alejandro Domínguez.

Now, those are the two of the main approaches to the idea in the context of linear systems theory. One favors analytical insight, and the other favors numerical solution. I think both are useful for a full picture of the importance of convolution. However, in the discrete case, neglecting linear systems entirely, there's a sense in which convolution is a much older idea.

Polynomial Multiplication

One good presentation of the idea that discrete convolution is just polynomial multiplication is given by Gilbert Strang in this lecture starting around 5:46. From that perspective, the idea goes all the way back to the introduction of positional number systems (which represent numbers implicitly as polynomials). Because the Z-transform represents signals as polynomials in z, convolution will arise in that context as well--even if the Z-transform is defined formally as a delay operator without recourse to complex analysis and/or as a special case of the Laplace Transform.


Facing problem in plotting histogram of DCT coefficient of an image

My guess is that the assertion

  "the histogram of DCT coefficients of the un-compressed natural image
  has the smooth histogram ... However, the histogram of DCT
  coefficients of the compressed natural image has peaks and gap in it."

is valid when one uses a compression method that uses a DCT (Discrete Cosine Transform), and quantizes (sufficiently) its coefficients. This might not happen for compression methods that differ from JPEG/MPEG principles, that use a 2D DCT. Moreover, JPEG applies DCT2 on $8\times 8$ blocks, and the DC coefficient undergoes a predictive compression scheme different than that of the other AC coefficients.

So I am not surprised that you don't observe quantization when you apply the DCT on the whole image instead of blocks (I obtained the same results).

What could be more successful would be to compare the DCT coefficient histograms (DC coefficient removed) on each of the  $8\times 8$ blocks, raw and compressed, and then combining the coefficients of all the blocks. Here is an attempt, in Matlab.

As you can see, you now have holes in the right-side histograms, due to quantization.

nBin = 512;
dctCoeffLimit = 80;

img_1 = double(rgb2gray(imread('lena.bmp')));
img_2 = double(rgb2gray(imread('lena.jpg')));

% Compute block DCT with mean removed
funDCT2 = @(block_struct) dct2(block_struct.data-mean(mean((block_struct.data))));
img_1DCT= (blockproc(img_1,[8 8],funDCT2));
img_2DCT= (blockproc(img_2,[8 8],funDCT2));

% Set DC coefficients to NaN
img_1DCT(1:8:end,1:8:end) = nan;
img_2DCT(1:8:end,1:8:end) = nan;

% Convert 2D to 1D
img_1DCTflat = img_1DCT(:);
img_2DCTflat = img_2DCT(:);

% Remove DC (NaN) values
img_1DCTflat(isnan(img_1DCTflat)) = [];
img_2DCTflat(isnan(img_2DCTflat)) = [];

subplot(1,2,1)
histogram(img_1DCTflat,nBin);axis([-dctCoeffLimit dctCoeffLimit 0 Inf]);
subplot(1,2,2)
histogram(img_2DCTflat,nBin);axis([-dctCoeffLimit dctCoeffLimit 0 Inf]);


How can I decide proper FFT length(size)?

But I want to know more about 21,53 Hz. exactly FR.

This represents the bandwidth (how wide in Hz a bin is) of a particular bin. The smaller this value is, the better FR you have. In order  to get the bin center frequency, you can use the following formula

bin_center_freq = FR * bin_number;

So, for, say, bin at index 20 , we get

bin_center_freq =   21,53*20 = 430,66 Hz (Note that the bin at index 0 represents DC (0 Hz)

When we have a bin center frequency, we can easily calculate the lower and upper frequencies for that bin.

bin_lower_freq = bin_center_freq-(FR/2) = 419,89 Hz

and

bin_upper_frequency = bin_center_freq-(FR/2) = 441,42 Hz

Note that the real frequency resolution (our ability to separate close-by frequencies) will still depend on the type of the window used for the FFT (Blackman, Gaussian, Hamming, rectangular etc) as well as on the spectral content present (whether we have components that are both close in frequency and amplitude etc) in the signal so you shouldn't assume
that a particular bin contains only frequencies belonging to that bin(due to the spectral leakage, every bin will also have energy contributions from neighbouring bins).

So, what is the right FFT size?

It really depends on what you want to do with your signal. If you have a signal containing 2 sine waves close in frequency and amplitude (e.g. one 220 Hz, and the other 225Hz), you should choose a relatively long FFT length such that
the FR is less then 5 (f1-f2). But since windowing also affects the result, make your FFT length even longer to compensate for the windowing effect. You can read up more on window functions here ( https://en.wikipedia.org/wiki/Window_function )


Acoustic waves propagation distance relation with carrier frequency?

For a sine wave of amplitude $A$ and frequency $f=1/T$, its energy from $-T/2$ to $T/2$ is $A^2T/2$. Thus the power depends only on amplitude. The energy over one period does depend on frequency, but in the same observation duration, two sines having same $A$ should have the same energy (assuming observation duration is multiple of their periods).

In communications, we don't transfer energy but signal. Then power, not energy, matters.

Propagation loss is characterized grossly by distance (pathloss model) and absorption that the latter explains the choice of low frequency in underwater acoustics communications. For example, frequencies >= 1 MHz are absorbed very quickly whereas the absorption of low frequencies is much weaker. This is physics and the explaination would be complicated.

FYI, you can use this calculator for sea absorption.

And no, $E=hf$ is for photon, I have not met people who uses quantum mechanics for underwater acoustic waves.


Values of angular frequency of discrete-time sinusoidal signal

Because the sequence

$$ x[n] = A \cos( \omega n) $$

will produce exactly the same samples for all discrete-time frequencies such as $\omega = \pi/3$ or $\omega = \pi/3 + 2\pi$ or $\omega = \pi/3 + 6\pi$ or $\omega = \pi/3 -8 \pi$... In other words every frequency of the form $w = w_0 + 2\pi k $ for all integers $k$ will yield the same set of samples; they are indistinguishable, and we chose to use the range $[0,2\pi]$ or $[-\pi,\pi]$, when referring to the discrete-time sinusoidal frequencies.

Note that the periodicity of the discrete-time frequency $\omega$ relies on the fact that the index $n$ is an integer. Otrherwise, like in the continuous-time case, the frequency $\omega$ would range from $-\infty$ to $\infty$.


Group delay of non differentiable unwrapped phase

If there are phase discontinuities after all $2\pi$-jumps have been removed, then these discontinuities are usually caused by zeros of the frequency response. The phase jumps by $\pi$ at these frequencies, and the group delay doesn't exist, or, if one prefers, is non-finite. Note that the group delay is meaningless anyway at frequencies where the frequency response is zero.


What are the differences between the DWT and the MODWT?

from "wavelet methods for time series analysis"


how do I derive the system equation for a simple delay with feedback?

In such cases it helps to define an additional signal $w[n]$ at the input of the delay. Now you can write down two equations describing the system:

$$\begin{align}w[n]&=x[n]+g\cdot w[n-N]\\y[n]&=x[n]+w[n-N]\end{align}\tag{1}$$

where $N$ denotes the delay (in samples) and $g$ is the feedback gain. From the second equation you can express $w[n-N]$ in terms of $x[n]$ and $y[n]$:

$$w[n-N]=y[n]-x[n]\tag{2}$$

Plugging this into the first equation gives

$$y[n+N]-x[n+N]=x[n]+g\cdot \left(y[n]-x[n]\right)\tag{3}$$

Rearranging and subtracting $N$ from all indices results in

$$y[n]=x[n]+(1-g)x[n-N]+g\cdot y[n-N]\tag{4}$$

Eq. $(4)$ is the (single) difference equation describing the system. However, it's more efficient to implement the system using the two equations given in $(1)$.


Matlab, like most computational packages, implement a pretty neat set of fast Fourier transform algorithms referred to as FFTW: fastest Fourier transform in the west. You can find more out about it here.

FFTW actually runs a preflight on your machine when it is first installed. At that time, it runs a bunch of test computations to figure out which prime factorization would be fastest given your processor speeds and cache sizes. Then, it uses a symbolic program to generate a bunch of FFT algorithms specific to your machine for a whole bunch of cases, like real valued, symmetric, etc. So, it's hard to say exactly what's going on under the hood except that if there's an efficiency to be had, FFTW will take advantage.


IFFT of a complex, conjugate-symmetric transfer function - Value at Nyquist frequency

If you have an exactly conjugate-symmetric frequency domain vector, the values  at DC and at Nyquist must be real-valued. The 'symmetric' flag of Matlab's ifft command makes sure that the result of the inverse FFT is real-valued, implying that the values at DC and Nyquist must also be real-valued.

Let $X[k]$ be the length $N$ DFT of a real-valued sequence $x[n]$:

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

The value at DC is

$$X[0]=\sum_{n=0}^{N-1}x[k]\tag{2}$$

which is clearly real-valued if $x[k]$ is real-valued. Similarly, assuming that $N$ is even, the value at Nyquist is

$$X[N/2]=\sum_{n=0}^{N-1}x[k]e^{-j\pi n}=\sum_{n=0}^{N-1}x[k](-1)^n\tag{3}$$

which is also real-valued for real-valued $x[n]$.

Conjugate symmetry means

$$X[k]=X^*[N-k]\tag{4}$$

which for $k=N/2$ (i.e., at Nyquist) gives

$$X[N/2]=X^*[N/2]\tag{5}$$

Of course, $(5)$ just shows in another way that the bin at Nyquist is real-valued for real-valued $x[k]$.


Pattern classifier experiment for determining sex from face image of persons using kNN?

At your first stage, which classifier to use is not that important, you may need to make your code run first. There was a very famous paper on Face Recognition and Gender Determination, in which a very complex transform was implemented to represent each node on the face. The used gabor kernal based wavelet transform to get set of coefficients for each node(refers to jet), and use sparse jets to represent the face. They also consider the template matching to ensure that each face can be decomposed of the jets. They built a look-up table to record the correlation between each jet and the gender. When a new image comes, the gender is determined by the pooling result of the jet sets.

But most approaches nowadays won't be that complicated. Basically the image intensity is used as the face input (or RGB color components but they are intensity equivalently). Each image matrix is reshaped into a long vector, and all the images from training set will form a large M X N matrix where M is the image pixel size and N is the sample number. Then Principal Component Analysis is used on the matrix to extract the core components, and project them onto a sub-space. Intra-class and inter-class scatter difference is calculated, then Linear Discriminant Analysis is applied to maximize the ratio of the determinant of Inter-class scatter and intra-class scatter. Please refer to this link to see more details (The algorithm is implemented by Python though).


Guess stability based on filter coefficients

In the general case there is no simple way to determine the stability of a discrete-time system (unless it is FIR, in which case it is always stable). In order to check stability you don't necessarily need to compute the poles of the transfer function, but you can use the Jury stability criterion, which does not explicitly compute the poles, because the latter is generally a numerically ill-conditioned problem.

For a second-order system the denominator coefficients have to be inside the stability triangle, which can be checked very easily.


Since your offset is 3/5ths (plus some integer), the truncation fractions will be different, resulting in a different sum of quantization noise.  You can get rid of this difference by using an integer offset (including/after scaling) before truncation or rounding.  (e.g. not (2^N)/5)

Also, if your signal is not exactly periodic in your FFT aperture, then the result magnitude of very low frequencies will have a small dependence on phase, relative to the window (edge or center), and thus not exactly match the input amplitude.


A simple example: consider a source $S = \lbrace 0, 1 \rbrace$ with probabilities $\lbrace 0.9, 0.1\rbrace$.

If you take one symbol at a time, no compression is possible: you need one bit per symbol.

Taking two symbols at a time, compression is possible. Now you have a source $S_2 = \lbrace 00, 01, 10, 11 \rbrace$ with probabilities $\lbrace 0.81, 0.09, 0.09, 0.01 \rbrace$. A possible encoding is:

Message | Codeword
--------|---------
00      | 0
01      | 10
10      | 110
11      | 111

Here, the average number of bits per message from source $S$ is $(1 \cdot 0.81 + 2 \cdot 0.09 + 3 \cdot 0.09 + 3 \cdot 0.01)/2 = 0.645$. It is clear, however, that by having only four combinations, the prefix code is not as short as it could be -- for example, $01$ and $10$ should have codewords of the same length, which in turn should be shorter than the codeword for $11$. If you group more $S$ messages, then you have more room for optimizing the code.

Note that the entropy of $S$ is around $0.47$; this means that the code above is still far from optimal. What Shannon says is that one may approach the entropy by increasing the number of symbols that are considered as input to the prefix code.


You would upsample by 441 and downsample by 640 since $11025 \cdot 640 = 16000 \cdot 441$

That would be a rather expensive operation, so in practice you would either fugde it a bit: 31/45 would get you to 11022Hz and 11/16 to 11000 Hz or you would use a suitable poly-phase irrational conversion algorithm


The interpolation lowpass filter computes a weighted sum of input samples, which results in the zero input samples being interpolated using the non-zero samples of the input signal. The $*$ sign is NOT multiplication but convolution!

EDIT: I'm adding a simple example to clear things up a bit. Let's consider the impulse response $h = [0.5, 1, 0.5]$ and a zero-padded input signal
$x=[1,0,3,0,5,0,3,0,1]$. The time index is denoted by $n$, and $n=0$ corresponds to the left-most values of $h$ and $x$. The output $y$ is the convolution of $(x*h)(n)$:

$$y(0) = h(0)x(0) = 0.5\\
y(1) = h(0)x(1) + h(1)x(0) = 1\\
y(2) = h(0)x(2) + h(1)x(1) + h(2)x(0) = 2\\
y(3) = h(0)x(3) + h(1)x(2) + h(2)x(1) = 3\\
y(4) = h(0)x(4) + h(1)x(3) + h(2)x(2) = 4\\
y(5) = h(0)x(5) + h(1)x(4) + h(2)x(3) = 5\\\vdots$$

This impulse response obviously performs linear interpolation of the zero-padded input signal.


Simplification after sampling an analog signal

Just for explanation: What you experience here, is aliasing. Look at this code and its output:

Fs = 80 # the sampling frequency
f = 50 # the frequency of the cos wave
t = np.arange(0, .1, 1./Fs)  # the sample times
t2 = np.arange(0, .1, 1./(10*Fs))  # higher sampling frequency illustration purposes

plt.stem(t, np.cos(2*np.pi*f*t), label='sampled points')  # plot sampled values
plt.plot(t2, np.cos(2*np.pi*f*t2), label='50Hz cos wave')
plt.plot(t2, np.cos(2*np.pi*(Fs-f)*t2), label='30Hz cos wave')
plt.legend()

What I show is a 50Hz cos wave (blue) and the sampling points when sampling with a frequency of 80Hz (blue dots). I also draw a cos wave with frequency 30Hz on top of the diagram (green).

As you can see, at the sampling points the green and blue line intersect. I.e. when your sampled values do not differ, if you sample a 50Hz or 30Hz wave. This is the aliasing effect: You cannot distinguish between waves that have a different frequency, since the sampling rate $F_s$ is too low. In general, given a sampling frequency $F_s$, waves with frequencie $f$ and $F-f$ create the same sampled values.


Conversion of transmitter power from Watts into dB

Thanks to @MBaz who pointed out my mistake: dBm and dBW are unitless, but they relate to  a reference of 1 mwatt, and 1 Watt.

The question has an error, as well as the answer. Watts does not convert to dB, it converts to dBW. (i) should be as follows:

find (i) the transmit power in dBm and dBW.

The reason behind that is that dB is a ratio, while dBm and dBW are ratios with a reference to 1 mWatt and 1 Watt, respectively, and you can use them to recover an absolute value of the power.

For example, the signal to noise ratio is a measure of how strong a signal is compared to the noise. In the linear scale it is written as:

$\frac{\text{Signal Power (in watts or miliwatts)}}{\text{Noise Power (in watts or miliwatts)}} = SNR (\text{unitless})$

or in logarithmic domain:

$\text{Signal Power (in dBW or dBm)}-\text{Noise Power (in dBW or dBm)} = SNR (\text{in dB})$

So the correct answer to (i) is

TX power = $10\log(50)=17\text{ dBW}=17+30=47\text{ dBm}$


Contradiction while using the convolution sum for a non-LTI system

You've some troubles about the LTI property and the impulse response. First of all let me give you the answer to your quiz question.

Given two LTI systems $\mathcal{T}_1$ and $\mathcal{T}_2$ with impulse responses :

$$h_1[n] = \delta[n+3]$$
$$h_2[n] = \delta[5-n] = \delta[n-5]$$

you can compute the output of the cascade $$y[n] = \mathcal{T}_2 \{ \mathcal{T}_1 \{ x[n] \} \} $$ for the input signal $x[n] = n$ as :

$$ y_1[n] = \mathcal{T}_1 \{ x[n] \} \implies y_1[n] = n \star \delta[n+3] = (n+3) $$

and

$$ y_2[n] = \mathcal{T}_2 \{ y_1[n] \} \implies y_2[n] = (n+3) \star \delta[n-5] = (n-2) $$

so your output is $$y[n] = n-2$$

Note that $h_2[n] = \delta[5-n]$ and $h_2[n] = \delta[n-5]$ are the same systems, and yield the same output.

Note again that both systems are declared to be LTI and are given their respective impulse responses $h_1[n]$ and $h_2[n]$ , therefore their cascade will also be LTI and furthermore the output does not depend on the ordering of the systems; i.e.

$$y[n] = \mathcal{T}_2 \{ \mathcal{T}_1 \{ x[n] \} \} = \mathcal{T}_1 \{ \mathcal{T}_2 \{ x[n] \} \}  $$

Now coming to your own questions:

1- You cannot use convolution sum to deduce LTI property. You cannot form a convolution sum if the system is not already LTI.

2- The question seems ok:

Note: your claim that $ h[n] = \delta[5-n] $ , signifies a non-LTI system is wrong...


Real time means that your algorithm can process the video faster than the time it takes to watch the video. In other words, if you were collecting the video live you would be able to finish processing it as soon as you stopped collecting it with no extra work needed to be done after.


Estimate the Filter Coefficients of 1D Filtration (Convolution)

This is a nice question.
I will try solving it using 2 approaches (Which are basically the same).

The solution is the Least Squares Solution:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

We assume data is given in finite discrete form (As it is in practice).
The convolution is done in valid mode (Like in MATLAB valid property).

Direct Solution

One could write the above as:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2}  $$

Basically using the commutative property of the convolution one could the above as $ x $ was the filter and hence build a Convolution Matrix from $ x $.

Now, the solution is the usual Least Squares:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2} = \left( {X}^{T} X \right)^{-1} {X}^{T} y $$

Pay attention that for long signals the matrix is huge.
Hence this method is practical only for small signals.

Iterative Solution

Well, one could use the Gradient Descent Method to minimize the cost function:

$$ f \left( h \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

The Gradient is given by (Correlation is the Adjoint Operator of Convolution):

$$ \frac{d}{d h} f \left( h \right) = \left( h \ast x - y \right) \star h $$

Where $ \star $ stands for Cross Correlation.

In MATLAB Code it is given by:

hObjFun = @(vH) 0.5 * mean( (conv2(vX, vH, 'valid') - vY) .^ 2 );
vG      = conv2(flip(vX, 1), (conv2(vX, vHEst, 'valid') - vY), 'valid');

Where all data is in Column Vector (The conv2() is used as it is faster since conv() is just a wrapper around it).

This method is fast and not limited by the size (Or less, to be accurate).

The full MATLAB Code is in my StackExchange Signal Processing Q18993 GitHub Repository.


What is the effect of noise on Shannon entropy?

The grayscale image is actually a discrete image. entropy calculates a histogram and from that extracts "empiric" probabilities, to be used in the common information formula

$$I= -\log_2 (P(g))$$

with $P(g)$ being the (empiric) probability of an intensity value $g$.

Now, entropy $H$ is the expectation of Information. For your whole image $V$ that means

$$\begin{align}
H(V) &= E\left\{ I(g_{x,y})\right\}\\
&= \sum\limits_{x,y\in \text{idx(V)}} P(g_{x,y})\cdot\left(-\log_2 (P(g_{x,y}))\right)\\
&= -\sum\limits_{x,y\in \text{idx(V)}} P(g_{x,y})\log_2 (P(g_{x,y}))
\end{align}$$

If you look at the product in the sum intensely enough, you'll notice that the sum will take a maximum value if all properties are equal – that is, if every intensity value is just as likely as every other. The histogram of that image would be a flat line.

Adding "well-behaved" noise will make almost any image look more like that – in the end, pixels with identical intensity will be relatively common in most natural images, and adding noise will "break" that.

Now, salt'n'pepper is different: you add extreme values to your image. That will make the intensity at many pixels clip at the maximum intensity (or the minimum one), skewing the histogram away from flat to a very "horned" thing, which has less entropy.


Determining Causality From Discrete Impulse Response

As soon as you've figured out how the given values for $h[n]$ relate to the time index $n$ (which is unclear in the current problem formulation), it should be easy to decide whether the given system is causal or not by considering the input-output relation of a discrete-time LTI system:

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{1}$$

where $x[n]$ is the input sequence, and $y[n]$ is the output sequence.

For the system to be causal, the output at time $n$ must only depend on the current value or on past values of the input $x[n]$. What does that imply for the impulse response $h[n]$?


Does image size impact compression efficiency differently in different formats?

JPEG images can't be that big; their maximum size is (2¹⁶-1)×(2¹⁶-1). (doesn't really matter)
These are fundamentally different file formats, and their compression stems from different properties: Where JPEG is a lossy format who is able to drastically reduce the file size of content that has photographic properties (mostly: low-pass behaviour in 8×8 pixel blocks) at different levels of lossiness, PNG is lossless and compresses especially well on images with similar rows.
Thus, the size of the images depend on the content. You can't directly extrapolate from a small image to a larger version: it depends on the content of the original large image whether it compresses better in JPEG (for a given lossiness factor) or PNG.

Hence, no statement of the sort you seem to be looking for can be made, in general.

Specfically, for images that are large and that results in large swaths of identical pixels, JPEG always operates on blocks of 8×8 pixels. Even if there's a lot of monocoloured blocks, these still have to be stored and later entropy-coded in JPEG; this leads to repeating bit patterns, but not necessarily constant 0-bytes, having to go through a Huffmann-alike compressor.

In PNG, identical runs of pixels will be converted to differences to neighboring pixels; thus, large portions of 0-runs exist in large pictures with constant areas. The DEFLATE algorithm that follows this differential encoding is very efficient in reducing consecutive constant values to very few bits.


recovering phase of sine signal from FFT

It makes sense after all: your phase estimate is an estimate of the unknown phase angle $\phi$ minus $\pi/2$, because you're using a sine function. So the phase of the FFT at the appropriate frequency index corresponds to $\phi-\pi/2$. So in your example, $\phi\approx -0.98+\pi/2=0.59$, which is a lot closer to the actual phase. If you used a cosine function, your way of estimating the phase would be correct.


zero-phasing frequency components while keeping the same magnitude, in Matlab

it is obvious that abs(real(fft(x)))!=abs(fft(x)) since you are removing the imaginary part while taking the real function.One way if you want a signal with all dft components being real.Is to normal to just take abs(fft(x)) as the spectrum itself rather than taking the real components.In this you maintain the magnitude and all the components have zero phase.I hope this answers your question


Sample-rate, filtering, digital-filtering and aliasing

An ideal square pulse - which I assume is a model for your up-and-down deflections - would have an infinite bandwidth, but the bulk of the energy is within a bandwidth of $1/T$, where $T$ is the pulse duration.  Roughly speaking, therefore, the 100 kHz Bessel filter will thus allow you to detect pulses of duration $10 \mu s$ and above.  It will also limit the sharpness of the up-and-down transition.

You want the Bessel filter to prevent aliasing, which means that it needs to attenuate signals in the 500 kHz +/- 100 kHz range, which would alias back into the passband of the filter.  By my calculations, the attenuation of a 4th order Bessel filter would be about 23 dB in this range (but check your datasheet), which isn't great but may be good enough in some applications.  It depends on how high of a signal-to-noise ratio you need, how noisy your signal is within this range, and how much spectral content your deflections have in this range (which will depend on their width and the sharpness of their transitions).

The rule-of-thumb you mentioned about sampling 10x above your filter settings may be specific to a particular filter.  In general, you only need to sample at twice the bandwidth of your signal, which is the well-known Nyquist Criteria.  In practice, we usually sample at higher rates because it allows us to use a more realistic, lower cost filter - such as your 4th order Bessel filter.  The required sampling rate, the filter specifications, and the spectral properties of the signal you are sampling are all related.

Having sampled at 500 kHz, you would want to filter your signal with a 35 kHz filter if you are only interested in spectral content less than 35 kHz.  This also allows you to reduce your sample rate to some value greater than 70 kHz, which could reduce the computational burden of whatever processing you are doing.  But you would not be concerned with the 10x rule-of-thumb anymore, because any aliasing due to the original sampling has already occurred.  You can implement a digital filter that is considerably tighter than the Bessel filter.


why if we pad any signal with enough zeros we can get the same result

Hint:  try computing the linear convolution of $~[1\quad 1\quad  1]~$ and
$~[0\quad 1 \quad 1]~$ by hand and the circular convolution of the zero-padded
versions
$~[1\quad 1\quad  1 \quad 0\quad 0\quad 0]~$ and $~[0\quad 1 \quad 1\quad 0\quad 0\quad 0]~$ also by hand (or via FFTs if you are desperate enough).


If you already know the - quite well-known :) - FT pair $$e^{-at}u(t) \longleftrightarrow \frac{1}{a+j2\pi f}$$ then by setting $a = 1/T$, you get $$e^{-\frac{t}{T}}u(t) \longleftrightarrow \frac{1}{\frac{1}{T} + j2\pi f}$$
Time delay property $$x(t-t_0) \longleftrightarrow X(f)e^{-j2\pi ft_0}$$ can be used now by setting $t_0 = T$, and then $$x(t-T) = e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{1}{\frac{1}{T} + j2\pi f}e^{-j2\pi fT}$$
Adding your $1/T$ constant in front, you get $$\frac{1}{T}e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{\frac{1}{T}}{\frac{1}{T} + j2\pi f}e^{-j2\pi fT}$$
Finally, multiplying both numerator and denominator by $T$ you get what you want: $$\frac{1}{T}e^{-\frac{t-T}{T}}u(t-T) \longleftrightarrow \frac{1}{1 + j2\pi fT}e^{-j2\pi fT}$$


Is there any connection between frequency response of a system and spectrum of a signal?

As mentioned in my comment, look at the Wold Decomposition Theorem :

https://en.wikipedia.org/wiki/Wold%27s_theorem

but actually showing how it relates to your question,

http://www.phdeconomics.sssup.it/documents/Lesson11.pdf.

In simpler terms, your AR system can be represented by a MA system (possibly infinite order) and as the second reference shows, an ARMA system.  The process has more than one representation so there isn't a one to one correspondence between the time series and some system that generated it.

So the answer is yes, there is a correspondence and often the system itself, is of interest, but the time series and system are not uniquely equivalent.

another way of thinking about it is that $H(\omega)$ has a phase and magnitude response that gives $|H(\omega)|^2$ but there are other possible phase responses that results in the same magnitude squared.


What is the gain (advantage) of oversampling and noise shaping in D/A conversion?

If you play 16-bit audio at 48kHz, you need the DAC analog reconstruction filter to pass 20kHz and attenuate 96dB at 24kHz, which is quite steep and requires complex multistage analog filter.
The advantage of using oversampling is moving the sampling rate much higher, for example oversampling by 4x means the DAC runs at 192kHz, and the analog filter only needs to pass 20kHz and block 96dB at 96kHz, which allows for much simpler filter.
The noise shaping is much like dithering after signal processing, but instead of adding white noise uniformly to the whole band, it weights the noise to high frequencies above the audio band, which can then be filtered away by the DAC reconstruction filter, so less noise remains on the audio band, thus increasing the apparent signal-to-noise ratio of the audio band.


Difference Between Contrast Enhancement and Texture Enhancement in Images

Usually Contrast Enhancement is Element / Pixel Wise Operation in Digital Image Processing.

It means it adjusts the value of a pixel based only in its value.
In order to increase contrast we usually push bright values to be brighter and dark values to be darker (S Shaped Contrast).

When people try to enhance Texture they usually refer to the term Local Contrast.
Local contrast means we want to enhance contrast but using a Neighborhood Operation.
In the sense of contrast, we will push a pixel to be darker even if its value is bright but it is darker than its neighbors. Namely we have the locality sense of the operation.
For instance High Pass will increase local contrast.


When to use the DTFT vs the DFT (and their inverses) in analysis?

DFT and DTFT are obviously similar as they both generate the fourier spectrum of time-discrete signals. However, while the DTFT is defined to process an infinitely long signal (sum from -infinity to infinity), the DFT is defined to process a periodic signal (the periodic part being of finite length).

We know that the number of frequency bins in your spectrum is always equal to the number of samples processed, so this also gives a difference in the spectrums they produce: the DFT spectrum is discrete while the DTFT spectrum is continuous (but both are periodic with respect to the Nyquist frequency).

Since it is impossible to process an infinite number of samples the DTFT is of less importance for actual computational processing; it mainly exists for analytical purposes.

The DFT however, with its finite input vector length, is perfectly suitable for processing. The fact that the input signal is supposed to be an excerpt of a periodic signal however is disregarded most of the time: When you transform a DFT-spectrum back to the time-domain you will get the same signal of wich you calculated the spectrum in the first place.

So while it does not matter for the computations you should note that what you are seeing there is not the actual spectrum of your signal. It is the spectrum of a theoretical signal that you would get if you periodically repeated the input vector.

So I would assume in the literature you were mentioning, every time it is important that the spectrum you're working with is actually the spectrum and disregarding the computation side of things, the author would pick the DTFT.


Inverse Z-transform mystic simplification

This is probably homework, so I'll just give you a few hints to get you on the right track. First of all, given an expression $X(z)$, the corresponding sequence is usually not unique. Your solution is correct if the region of convergence (ROC) of $X(z)$ is $|z|>4$, which gives you an exponentially increasing sequence. If this were the impulse response of a discrete-time system, then the system would be causal but unstable.

The other solution assumes that the ROC is $\frac14<|z|<4$, which includes the unit circle. This corresponds to a two-sided decaying sequence. As an impulse response it would correspond to a non-causal but stable system.

So, to answer your question, both solutions are correct but they are not identical, because they are inverse transforms of $X(z)$ assuming different regions of convergence.


In time domain, OFDM breaks one serial fast bit stream into many parallel slow bit streams.
Then, these parallel slow bit streams are multiplied with orthogonal sinusoids, where orthogonality between two sinusoids is defined with a summation over a certain time interval as

$$\sum \limits _{n=0} ^{N-1} \cos 2\pi k f_0 n \cdot \cos 2\pi k' f_0 n = 0$$ when $k \neq k'$.

This process is illustrated in detail in Figure below with the help of an example.

Now for each such individual signal, the multi-path components arrive in not
too distant future (just like a low rate stream). Hence, the
equalizer design is easy having less spread paths and consequently less interference with future symbols, provided that we find a way to separate the subcarriers at the Rx.

Separating these bits at the Rx is easy: we can correlate (multiply sample by sample and sum) the composite signal with just one subcarrier.  Utilizing their orthogonality property, contribution from all other subcarriers will cancel out to zero, while the contribution from the subcarrier with that frequency will pop out, scaled in amplitude by our modulation signal. The same procedure can be repeated for all other subcarriers.

Essentially, this is an operation of Discrete Fourier Transform (DFT) as

$$X[k] =\sum \limits _{n=0} ^{N-1} s[n] e^{-j2\pi \frac{k}{N}n}$$

which, for each $k$, will generate our modulate data.

With respect to frequency domain, the signal at the Rx is a product of the spectra of the Tx signal and the wireless channel. Thus, the channel will allow some frequencies to pass through unharmed while suppressing some others. The low data rate signal needs less manipulation by the Rx to get the original data back due to a narrower spectrum. Essentially for this kind of signal, the channel acts just as a single multiplier that can be equalized through estimating that number and dividing the Rx signal by that number. So the equalization reduces to a single division operation. On the other hand, a high data rate signal needs a lot of Rx processing to equalize the channel.

To solve this problem, what OFDM does in frequency domain is fairly simple. It just segments the available bandwidth into many parallel almost flat channels through utilizing those sinusoidal subcarriers. This slicing of the spectrum is drawn in Figure below, just like slicing a bread.

Hence, equalization for each narrow slice requires just a single division operation, rendering the computational load of the equalizer to a total of $N$ divisions.


There is no "best method", it will vary according to application. There are many trains of thought regarding this subject.

I've done something in that area using Cepstral Analysis. Take a look at this:

https://8b14d84720d60d831e70ec440417feb774595227-www.googledrive.com/host/0B4-fFSwiy1mGdS00SHB2Nlo3bkU/imreg_translation/imreg_translation.html

Note: I am the author of the above.

Basically, you do the following:

Sum the target image and the object
Calculate the Cepstrum of said sum.
Find the peaks (there will be 2, you have to resolve the ambiguity using something like SAD)
The location of the peak will tell you the translation you must apply to the image.

I based it off this article: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6128658&tag=1

I'm actually in the middle of creating a barcode (1D and 2D) open-source program using this kind of processing. The link to my bitbucket acount is here: https://bitbucket.org/DSevero/

Although the repo is private for now, I should be releasing a first version around January

:)


Translation invariance in max-pooling and cascading with convolutional layer

The idea of complex cell layer is to "pool" a set of simple cells and acquire the same data from those simple cells (e.g invariant translation, rotation invariant, etc). One nice illustration is the computational model of Neocognitron  https://youtu.be/Qil4kmvm2Sw?t=6m9s

As for the example of translating 2x2 and 3x3 matrix by one pixel in deeplearning.net, I think one possible explanation is that the Pooling map returns the same value for the subsequent Convolution map. Considering the following 4x4 image with our significant data in 2x2 region (row = 2, col = 2) to (3,3):

0 0 0 0
0 1 6 0
0 2 3 0
0 0 0 0

The Max-Pooling map always look at the region (2,2) and (3,3) and it returns 6. If our data were translated (moved) one pixel, there would be 8 different ways to do that. Among those there are 3 ways that keep region ((2,2); (3,3)) has 6 as the maximum value:

drow = 1, dcol = 0
0 0 0 0
0 0 0 0
0 1 6 0
0 2 3 0

drow = 0, dcol = -1
0 0 0 0
1 6 0 0
2 3 0 0
0 0 0 0

drow = 1, dcol = -1
0 0 0 0
0 0 0 0
1 6 0 0
2 3 0 0

That said, if our data is represented by the maximum value then the recognition will work for all those 3 cases.

When applying the same translation process (by 1 pixel) for 3x3 matrix

1 2 3
4 5 6
7 8 9

The the maximum value at the corner (1, 3, 7, 9) will be retained in 3 translations, the value in the middle of each edge (2, 4, 6, 8) will be retained in 5 translations. The simple cells are responsible for line/edge detection, so probably the center position (5) is left out.

That's how I understand translation invariance. I'm not completely sure though, since the author stated that "[...]produce exactly the same output at the convolutional layer." It would be great if the author of deeplearning.net could explain a bit more on that.


Is it possible to lose half the samples, and yet see the full spectrogram as the original signal?

To rephrase your question: Is there a process that halves the number of samples of a full-bandwidth signal without significantly affecting the appearance of the spectrogram, except for scaling of time or frequency axes?

Resampling does not do that. The only such process I know of is time stretching / pitch scaling followed by decimation. Actually, one definition of time stretching / pitch scaling would be to "modify the signal so that the spectrogram is stretched".

Time stretching and pitch scaling are equivalent from the point of view that the number of samples needed to represent the processed signal changes. Time stretching to half the original length requires only half the number of samples, and pitch scaling to half the frequency allows to discard every second sample without causing aliasing, assuming that the empty upper part of the frequency spectrum is not filled by something made-up by the algorithm used.

Time stretching / pitch scaling does not have a unique definition or a unique mathematical formula. Well-done time stretching should match one's expectation of, in your case, "speeding up" the process that generated the audio. Speeding up natural processes is not uniquely defined. For example, should the pianist be made to move their fingers faster, resulting in changes in note attacks in the least, or should we be more agnostic of the actual process taking place and just isolate the notes in the audio and move them around into a more compact pattern? And should the note decay times be made shorter, too, because we want to make everything shorter? Depends on one's expectations or needs.

There are various algorithms though. Here are before and after spectrograms of the phrase "hello" (freesound), processed by rubberband -t0.5 speech16kHz.wav out.wav with Rubber Band v1.8.2 command line and analyzed in Adobe Audition 3.0 with 512 bands and Blackman–Harris window:

Figure 1. Top: original speech spectrogram, bottom: spectrogram of the speech time-stretched to half the original length and half the number of samples. The spectrograms are zoomed out to full view for easier comparison.

The spectrogram of the time-stretched signal looks essentially the same, but has reduced detail and resolution.


What are the advantages and disadvantages of the Ridgelet transform?

Because of their ability to better represent signals, Ridgelet transforms are preferred over Wavelet transforms in general. However, to be used in practice it has to be discretized. Now, the discretization of Ridglet transform is a challenging task on its own as it involves interpolation in polar coordinates which makes perfect reconstructions (inversions) difficult. For example, Minh N. Do and Martin Vetterli propose an orthonormal version of Ridgelet Transform to achieve better invertability [1].

Ridgelet transform enhances the idea of point-to-point mapping of singularities to point-to-line mappings, which is more effective in handling directions. Yet, even Ridgelets might be ineffective for certain applications such as image processing, where edges are curved rather than being straight lines. In such cases, better transforms such as Curvelets are already proposed and proven to be more powerful.

In this regard, the method of choice really depends on your application. As wiki quotes:

  Continuous Wavelet Transform (CWT) is very efficient in determining the damping ratio of oscillating signals (e.g. identification of damping in dynamical systems). CWT is also very resistant to the noise in the signal.

If your problem is suitable to CWT, then you will be perfectly fine with using it. The Academia involves a lot of works in comparing and demonstrating the strengths of all of these transforms. Some of these would be [2,3,4].

Under the shade of the comments below, I now elaborate more on the theory.

First of all, singularity is simply the concept of point where the function or its derivative is not defined, the point of cusp, double point or any other simlar undefinedness, which I won't go into details here. You can check the wiki "Mathematical Singularity". A good measure of singularity strength would be Lipschitz exponent, which can well be characterized by Wavelet transforms. For more information on that take a look at [5].

But, as you ask, why are you interested in the singularities anyway? Well, because such irregularities basically play an important role in carrying the information that characterizes the signal [6].

Having explained the singularities a bit, I certainly agree that singular lines are a concern of multidimensional structures, more specifically $R^2$ ($L^2$ etc.). This is why transforms like Ridgelets and Curvelets are good at image denoising than Wavelets. Yet, this doesn't violate the fact that 2D Wavelets would also be useful in image analysis.

Regarding the computational complexity, even before analyzing the equations, you could rely on your knowledge that computation of the Ridgelet transform involves a Wavelet transform in the Radon domain. That is why it should definitely be heavier in terms of computational complexity, compared to a Wavelet transform.

.  REFERENCES

http://www.ifp.illinois.edu/~minhdo/publications/FRIT_trans.pdf
http://jstarck.free.fr/curvencyclop09.pdf
http://authors.library.caltech.edu/1381/1/STAieeetip02.pdf
http://www.iro.umontreal.ca/~mignotte/IFT6150/ComplementCours/CurveletTransform.pdf
http://www.iis.sinica.edu.tw/papers/whwang/5114-F.pdf
http://www.cmap.polytechnique.fr/~mallat/papiers/MallatSingDet92.pdf


Autocorrelation function between received and training signal used for channel estimation

I think it's just a trick which makes $R(t)=0, \mathrm{when}\quad t<0$.

First, $-N+1$ is just a time offset and not affect the essence of function.

$r_{CE}(t)$ is bounded to $[0, 2N+T_{CH}-2]$, and $a_N(t)$ to $[0, N-1]$. In the form you gave, $R(t)$ is bounded to $[-N+1, 2N+T_{CH}-2]$. May be the author want the non-zero value of $R(t)$ also starts at 0, so he just move it right with $N-1$ length, which causes the $-N+1$ term in the equation.

--------- Why $R(t)$ is bounded as that -------------

This is just a explanation to give some intuition, not a formal proof. Consider two function $f(t)$ and $g(t)$, which $f(t)=0, t \notin [a, b]$ and $g(t)=0, t\notin [c, d]$. Now we want to calculate the sum of products, i.e., $S=\sum_t f(t)g(t)$, $[a,b]$ and $[c,d]$ must be overlapped, right? Thus we get $a \le d $ and $c\le b$.

Get back to the $R(t)=\sum _d r_{CE}(t+d)a^*(d)$(note here the argument is $d$), $a^*(d)$ is bounded to $[0, N-1]$ and $r_{CE}(t+d)$ to $[-t, L-1-t]$, in which $L=2N+T_{CH}-1$ is the length of $r_{CE}(d)$. Here we get $a=0, b=N-1, c=-t, d=L-1-t$. Using $a \le d $ and $c \le b$, we get $-(N-1) \le t \le L-1$.

As to why $r_{CE}$ is bounded as that, you can make a same process to understand it.


Looking for a way to monitor raw microphone digital output?

So my question is is there a way that the digital output of the microphone can be directly monitored and recorded in units of Volts. Hence, I would use these values to calculate Pascals. Can anyone point me to a (preferably free) program that I can use to do this direct monitoring?

So, in fact, your sound card is a Analog-to-Digital Converter, which does exactly that: convert the voltage coming from a microphone into a stream of numbers. Congratulations, you already have your data acquisition output :)

That stream of numbers has no physical unit (it is "best-case proportional to input Voltage", but it is a number, not a physical entity anymore)– you will need to calibrate it with a couple of known signals.

To make matters more complicated, the amplitude response of a soundcard can't, technically, be flat. Due to the [Shannon-Nyquist theorem], the microphone voltage signal must first go through a low-pass filter that limits the maximum frequency contained in the signal (you're a physicist: apply Fourier transform to time signal $s(t)$, giving you $S(f)$, find $f_{max}$ for which $S(f>f_{max}) \approx 0$).

Since perfect analog filters don't exist, the suppression "edge" can't be infinitely steep, so you'll see changing effect of voltage on digital numbers with changing frequency.

If you work with only one (or few) frequency(ies), and you have a multiple reference pressure wave source for this(those) at different powers, then go ahead and use the sound card, calibrate using the known powers, and calculate correction factors.

If you need something that actually promises linearity, sound cards with their analog filters optimized for audio perception might not be the tool of choice. Verify your approach with an oscilloscope, then use a suitable data acquisition unit – some digital oscilloscope can record signals for you, there's even digital oscilloscopes that have no display and are only meant for usage with a PC, other people I know use Software Defined Radio (SDR) devices to sample audio/sonar frequency signals instead of RF, yet others have just built their own microcontroller-based data acquisition devices based on 10€ evaluation boards.


Consider converting the image into L*a*b* color space — colors with subtly-different component values in RGB color space will often have significantly-different a*b* component values in L*a*b* color space.


dealing with Huge Length (Order) FIR filters

There are three different implementations of FIR filters that I am aware of.

The first is the basic "transversal" FIR in which you sum up the terms of the summation directly. It gets very expensive if $L$, the FIR length, gets long.

The second is the so-called "fast convolution" in which the FIR impulse response is zero-padded to a longer length and the DFT is applied to it once. Then there are two slightly different procedures, "overlap-add" (OLA) or "overlap-save" (OLS, sometimes called "overlap-scrap" which I think is more descriptive) to process blocks of the input into blocks of the output.

Both of these methods are well covered in textbooks.

A third method, not really covered in textbooks at all, is the "Truncated IIR" filters (TIIR) and this might apply in your case.  The simplest example of a TIIR filter is the moving-sum or moving-average that you might find in cascade-integrate-comb (CIC) filters. They have recursion in them, but they are not, mathematically, IIR. Their IIR impulse response gets cancelled or truncated.

If you're trying to implement a long FIR that is a sliding rectangular window (that's a moving sum) or a sliding Hann window or sliding Hamming window, you can do that very efficiently with a TIIR.  A sliding Blackman is a little harder and I don't think TIIR will work for a sliding Kaiser.

I have to go now, but I will return with some references. Or send me an email (it's easy to find) and I will send you copies of the references I have including a "how-to" paper I wrote about how to do biquad TIIR filters in general.


How might I go about detecting phase changes in a sine wave of a certain frequency?

The easiest way to do this would be to apply a digital FM discriminator to the signal. At a high level, you would do the following:

Form the analytic signal of your input. This also can be termed as ensuring it is at complex baseband.
Calculate the phase versus time of your analytic signal $x[n]$:
$$
\phi[n] = \text{atan2}\left(\frac{\text{Im}\{x[n]\}}{\text{Re}\{x[n]\}}\right)
$$
The instantaneous frequency of the signal versus time can be estimated by differentiating $\phi[n]$ using any discrete differentiation method (first-order difference is usually used, i.e. $f[n] = \phi[n]-\phi[n-1]$).

Why does this work? If your signal is typically a constant frequency, your FM discriminator's output $f[n]$ should be approximately constant. If there is a phase jump at a particular sample index, this would appear to the discriminator to be a sudden change to a much different frequency. After the jump, the discriminator output $f[n]$ will again track the sinusoid's frequency.

Therefore, you can run a discriminator continuously on your input signal and apply a threshold, looking for large changes in frequency from sample to sample. That would correspond to a phase discontinuity.

Edit: The above is subject to spurious outputs when the phase value output by the arctangent function wraps around $\pm \pi$ (see the comment below). A more robust way of implementing the above would instead be:

$$
\hat{x}[n] = x[n] x^*[n-1] \\
f[n] = \text{atan2}\left(\frac{\text{Im}\{
\hat{x}[n]\}}{\text{Re}\{
\hat{x}[n]\}}\right)
$$

That is, form the product of $x[n]$ with the conjugate of $x[n-1]$. At each sample instant $n$, the result will have a phase that is equal to the difference in phase between $x[n]$ and $x[n-1]$ (which is what you're trying to measure). The subtraction alleviates any wrapping issues when the absolute phase of $x[n]$ is near $\pm \pi$. Just find the phase of the resulting product to get $f[n]$.


Z transform - Inverse System function - Why number of poles and zeros myst be equal?

Regardless of causality and stability, if you count poles and zeros at the origin and at infinity, the total number of poles always equals the total number of zeros. I'll show a few examples to make this obvious.

First, take a polynomial

$$P(z)=(z-z_0)(z-z_1)\ldots (z-z_M)\tag{1}$$

It has $M+1$ zeros and no (finite) poles. However, it must have $M+1$ poles at $z=\infty$ (because of the term $z^{M+1}$).

Next consider an "all-pole" function

$$A(z)=\frac{1}{(z-z_0)(z-z_1)\ldots (z-z_N)}\tag{2}$$

with $N+1$ poles and no (finite) zeros. Clearly, due to the $N+1$ poles, the function $A(z)$ in $(2)$ has $N+1$ zeros at infinity, because of the term $z^{N+1}$ in the denominator.

In general, if you have

$$H(z)=\frac{(z-z_{0,0})(z-z_{0,1})\ldots (z-z_{0,M})}{(z-z_{\infty,0})(z-z_{\infty,1})\ldots (z-z_{\infty,N})}\tag{3}$$

and if $M>N$, you get $M-N$ poles at infinity. For $M<N$, you have $N-M$ zeros at infinity. Consequently, in all cases the number of poles equals the number of zeros if poles and zeros at infinity are included.


you sample the interval and build a filter for each frequency shift, unless you don’t need to, like in typical carrier/phase recovery schemes.

You make a bank of filters, particularly if Doppler is of itself, of interest.  In SONAR, tones have better Doppler resolution, and the DFT provide a bank of filters.

The ambiguity function function for the waveform will show how much loss you can tradeoff for the frequency sampling you choose.

Each waveform will have different loss for frequency miss-match.  The ambiguity function also accounts for delay resolution.

A waveform like hyperbolic FM is robust to Doppler so a single waveform covers a wide range of possible shifts.

Both miss-matches can couple.

In comms, Doppler compensation is usually a matter of compensation of the frequency shift, not full time dilation.  Shifting a matched filter is a matter of frequency shifting and can be accommodated in carrier/ phase recovery.


strange artifact in depthmap of a mirroring glass

Structured or not, light gets reflected and the sensor has no way of telling if the scene is imaged directly or by reflection.

To try and reduce the effect of reflections from the glass, you might want to use a polariser. Photographic polarising filters have lots of useful applications and one of them is reducing glare when you are trying to shoot artifacts within a museum glass case (for example).

You can find really cheap polarising filters in the second-hand boxes in camera shops, or buy polarising sheets which you can cut into various shapes.

One thing that would be good to note however is that the camera may already be using a polarising element internally which might interfere with an additional polariser which will worsen the situation.

Hope this helps.


Resources for the Basics on Signals and Electromagnetic Waves in the radio frequencies

A classic reference is "The Science of Radio", by Nahin; it's one of my favorite books on any topic. It does have a little bit of math -- there's just no way around it, but it is very well explained and makes few assumptions on your math background.


How to choose suitable segmentation process for image processing

The thresholding techniques that you would find here should be ideally suited to OCR applications. You could read and easily implement. They are also very efficient.


If you're working in the MATLAB environment, the STFT of a 1D signal can be obtained by using the Spectrogram() command. Refer here

And as for PSD, you can use the FFT to obtain the same. Refer here


Question about correlation of two signals

You are correct in noting that cross-correlation is not commutative. If the question says "correlation between two signals $x$ and $g$" then it is usually assumed that you go in that sequence and compute $\Phi_{xg}$.

Or you can choose to say that the question and ambiguous compute both $\Phi_{xg}$ and $\Phi_{gx}.$


The Instantaneous Frequencies from IMF Aren't Consistent with Its Spectrum

This is a major bug in this library. see http://github.com/jaidevd/pyhht/issues/43#issuecomment-398388106 I think it's better to use more robust libraries that are well-tested.


Multipath fading in wireless channel - parameter needed

You don't need to compute the parameter $|\alpha|$, because it is the antenna gain, which is part of the antenna's specification. This is also mentioned on p.21 of the book you quoted (above Eq. $(6)$).


Why is beginning of this FIR filtered signal attenuated?

The filter you designed is a linear phase filter, therefore, response of the filter will start at $N/2$th sample, where $N $ is the filter order. As $N$ increases so does the delay of the response. You can find this delay value by computing group delay of the filter.

In case you want to reduce the delay, you can design a minimum phase filter.


Detecting objects in real world video feed?

There is no straight-forward answer to this. But I can tell you the things you need to consider to choose a method.
Before getting into the methods to detect the object in the reference image, let me tell you about invariance.

If the object in the video is of the same scale as reference image? ( I am guessing NO), then you have to consider detecting at various scale.

If the object in the video is of the same orientation as reference image? (I am guessing NO), then you have to choose a feature that is rotation invariant)

Now about various method that you can use. Remember this is non exhaustive and there are many other methods.

If you have only one reference image, then it is difficult to say without looking at the reference image and the kind of video you are referring to. But a simple answer would be template matching. Other method could be to compute feature descriptors like SIFT/SURF or even color of reference image and then try to search for it in the video. Or you can also use image processing techniques like thresholding, canny edge detection etc.

If you have 100s of images of same reference image, you can train a SVM classifier. The choice of feature depends on the object in the reference image.

If you have lots of different images(like in 1000s) of the reference image class, then try deep learning approach. If you have few few images, it is still possible with this dlib library for deep learning using MMOD(check this), but you need check the accuracy.

So it all depends on the problem in hand specifically.


how does the Frequency spectrum of a signal change in the logarithmic space?

Plus is the simplest operation : $z = x+y$. Fourier is inherently linear, and good at addressing it. However, most processes and data combination are nonlinear, and they should be dealt with. The second simplest operation is multiplication. Homomorphic filtering deals with $z = x\times y$. It is more complicated, especially because of zeroes.

Logarithms were invented to linearize products. Because $\log z =  \log x + \log y$, and then we get back to linear. The problem is that logarithms are not defined everywhere in a simple way, especially with negative numbers. So, it is necessary to deal with that, and I do not know of a naturally sound method. But classically, one offsets and scales data, such as with a  modified logarithm:

$$ \mathrm{l}_m (x) = \mathrm{sign}(x) \frac{\log(1+a|x|)}{\log(1+a)} $$

related to companding. When $a$ and $|x|$ are small, this is close to $$\mathrm{sign}(x) \frac{\log(1+a|x|)}{\log(1+a)}\approx \mathrm{sign}(x) \frac{a|x|}{a}\approx x$$ thus back to the linear case.

You can also check: Why do we substract a background image and not divide it?, talking about LIP (Logarithm Image Processing).


For some context, let's go back to the Kalman Filter equations:

$\mathbf{x}(k+1) = \mathbf{F}(k) \mathbf{x}(k) + \mathbf{G}(k) \mathbf{u}(k) + \mathbf{w}(k) \\
 \mathbf{z}(k) = \mathbf{H}(k) \mathbf{x}(k) + \mathbf{v}(k)$.

In short, for a plain vanilla KF:

$\mathbf{F}(k)$ must be fully defined. This comes straight from the differential equations of the system. If not, you have a dual estimation problem (i.e. estimate both the state and the system model).  If you don't have differential equations of the system, then a KF isn't for you!

$\mathbf{x}(k)$ is, by definition, unknowable. After all, if you knew it, it wouldn't be an estimation problem!

The control vector $\mathbf{u}(k)$ must be fully defined. Without additional system modelling, the only uncertainty on the control vector may be AWGN, which may be incorporated into the process noise.   Known matrix $\mathbf{G}(k)$ relates the control input to the states - for example, how aileron movement affects the roll of an aircraft. This is mathematically modelled as part of the KF development.

The system process noise $\mathbf{w}(k)$ is also, by definition, unknowable (since it is random noise!). However, the statistics of the noise must be known which, for a plain vanilla KF, must be zero mean AWGN with known covariance $\mathbf{Q}(k)$.  Sometimes, the covariance of the noise may change between samples, but in many cases it is fixed and therefore $\mathbf{Q}$ is a constant.  In some instances, this will be known, but in many instances, this will be "tuned" during system development.

Observations are a similar story. The matrix that relates your measurements to the states $\mathbf{H}(k)$ must be fully defined. Your measurements $\mathbf{z}(k)$ are also known because that's the reading from your sensors!

The sensor measurements, however, are corrupted by AWGN  $\mathbf{v}(k)$, which, being random noise, is by definition unknown. The statistics of the noise must be known, which is zero mean with covariance $\mathbf{R}(k)$. Once again, the covariance may change with time, but for many applications, it is a fixed value.   Often, your sensors will have known noise characteristics from the datasheet.  Otherwise, it's not too hard to determine the mean and variance of your sensors that you need to use.  Yes, this can also be "tuned" empirically.

There are a huge number of "tricks" that can be done to work around the restrictions in a plain vanilla KF, but these are far beyond the scope of this question.

Afterthought:

Whilst googling for "Kalman Filter" results in a million hits, there are a couple of things that I think are worth looking at.  The wikipedia page is a too cluttered to learn from effectively :(

On AVR Freaks, there is an "equation free" intro to the Kalman Filter that I wrote some time ago to try to introduce where it is used for real.

If you're not afraid of maths, there are several books worth reading that are at the senior undergraduate/early postgraduate level. Try either Brown and Hwang which includes all the theory and plenty of example systems. The other that comes highly recommended but I have not read is Gelb, which has the distinct advantage of being cheap!


Intersymbol Interference -- Root raised cosine filter redundant if doing impulse response equalization?

There are three reasons to avoid ISI through pulse shaping rather than correcting it via channel correction methods.

As Bruce pointed out, it is a simpler solution and requires fewer computations.
Channel correction usually amplifies the received noise.
Received noise will corrupt the channel model, making the ISI removal imperfect to some degree.


Minimize the Cost Function of Values of Vectors Based on Their Amplitude

Since there is no prior at the Vector level this is basically element wise problem.
Moreover, if we assume the noise to be White Noise with zero mean then the answer can be very simple.
Since the phase difference is always a multiply of 180 [Deg] we can, without loss of generality, assume they are on the real axis.
So what we have can be modeled as:
$$ {z}_{1} = {w}_{1} + {n}_{1}, \; {z}_{2} = {w}_{2} + {n}_{2}, \; \left| {w}_{1} \right| = \left| {w}_{2} \right| $$
So the answer will be $ \hat{w}_{1} = \operatorname{sign} \left( {z}_{1} \right) \frac{\left| {z}_{1} \right| + \left| {z}_{2} \right|}{2} $ and $ \hat{w}_{2} = \operatorname{sign} \left( {z}_{2} \right) \frac{\left| {z}_{1} \right| + \left| {z}_{2} \right|}{2} $.
This comes from a simple intuition. Let's assume bot are positive and $ {w}_{1} = {w}_{2} $. Then of course the answer will be, for any white noise, to average the results. So we are doing the same, averaging them, just taking care of the case they have opposite phases (In the complex plane they are on the same circle just could be on opposite directions).
If you assume AWGN and build the Maximum Likelihood Estimator you'll get the same result.


Zero-Padding seems to CHANGE FFT of signal!

Your problem is these two lines:

t=(1:L)*(1/Fs);
f=(1:L)*df;

They should read:

t=(0:L-1)*(1/Fs);
f=(0:L-1)*df;

If I modify your script as below (for scilab, not matlab), I get the following plot.  The peaks all line up correctly.

And the peak of the FFT is at 0.25 not 0.3125.

for ZP=0:5       // of times of original length of signal to ZP with

    x=[0 1 0 1 0 1 0 1];
    x=[x,zeros(1,length(x)*2^ZP-length(x))];
    X=fft(x);

    Fs=1;
    L=length(x);
    df=.5*Fs/L;
    t=(0:L-1)*(1/Fs);
    f=(0:L-1)*df;

    figure;
    subplot(211), plot(x);
    subplot(212), plot(f,abs(X)); //xlim([-.02,.52])

    test_pnt=.3125;
    I=find(f==test_pnt);
    title('abs[X(' + string(test_pnt) +']='+ string(abs(X(I))));

end;


number of multiplications per second required by the system

One potential option would be to split the lowpass filter into the FIR part and the pureley recursive IIR part.
$$H(z) = H_1(z) \cdot H_2(z) = \frac{1}{1 - a1z^{-1} - a2z^{-2} - a3z^{-3} - a4z^{-4}} \cdot \left [ b0 + b1z^{-1} + b2z^{-2} + b3z^{-3} + + b4z^{-4} \right ]$$
You still have to apply the recursive part to each sample in the upsampled domain, but you only need run the FIR when you need an actual output sample. So this could be
$$N = 5000*4 + (5000*5)/7$$
That's not a great idea in terms of stability and noise performance, but would work if you have enough precision and head room.
UPDATE:
There are further ways you could potentially reduce multiplactions. The FIR section of a butterworth can be split into two second order sections which are simply $[1 2 1]$ so you could potentiall write this as
$$H(z) = H_3(z) \cdot H_4(z) \cdot H_4(z) = \frac{b_0}{1 - a1z^{-1} - a2z^{-2} - a3z^{-3} - a4z^{-4}} \cdot \left [ 1 + z^{-1} + z^{-1} + z^{-2}  \right ] \cdot \left [ 1 + z^{-1} + z^{-1} + z^{-2}  \right ]$$
So in this case it's down to 25000 multiplies per second, although it's not a particularly useful implementation. At the end there are multiple different ways of implementing them and which one is better depends on the properties of the platform.


To show that a system

$$y[n] = T\{x[n]\} $$

is linear, you have to show that it's additive and homogeneous which is codified in

$$T\{ a ~x_1[n] + b ~ x_2[n] \} = a ~T\{x_1[n]\} + b ~T\{x_2[n] \} $$ for all complex $a,b$ and all $x_1,x_2$.

The property set above is satisfied for a system given by :
$$y[n] = T\{x[n]\} = \sum_{k=-\infty}^{\infty} h[k] x[n-k] = h[n] \star x[n] $$ where $\star$ is the convolution operator and $h[n]$ is some sequence to be defined (incidentally it will be the impulse response of the given LTI system).

Of course you didn't have to prove it, as it's already known that if the input-output relation of a system is given by convolution sum; i.e. $y[n] = h[n] \star x[n]$, then it's by definition an LTI (linear and time-invariant) system. Indeed input output relation of any LTI system is given by the convolution sum...


Unexpected result after calculating the response of a system in frequency domain

That's not strange – that's expected!

The DFT is a circular operation – by shifting to the right, your signal "wrapped" around.

What you need to do is called fast convolution based on the DFT/FFT. The Wikipedia has a nice article on the Overlap-Add and the Overlap-Save methods, so I'd recommend starting with these, and then going for the literature that suits your understanding best.