Time-shifting operation post the time-reversal operation when performing convolution

Just think of the value of $t$ where $x(0)$ appears. For $x(t+t_0)$ it is at $t=-t_0$, which corresponds to a left shift if $t_0>0$. However, for $x(-t+t_0)$ the value $x(0)$ occurs at $t=t_0$, which is to the right of its original position if $t_0>0$.
You can think of deriving $x(-t+t_0)$ from $x(t)$ in two different ways:

invert the time axis: $x(-t)$
replace $t$ by $t-t_0$, i.e., shift to the right if $t_0>0$: $x(-(t-t_0))=x(-t+t_0)$

Or

shift $x(t)$ to the left (for $t_0>0$): $x(t+t_0)$
invert the time axis: $x(-t+t_0)$: this flips the shifted function from the left to the right of the original function (again, if $t_0>0$).


Retrieve DFT coefficients from DCT coefficients

Behind an FFT or a DCT operator, which can be implemented as a matrix product. Depending on the shape of your discrete vectors, you may have  $b=D\times a$ and  $c=F\times a$. So you can find your answer with a suitable matrix product. You can also find fast algorithms by decomposing each of the above matrices in simpler matrix product, or using lifting-like implementation. The very fastest can be more difficult to find.


Using the given identities, find the inverse DTFT

As pointed out in the comments, it is important for you to know and understand that multiplication in one domain (time or frequency) corresponds to convolution in the other domain. And, as a consequence, multiplication with $e^{-j\omega}$ in the Fourier domain corresponds to a delay of one sample in the time domain (i.e., convolution with $\delta[n-1]$). So you have

$$\begin{align}\mathcal{F}^{-1}\left\{B\frac{e^{-j\omega}}{1-ae^{-j\omega}}\right\}[n]&=B\cdot\mathcal{F}^{-1}\left\{\frac{e^{-j\omega}}{1-ae^{-j\omega}}\right\}[n]\quad (\textrm{due to linearity of }\mathcal{F}^{-1})\\&=B\cdot\mathcal{F}^{-1}\left\{\frac{1}{1-ae^{-j\omega}}\right\}[n-1]\\&=B\cdot a^{n-1}u[n-1]\end{align}$$

where I used $\mathcal{F}^{-1}$ to denote the inverse discrete-time Fourier transform.


Let me put in an unconnected form things about audio codecs.

Audio can be music or speech (discluding ultrasound, sonar etc).
Music is wideband and requires high-fidelity. Speech has a lower bandwidth and does not require fidelity but intelligibility.

A Codec can be lossy or lossless, the choice depends on the source type, the purpose of application and the resources available.

Codecs are also divided into two as waveform coders (such as AC-3, AAC, MP3, OGG vorbis, (lossless) FLAC, ADPCM) and source coders (such as vocoder, LPC, CELP etc).

A waveform coder, compresses (transforms) the data trying to preserve the waveform shape. On the other hand, a source coder tries to compute parameters of a model that's assumed to produce the given data.

Waveform coders are preferred for high quality multimedia audio, broadcast audio and music compression. Source coders are/can be preferred for very low bitrate speech compression for communication systems such as GSM, military radio etc.

As technology advances, bandwidth limitations have disappeared for certain (but not all) applications. In addition, the compression efficiency of lossy waveform coders improved to the point that medium quality speech compression comparable to that of source coders can be achieved.

So, technically speaking,  your VoIP (voice over internet protocol) application can select from a wide range of speech codecs available for different applications needs and available resources. However you must consult to your VoIP protocol standards to see whether it leaves you free to choose your speech codec or dictates a particular type for use. In short, it's more of a standards issue than a technical one.


Determing the frequency of a signal with a sample at a time

A single sample tells almost nothing about the frequency of a signal.  The only hypothesis that is eliminated is a DC signal of a different value from that of the sample.

The more samples one accumulates over a longer period of time, the better one might be able to fit a maximum likelihood estimate of higher degree models of what the spectrum of the signal might be to those samples, assuming non-aliased samples of a high enough density relative to the bandwidth the signal.


How to determine if a pixel is inside a contour?

I am not in the image processing field, but here is an approach from my control loop background that I believe would apply as an elegant solution for determining if a pixel is inside a contour.  This would be to treat the image as a complex plane and then use the Cauchy Argument Principle (from 1831!!!) to determine if the pixel is inside or outside the contour.

The details and background are below, but as an algorithm you would do the following:

Treat your pixel of interest and your contour as a complex number, meaning a pixel at an x, y location is represented as x+ jy (I am an engineer so trained to you j instead of i).  Using complex number subtraction, subtract the pixel from every point on your contour (which is a new vector with magnitude and phase) and measure the phase of this new vector. Continue along the contour while keeping track of the minimum and maximum phase value in the results after the subtraction. Once the contour is complete (you have completely enclosed your pixel) determine the maximum - minimum phase. If it is 2$pi$ the pixel is inside your contour, if not then it is outside.

See figure below for an introduction to the Cauchy Argument prinicple, and then I will try to explain how this becomes the algorithmic solution above. From this you may define other similar algorithms depending on the details of your contour, such as you will notice that the phase in the example given where the contour phase keeps increasing will also keep increasing when the pixel is inside the coutour, but will increase and then decrease when it is outside.

To explain the principle I will need to explain (in simplest terms) the idea of poles and zeros and encirclement on the complex plane.  For your case we could simply identify the one pixel as a "zero" and proceed from there, which makes the explanation even simpler, so I will avoid any details on encircling poles.
(My plot is titled "Nyquist Plot Background" as this is the theory behind the Nyquist Plot that is used in classical control theory).

Given a ratio of polynomials, the poles are the values from a complex plane that would result in the ratio being infinity, and the zeros are the values that would result in the ratio being zero. Typically the polynomials are factors into first order functions to readily reveal the pole and zero locations. For example:

$$H(s) = \frac{(s+2)}{(s+3)(s+5)}$$ would have a zero at s= -2 and two poles at s=-3 and s=-5. These are the values that would make the ratio go to zero and infinity. At this point for a test of a single pixel, we can proceed with the idea of a single zero so will use that going forward.

In terms of an image with two axis x and y, x would correspond to the real numbers in this case, and y would correspond to imaginary numbers (the solution in the end will be very simple and elegant, but we need to use complex numbers to use it). So in the case above, the zero would be on the image at x= -2 and y = 0.  To position a zero at an arbitrary location with x and y values we would use $(s - (x+jy))$ as the zero.

That said, observe how the Cauchy Argument Principle and complex mapping works in the case of a single zero:

In the following figure, we have a mapping of a zero at location x= -2, y=0 defined by the function $$H(s) = s+2$$

The mapping of H(s) for the example contour shown is on the right. What the mapping is specifically is to take each point on the contour and subtract it from the zero. So for each point $(x_n, y_n)$ on the contour, and the pole at $(x_z, y_z)$ we would get the result $(x_n+j y_n) - (x_z+j y_z)$ which is a complex vector with the magnitude and phase just as shown by the blue vector on the left in the plot below. Note that the magnitude is simply the geometric length between these two points in 2D space, and the angle is the angle found by moving the zero to the origin; by subtracting the two, we are simply moving the zero to the origin which is the complete mapping in the case of the zero (simple!). The Cauchy Principle reduces to saying if you encircle the zero that is inside the contour, the mapping will encircle the origin in the same direction.

In contrast, if the zero is outside the contour, as in the plot below with the zero at x= -4, y =0, the mapping will not enclose the origin.

So using this, we can position your one pixel to be the zero as a complex number at the x,y position, and complete the mapping by subtracting each x,y position on your contour to each x,y position on your pixel. Now to algorithmically test that the result encircles the origin, note the phase for each of the two cases above. For the case of encirclement, the phase will go from 0 to $2\pi$ always. For the case when we don't encircle the origin the phase can never extend completely from 0 to $2\pi$.


Effect of Gaussian Blur on Different Frequency Components of an Image

Gaussian Blur is Spatially Invariant Linear Filter.
Hence it can be analyzed in the Frequency Domain which in fact shows its Low Pass properties.
Namely it attenuates High Frequency Energy.

In Image, Edges, which are abrupt change from one color to other, requires high frequencies in order to be local.
Gaussian Blur attenuates and smear it.
As written above since it is Spatially Invariant, it behaves the same on DC zone or Edge.

Now, there are families of "Edge Aware Filters" which are adaptive to the neighborhood they operate at:

Bilateral Filter.
Non Local Means Filter.
Median Filter.
Guided Filter.

All of those react to the properties of the pixels they operate at hence can not be defined in Linear Convolution (Spatially Invariant).


Why does the specgram of this simple sine wave look strange?

You are computing the spectrogram of the DFT as if it were a time signal. Try this instead

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

t = np.linspace(0,100,1000)
omega = 2*np.pi*1
y = np.sin(omega*t)

spacing = t[1] - t[0]
plt.specgram(y,NFFT=100,noverlap=0,Fs=1./spacing,window=mlab.window_none)
plt.show()


Is there a method to combine speech signals recorded by different microphones?

The easiest approach would just be averaging – since noise on these three microphones should at least partially be uncorrelated, that would improve SNR.

However, that requires that the three microphones are

synchronously sampled (i.e. the $n4th sample from each microphone happened "exactly" at the same time) and
in the same distance from the desired source of sound.

Since both 1. and 2. are usually not fulfilled, you'll have to look into adapting to that – the result is a beamformer. This might be a term often linked to regularly spaced microphone arrays, but that's not a mathematical requirement – you'd just start by detecting the phase/time offset between the strongest (or, better: the desired) source of sound on each of your microphone channels, and then align (read: fractional delay) the channels in a manner that maximizes SNR for that. This inherently leads to a geometric beam shape, though it usually is but a singular value decomposition of a receive signal matrix (that's where "beamforming" and "MIMO processing" meet).


Centering zero frequency for Discrete Fourier Transform

Oh!  What a cool trick!  It works because of the convolution theorem (i.e., multiplication in the spatial/time domain is equivalent to convolution in the frequency domain.)

It's not flipping across the $x$ and $y$ axis, it is rotating the image of the Fourier transform (think of shifting halfway around a cylinder).  The trick here is that alternating -1,1 in the spatial domain is a signal with the highest frequency.  So the Fourier transform of that image is a single point in the frequency domain.  Convolving by a single point is equivalent to shifting (rotating) the image by the offset of the point from the zero frequency.

Here is a test image:
.  It's Fourier transform looks like:

If you take the Fourier transform of the alternating image (), it results in a single point right at the center of the Fourier transform: .  (Recall we haven't done our rotation yet, so the center of the fourier transform is the high frequencies and the low frequencies are at the corners still.)  But this is the "rotation kernel!"  Convolving with this rotation kernel moves everything down and to the right (but things that fall off the bottom right rotate into the top left.)

Convolving the original image with the rotation kernel (in the image domain) gives you: , while convolving the fourier transform image with the rotation kernel (in the frequency domain) gives you: .

And we can check that multiplying the testimage by the checkerboard in the image domain gives , which has a fourier transform of: .


How to calculate the variance of the noise samples $n[j]$ in terms of $N_0$ and $B$, where $n[j]$=$n_f(jT_s)$ and $T_s$ is the sampling period?

Do you know how to calculate the variance of the process
$\{n_f(t) \colon -\infty < t < \infty\}$?
No? Hint: it is the area under the power spectral density curve of $\{n_f(t)\colon -\infty < t < \infty\}$ provided that you are not using radian frequency but are using Hertzian frequency to express the PSD.  Well, if you know the process variance, then be assured that each random variable $n_f(t)$ in the process has the same variance as the process variance. A more interesting question for you to ponder might be: "Are $n[k]$ and $n[\ell]$ independent random variables or correlated random variables?" which will bring in that nuisance parameter $T_s$ into play.


Setting up RRC filter and useful equations

First of all, the correct equation for an RRC pulse is given here.
You are correct that you need to define the sampling frequency Fs and the symbol interval Ts. The number of samples per symbol interval is then Ts*Fs, which I assume to be an integer.
You also need to define beta, and the pulse duration D. For simplicity, let's assume D is given as a multiple of the symbol interval Ts.
Then, you can easily define the time vector. Assuming symmetry around zero:
t = (-D/2)*Ts:1/Fs:(D/2)*Ts;

Then, if you have a function rrc(t) that evaluates the RRC formula for a given time vector t, you're done (all you need to do is transcribe the formula I gave above into Matlab) or whatever language you use).
Note that you don't need to worry about upsampling. The symbols that go into the RRC filter have to be upsampled to space them out by Ts, but the filter itself doesn't have to be.
For more details, see Chapter 11 in the book linked at the bottom of the page here. That book is actually quite practical and may be a good complement to your current textbook.


How to do smoothing without affecting phase

Any smoothing in the frequency domain is similar to convolution filtering.  And convolution filtering (with anything other an impulse or zero) will spread the width of a signal.

If you want a real result, you may need to check to make sure the input to your IFFT is conjugate symmetric.

If you don't want any phase shift, you could try copying the phase from before your filter into the result after your filtering and before the IFFT.


partial fractions expansion inverse Z-transform, help

First convert your input-output relation to Z-domain to see:

$$ Y(z) = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) X(z)$$

From which you find the transfer function of the LTI system as:

$$ H(z) = \frac{Y(z)}{X(z)} = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) $$

Now this is an FIR system with 4 taps. As MattL. tried (without luck) to convince you that for such a system a partial fraction expansion is not the proper way to compute the output. Indeed for any FIR system, the transfer function is of the form of a $H(z) = \sum h[k]z^{-k}$ and for such a system the output is defined through a convolution as the following:

$$ Y(z) = H(z) X(z) =  \left( \sum h[k]z^{-k} \right) X(z) \implies y[n] = \sum h[k]x[n-k] $$

so for your example the proper way to describe the output for the input $x[n]= u[n]$ is

$$ \boxed{ y[n] = 0.25 ~( u[n] + u[n-1] + u[n-2] + u[n-3]) } $$

Nevertheless, assuming you want to insist on a partial fraction expansion to compute the output $y[n]$ then you would do the following:

Given $x[n] = u[n]$ and its Z-transform $X(z) = \frac{1}{1-z^{-1}} $, then the output is:

$$Y(z) = \frac{1}{4} ~\frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} } $$

Now, ignoring the linear scale factor $1/4$ to the end, since the degree of numerator polynomial is larger (in negative) than that of the denominator polynomial, you should first apply a long division of numerator polynomial into denominator polynomial to get :

$$ \frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} }  = -(3 + 2 z^{-1} + z^{-2}) + \frac{ 4} {1 -z^{-1} } $$

hence you have this (including the scale factor):
$$ Y(z) =  -(\frac{3}{4} + \frac{2}{4} z^{-1} +  \frac{1}{4} z^{-2}) + \frac{1} {1 -z^{-1} } $$

And finally you get the output from the PFE method via inverse Z-transformas:
$$ \boxed{ y[n] =  u[n] - \frac{3}{4} \delta[n] - \frac{2}{4} \delta[n-1] -  \frac{1}{4} \delta[n-2] } $$


Finding the Sections of Signals that are the Most Similar

The best way to achieve this, to my knowledge, is a matched filter or cross-correlation. This will tell you not only how well $A$ and $B$ match but will also give you a time offset for it. Just make sure to normalized based on the energy of both signals so you aren't getting false positives if that's a major concern.
If you're looking for a more strict time match (the signal similarity needs to be in-time) then you only want to look at $R[0]$ where $R[n] = A[n] \star B[n]$.
Technically, cross-correlation of non Wide Sense Stationary signals is two-dimensional. So if you're interested in the particular time when the signals best match up with no offset, your process may look something like this over $Z$ trials.
$$
\begin{align}
R[n,m] &= E\{A[n]B^*[m]\}\\
m &= n\\
R[n] &= E\{A[n]B^*[n]\}\\
R[n] &= \frac{1}{Z}\sum_{\zeta=0}^{Z-1}{A_\zeta[n]B_\zeta[n]}
\end{align}
$$
If you clarify what sort of signals these are and what about their matching is important, I think it would help narrow down what solution would best fit.


Need help understanding 16QAM constellation diagram

Yes and no.

The mapping is arbitrary as long as the receiver correctly determines which constellation point a symbol is.  If the receiver makes a mistake, though, it is most likely going to pick a "neighbor" constellation point (i.e. a constellation point that is only one spot away).  It is highly unlikely that a correctly implemented receiver will pick a constellation point far away from the transmitted point.  Communication systems try to mitigate the impact of the "neighbor" mistakes by minimizing the Hamming distance between neighbors.  A map where all of the nearest neighbors have a Hamming distance of 1 is called a Gray map.  If you use a Gray map almost all of the symbol errors will only result in one bit error.


Maximum Likelihood carrier phase estimation of $A\cos(2\pi f_ct)$

for the last line
$$\sin(2 \pi  f_{c}t + \phi_{ML} ) =\sin(2 \pi  f_{c}t)\cos (\phi_{ML}) +\cos(2 \pi  f_{c}t)\sin(\phi_{ML})
$$
so the integral could be written as:

\begin{align}
 \int_ {T_{0}} r(t)(\sin(2 \pi  f_{c}t)\cos (\phi_{ML}) +\cos(2 \pi  f_{c}t)\sin(\phi_{ML}))&= 0\\
\Rightarrow \cos(\phi_{ML})\int_ {T_{0}} r(t)\sin(2 \pi  f_{c}t) + \sin(\phi_{ML})\int_ {T_{0}} r(t) \cos(2 \pi  f_{c}t) &= 0
\end{align}

Therefore,
$$
\tan(\phi_{ML})= -\frac{\int_ {T_{0}} r(t)\sin(2 \pi  f_{c}t)}{\int_ {T_{0}} r(t) \cos(2 \pi  f_{c}t)}
$$

Hence the result.


Unlike a standard high pass filter where you set a cut-off frequency and other design parameters for a fixed filter result with a pass band ripple, stop band rejection, phase response etc.. the "recursive least squares filter" is an adaptive filter commonly used for channel equalization.

As a high level working example of a LMS equalizer, please see my post here where I derived and used the Wiener-Hopf equations to solve for the LMS filter coefficients. This would be an approach to equalization when the coefficients are solved for a block of data and then applied; most applicable to fixed filtering solutions such as compensating for distortions introduced in the RF portion of a receiver (compensating for analog filters for example).

Compensating Loudspeaker frequency response in an audio signal

To adaptively solve for the same coefficients, such as compensating for changing channel effects, two common algorithms are referred to as "LMS: Least-Mean Squared" and "RLS: Recursive Least Squares". I compare the two in my slide copied below, and for more details please refer to the following references

Rappaport, T.S. Wireless Communications Principles and Practice, Second Edition, Prentice Hall 2002

Proakis, J.  “Adaptive Equalization for TDMA Digital Mobile Radio,“ IEEE Transactions on Vehicular Technology, Vol. 40, No. 2 pp 333-341, May 1991

Sayed, Ali, Fundamentals of Adaptive Filtering, Wiley, 2003

Bingham, J.A.C., The Theory and Practice of Modem Design, 1988

The comparison of the implementation for the LMS and RLS adaptive equalizer is shown in the figures below, where vectors are indicated by single bars and matrices with double bars:

LMS: Easy!

RLS: Fast!

Finally to note that a least-squares equalizing filter is NOT the best choice for a channel with severe frequency selective fading as noise enhancement will result where deep fading nulls occur. For this condition, a decision feedback equalizer is often a better choice if time-domain equalization is to be used.


Re-recording audio vs. convoluting old audio data - differences?

Would there be any difference in the resulting audio files using one of the following two aproaches?

There would be a difference - probably negligible - which would represent the non-linearity or time-dependence in the system.


Can discrete wavelet transform for denoising purposes be implemented in real time?

Even though I'm not much of a wavelet expert, I can answer with a definitive YES!!
Now, do you have any messy details to complain about, like the size of the box, the power that it consumes, or the number of frames of delay between the taking of the image and its display?  Because that's your real problem.
Take the particular wavelet transform you want to use.  Sit down and with pencil and paper figure out how many calculations (floating point operations, basically multiplies and adds) that need to be performed per frame.  Look to see if the transform you want to use works with the fast wavelet transform.  Multiply the FLOP per frame by your frame rate -- you should have an answer in MFLOPS or GFLOPS.  Now hold that up against available processors or FPGAs, and figure out what needs to be in your box.  Then multiply by three, because it always takes more real estate to do a calculation than you think.


What is the sequence type if ROC is an annulus in the z-plane

if the ROC of a Z-transform is a ring on the z-plane, then the associated time domain signal is :

two sided infinite length (hence non-causal)

and

stable if the ring includes unit-circle or unstable otherwise.


Where is the Nyquist frequency in an FFT spectrum?

You clearly haven't digested the information from the other answers I gave you.  With an odd number of points, the Nyquist frequency lies cleanly between two bins and does not occur at a bin.
For conceptual purposes, an odd number of points let you set zero at the center flanked by the same number of points on each side.  Likewise, the DFT can be centered at DC flanked by the same number of bins on each side, positive one way, negative the other.
Zero padding either domain means zero padding at the extremes.
More reading, more math, more popcorn:
Using fourier coefficients to reconstruct data in matlab


The problem in your implementation is that it returns the tensor of float, while an image must be a tensor of int. Because of that, your rendering library, which I assume is matplotlib, cannot correctly plot an image.

To fix that you need to specify the type of final_img explicitly. That is, you need to add a parameter dtype=np.int when creating final_img, giving us the following:

final_img = np.zeros((maxx - minx + 1, maxy - miny + 1), dtype=np.int) \
        if len(img.shape) == 2 else np.zeros((maxx - minx + 1, maxy - miny + 1, img.shape[2]), dtype=np.int)

You can also refer to this python implementation of cv2.warpPerspective for more details.


Why is deconvolution used in telescopes, CT scans and seismographs?

As long as the "distortion" is caused by a linear system, its effect can be modeled by a convolution operation with some "filter". The filter is different for different applications and perhaps even different scenarios and is usually measured / estimated before the deconvolution operation.

The deconvolution itself attempts to recover the original signal given the measured output which is the convolution of a signal with a filter (essentially the signal was filtered and we use the result to try and recover the original signal before filtering).


Terminology for instantaneous phase of magnitude component of complex signal

and then discard the phase data

What's gone is gone. You can't reconstruct the phase of the original data unless you have some additional information.

in which the Hilbert transform is applied to this magnitude-only data to estimate the analytic signal, from which instantaneous phase is extracted.

That's a just mathematical process. You shove data in, and data comes out. Whether the result is useful or not depends on what you want to do with it and why you chose this mathematical process in the first place. If you want to associate any physical meaning with the result you need to justify this with some physical property or law.

To clarify, the Hilbert transform applied to the real component of a complex signal is not a good estimate of the original complex signal?

Correct. In general most signals are NOT analytic unless there is a strong reason why they should. If there is a reason to interpret the results that way, it should be stated in the paper. In order to "estimate" a discarded phase, you need have some additional information or (justified) assumptions about your original signal otherwise you might as well just toss the dice. If there is a additional information, it should be stated AND the phase estimation method should be properly derived from it.


Derivation of downsampling in the frequency domain

The expression

$$s_d[m] = \frac{1}{D}\sum_{k=0}^{D-1}\exp(j2\pi \frac{km}{D})$$

is the Fourier series expansion of the periodic function $s_d[m]$ in the discrete case (i.e. only up to a certain amount of frequencies).

$s_d[m]$ is a also called comb function, and I have written more about it (including the relation you ask for) in one of my articles.

Furthermore, here's some example Python code that verifies the correctness:

m = np.arange(64)
D = 16
sdm = 1/D*sum(np.exp(2j*np.pi*k*m/D) for k in range(D))
plt.plot(m, sdm.real)

The picture shows the Comb function, where every 16 samples the values becomes one, else zero.


(Not an answer, but can't post code in a comment)

Works fine for me:

from

 __future__ import division
from numpy.random import randn
from numpy.fft import rfft
from scipy.signal import butter, lfilter
from matplotlib.pyplot import loglog

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

fs = 250
lowcut = 1.0
highcut = 50.0
x_Vtcr = randn(10000)
x2_Vtcr = butter_bandpass_filter(x_Vtcr, lowcut, highcut, fs, order=4)
loglog(abs(rfft(x_Vtcr)))
loglog(abs(rfft(x2_Vtcr)))

What do you get that's different from this?


The time lag is equal to the index difference times the sample period.  The sample period is:
$$
T_s = \frac{1}{f_s} = \frac{1}{16000} = 62.5\mu s
$$
Thus, the time lag from the first peak to the second is $(12 - 2) * T_s = 625\mu s$.  The time lag from the second peak to the third is $(25 - 12) * T_s = 812.5\mu s$.


Question about linear combination of non-stationary signals

No, try with $\alpha_j = 0$, $\forall j$.

I suspect this to be a tricky question, to deceive fast thinking. For instance, having the "same length" does not mean a lot, for a continuous variable $t$. Misleading fake clues to make the question look serious.

Without more hypotheses on the notion of stationarity (in law, wide-sense) and the statistical interrelations between the so-called signals (it would be wiser to talk about processes). See Marcus' answer for more examples.


You could try detecting line segments in the image for example using hough transform, or maybe lswms, edlines or lsd.

After detecting the lines, you can compute their length and angle and find the longest almost horizontal line. If you are lucky, this line will be the horizon and you can rotate the whole image in the opposite direction to get the aligned image.

This might work but it will depend on how good the horizon is visible and how well it is detected.

Maybe you could find a number of potential horizon lines and have the user select the best one.

On second thought, if you involve users into this, just let them draw the horizon and then rotate the image. That's not really what you were asking for, but it will probably work better than what I suggested above. :)


Your pole-zero plot depicts a linear phase, real valued, three coefficient FIR filter, if we assume the zeros are complex conjugates on the unit circle.

The frequency response magnitude is exactly zero at the angle of the zeros on the unit circle. Because of these exact zeros, the log magnitude of the frequency response tends to minus infinity at those frequencies, which gives an impression of notch, however the frequency response is not narrow enough for a typical notch filter. Adding two poles close to those zeros would make it a strong notch filter.

Because of the location of zeros, the frequency response magnitude is larger towards frequencies $\pm \pi$ and smaller at low frequencies and for this reason it is classified as a highpass filter. However it's not a very successful highpass filter either.


How can I measure noise and sharpness of an image?

Simply, you can not 'calculate' the noise from single observation. Rather, you can estimate the target signal (denoised form), and then by using this estimation, you can determine the noise level. If the observation $y = x + n$ where the $x$ is target and $n$ is noise, after estimation of $x$ (via standart denoising methods), you can calculate the MSE, the SNR or other measures for noise level by using estimation $x$ and observation $y$. Of course these values will depend on the denoising method you use to estimate $x$, which can critically change the result.


For most in signal processing, a digital signal (source: wikipedia)

  is a representation of a physical signal that is a sampled and
  quantized

So it represents  data as a finite sequence of  discrete values. Often, digital data is embedded into a specific format to record/store values (for instance 16-bit words), possibly with a specific coding of symbols (think about sign conventions, or floating point standards), to end up in a concrete bit-wise or symbol-wise representation.

A discrete signal is a more abstract thing. It is an abbreviation for a discrete-time signal. It may possess a countable number of samples, and their values are not restricted to be discrete. Most often, the time interval between samples is constant (or uniform). To me, lacunary signals that are uniformly sampled, but for which some samples are missing,  are still considered as discrete.


Fast algorithm for the MATLAB's $\tt impz$ function?

the impz function of matlab, simply computes N samples of the impulse response of a given LTI system from its LCCDE description coefficients a nd b.
As far as I know its efficiency is a secondary issue, since computing an impulse response in real time is not a frequently encountered task.

One equivalent code is simply this:

  % given a,b vectors and number of samples as N
  x = [1 zeros(1,N-1)];
  h = filter(b,a,x);

I guess you know how the filter function is implemented and what its output means. In the simplest terms, the filter() matlab function computes the output vector y for an input vector x of an LTI system: $ \sum_{k=0}^{k=N}a_k y[n-k] = \sum_{k=0}^{k=M} b_k x[n-k]$ , by the recursive approach:

$$ y[n] = - \sum_{k=0}^{k=N} \frac{a_k}{a_0} y[n-k] + \sum_{k=0}^{k=M} \frac{b_k}{a_0} x[n-k] $$

Where the initial conditions of the system is assumed to be zero, inline with impulse response definition and LTI system with initial rest conditions.

Note that in the impulse response computation, the signal $x[n]$ has only a single nonzero entry at $n=0$ and the rest of the computation is a recursion on the signal $y[n]$ alone.


How to plot magnitude and phase response of 2 cascaded filters in Matlab?

Using MATLAB/Octave as the tool, the following approach lets you plot the magnitude & phase samples of the DTFT of the cascade of the two discrete-time LTI filters using their LCCDE coefficient vectors $b[k]$, and $a[k]$, assuming they have LCCDE representations..

Since the cascade LTI system is described as:
$$h[n] = h_1[n] \star h_2[n]$$ and consequently
$$ H(z) = H_1(z) H_2(z) $$

Where the $h[n]$ is the impulse response of the cascade system and $H(z)$ is its transfer function (Z-transform of impulse response). Those $h_1[n]$ and $h_2[n]$ refer to the impulse responses of the individual systems that make up the cascade.

Assuming that individual systems do have rational Z-transforms, then we have the following expression for the composite (cascade) system Z-transfom:
$$
\begin{align}
H(z)  &=\frac{ \sum_{k=0}^{M} B[k]z^{-k} } {\sum_{k=0}^{N} A[k]z^{-k}} = H_1(z) H_2(z) \\ \\
&=
\left(  \frac{ \sum_{k=0}^{M_1} b_1[k]z^{-k} } {\sum_{k=0}^{N_1} a_1[k]z^{-k}} \right)
\left( \frac{ \sum_{k=0}^{M_2} b_2[k]z^{-k} } {\sum_{k=0}^{N_2} a_1[k]z^{-k}} \right)
\end{align}
$$

Where individual coefficients $b_1[k],a_1[k],b_2[k],a_2[k]$ are those that represent the cascaded systems.

We want $B[k]$ and $A[k]$ that represent the cascade system, from which we can use the following Matlab/Octave function to plot the DTFT magnitude and phase responses, assumimg that the cascade is a stable LTI system so that it will have a frequency response function:

figure,freqz(B,A);

Now it can be shown by polynomial manipulations that the coefficients $B[k]$ and$A[k]$ are related to the individual coefficients $b_1[k],a_1[k],b_2[k],a_2[k]$ as the following:
$$ B[k] = b_1[k] \star b_2[k]$$ and
$$ A[k] = a_1[k] \star a_2[k]$$

Hence the required frequency response plot can be obtained in the combined call as follows:

figure,freqz(conv(b1,b2),conv(a1,a2));

where b1,a1 and b2,a2 are the coefficient vectors that represent the individual systems. The method can be generalised to N systems in cascade.


PSS convertion from Frequency to time domain

The PSS signal is generated based on Zadoff-Chu sequence and is mapped to the central 62 sub-carriers apart from DC sub-carrier as explained here. If you want to view the PSS alone, leaving aside other subcarriers in the PSS containing OFDM symbols then yes, you'll have to do a 64 IFFT on the length 63 Zadoff-Chu + DC symbols. It is worth nothing that on the receiver side, synchronization is achieved by locking on to the central 63 subcarriers using a 64 point FFT.


A few points:

Filtering in the frequency domain is complicated and needs a real algorithm like overlap add to avoid time domain aliasing. It's much easier to filter directly in the time domain: [b,a] = butter(1,100/(SampleRate/2); y_filtered = filter(b,a,y); is much better
You probably want to put notch filter on the drone frequency
For broad-band stationary noise a good method is Wiener filtering or spectral subtraction. Many papers have been published on this.


Underdetermined deconvolution of windowed output

Using more conventional notation, let $x_k$ and $y_k$ denote the $k$-th input and $k$-th output, respectively.

$$\begin{array}{rl} y_1 &= 3 = \frac 12 x_1 + \frac 12 x_0\\ y_2 &= 4 = \frac 12 x_2 + \frac 12 x_1\\ y_3 &= 4 = \frac 12 x_3 + \frac 12 x_2\\ y_4 &= 6 = \frac 12 x_4 + \frac 12 x_3\\ y_5 &= 7 = \frac 12 x_5 + \frac 12 x_4\end{array}$$

We have an underdetermined system of $5$ linear equations in $6$ unknowns. Let $x_0$ be a parameter.

$$\begin{bmatrix} \frac{1}{2} & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0\\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0\\ 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2}\end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix} = \begin{bmatrix} 3\\ 4\\ 4\\ 6\\ 7\end{bmatrix} - \frac 12 x_0 \begin{bmatrix} 1\\ 0\\ 0\\ 0\\ 0\end{bmatrix}$$

Using SymPy, we can solve the linear system

>>> from sympy import *
>>> from fractions import Fraction
>>> A = Fraction(1,2) * Matrix( [[1,0,0,0,0],
                                 [1,1,0,0,0],
                                 [0,1,1,0,0],
                                 [0,0,1,1,0],
                                 [0,0,0,1,1]] )
>>> y = Matrix([3,4,4,6,7])
>>> x0 = Symbol('x0')
>>> A**-1 * (y - Fraction(1,2)*x0*Matrix([1,0,0,0,0]))
Matrix([
[-x0 + 6],
[ x0 + 2],
[-x0 + 6],
[ x0 + 6],
[-x0 + 8]])

Thus, the solution set is a line parametrized as follows

$$\mathrm x \in \left\{ \begin{bmatrix} 6\\ 2\\ 6\\ 6\\ 8\end{bmatrix} + x_0 \begin{bmatrix} -1\\ \,\,\,\, 1\\ -1\\ \,\,\,\, 1\\ -1\end{bmatrix} : x_0 \in \mathbb R \right\}$$

Note that if we choose $x_0 = 1$ then we recover the original input vector. You chose $x_0 = 3$ instead.


The spectrum of the sum of two periodic signals

the Fourier Series for a single periodic function, having period $P$,

$$ x(t) = x(t+P) \qquad \forall \ -\infty < t < +\infty $$

is

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \, e^{j(2\pi k/P)t} $$

each harmonic $c_k$ is labelled with its harmonic number, $k$, in the above. the actual frequency of each harmonic, in Hz or units of 1/time, is $\frac{k}{P}$.  the fundamental frequency is the reciprocal $\frac{1}{P}$ of the period.  in angular frequency (radians per unit time), multiply both by $2\pi$.  but, in the harmonic number $k$ there lacks enough specific information to tell you what the actual frequency is of the harmonic.  e.g. the second harmonic of a periodic tone of 220 Hz is the same frequency as the fourth harmonic of a periodic tone of 110 Hz.

so the spectrum of $x(t)$ is about what actual frequencies the energy of a signal is.  it doesn't matter whether $x(t)$ is periodic or not (actually it does with the esoteric mathematics, but we won't go there).  and we use the Fourier Transform, not the Fourier Series to visualize it.

So each harmonic in the Fourier Series above is attached to a dirac delta function, $\delta(f)$.

$$\begin{align} X(f) &\triangleq \mathscr{F} \Big\{ x(t) \Big\} \\
 &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2\pi f t} \, \mathrm{d}t \\
 &= \mathscr{F} \Big\{ \sum\limits_{k=-\infty}^{\infty} c_k \, e^{j(2\pi k/P)t} \Big\} \\
 &= \sum\limits_{k=-\infty}^{\infty} c_k \,\mathscr{F} \Big\{  e^{j(2\pi k/P)t} \Big\} \\
 &= \sum\limits_{k=-\infty}^{\infty} c_k \, \delta(f-k\tfrac{1}{P}) \\
\end{align}$$

If you express the spectra of your periodic signals in this form, the spectra sum to a single spectrum.


Question about symbol detection in ISI channels (use of APP decoder in Matlab)

As you mentioned, applying linear equalization (e.g. filtering with the inverse of the channel) will not yield the optimal performance. Though, what you are searching is not the comm.APPDecoder (which is solely for convolutional codes, i.e. bits with binary $\Gamma$), but you are looking for maximum likelihood sequence estimation (MLSE).

Here, MATLAB offers the comm.MLSEEqualizer object, which you can e.g. use like this:

Constellation = [1+1j, 1-1j, -1+1j, -1-1j];
Channel = [1 0.7, 0.3, 0.2]';
sigma2 = 0.3;

tx = Constellation(randi(4, [100,1])).';

noise = sqrt(sigma2/2)*complex(randn(size(tx)), randn(size(tx)));
rx = filter(Channel, 1, tx) + noise;

mlse = comm.MLSEEqualizer('Channel',Channel,'Constellation',Constellation);

estimated = step(mlse, rx);

[estimated tx abs(estimated-tx)]

program output (only a few rows of it):

   1.0000 - 1.0000i   1.0000 - 1.0000i        0
  -1.0000 + 1.0000i  -1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
  -1.0000 - 1.0000i  -1.0000 - 1.0000i        0
  -1.0000 - 1.0000i  -1.0000 - 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0

as you see, the MLSE equalizer completely recovers the transmitted sequence.


Initial conditions for the LTI systems described as a difference equations

The problem is that non-zero initial conditions cause a term in the output signal that does not depend on the input signal. This explains why a system with non-zero initial conditions can neither be linear nor time-invariant. A linear system must have a zero output for zero input. With non-zero initial conditions the output will generally be non-zero, even for a zero input signal. Alternatively, think of scaling a given input signal. A linear system will have a response that is scaled in the same way. However, the part of the output signal caused by non-zero initial conditions will not scale accordingly, because it's independent of the input signal.

The same is true concerning time-invariance. For a time-invariant system, a shifted version of the input signal must result in an output signal with the same shift. However, the output term caused by non-zero initial conditions will not shift accordingly, as it is independent of the input signal.

Consequently, a system described by a linear difference equation with constant coefficients plus initial conditions is only linear and time-invariant if the initial conditions are zero.


How to do an FFT of a line profile obtained by imagej?

It appears there is no straightforward way to do it from the UI. After looking through ImageJ docs and finding out how to start the interpreter, I've come up with something that is essentially this:

macro "FFT1D Tool - C099A8877099L0088" {
    profile = getProfile();
    fft=Array.fourier(profile,"None");
    s=newArray(lengthOf(fft));
    for (i=1;i<lengthOf(fft);i++) {
        s[i]=2*PI/i;
    }
    s[0]=0
    Plot.create("Fourier spectrum of the line profile","wavenumber","FFT(line)",s,fft);
    Plot.show();
}

If placed in StartupMacros.fiji.ijm in the fiji folder, on the next start of fiji, the FFT1D tool will show up in a toolbar menu and, after a line was drawn over the image, a click on the toolbar button will execute the macro.

A caveat: DFT works great for periodic functions and input line sizes that are powers of two. Manually drawing over the image is likely going to result in a non-power-of-two, aperiodic array. Running this through a DFT will result in Gibbs effect which will ruin the accuracy of the transform.    Window functions (substitute "None" for one of the available window functions) can help ameliorate this effect. It helps if the intensity of the image decays towards the endpoints of the profile line.

Also, check this example of using the Array.fourier function.


How does MATLAB's $\tt upfirdn$ select extra points when downsampling?

The function performs upsampling, filtering and downsampling. The reason for the extra samples is due to the fir filter delay. If you want to reproduce this behaviour for downsampling only you should do the following steps:

Take your input signal, xin and add length(h)-1 zeros at the end. This is done because of the filter's delay.
Use filter(h,1,xin_padded); to filter the signal.
Downsample the filtered signal.


Ratio of the power in the sidebands to the total power in conventional AM

this is John BG jgb2012@sky.com

In AM modulation, as shown in

https://en.wikipedia.org/wiki/Amplitude_modulation#Modulation_index

the modulation index is the ratio between the amplitude of the modulated and the amplitude of the carrier.

So the sought modulation index $a=\frac{1}{Ac}$

modulating signal power: $\frac{1}{2}$

carrier power: $\frac{Ac^2}{2}$

side bands to total power ratio:   $\frac{1}{(1+\frac{Ac^2}{2})}$

solving:  $\frac{1}{(1+\frac{Ac^2}{2})} =\frac{1}{3}$ gives $Ac=2$

so the answer is $a=\frac{1}{2}$


What's the difference between the Hough and Radon transforms?

The Hough transform and the Radon transform are indeed very similar to each other and their relation can be loosely defined as the former being a discretized form of the latter.

The Radon transform is a mathematical integral transform, defined for continuous functions on $\mathbb{R}^n$ on hyperplanes in $\mathbb{R}^n$. The Hough transform, on the other hand, is inherently a discrete algorithm that detects lines (extendable to other shapes) in an image by polling and binning (or voting).

I think a reasonable analogy for the difference between the two would be like the difference between

calculating the characteristic function of a random variable as the Fourier transform of its probability density function (PDF) and
generating a random sequence, calculating its empirical PDF by histogram binning and then transforming it appropriately.

However, the Hough transform is a quick algorithm that can be prone to certain artifacts. Radon, being more mathematically sound, is more accurate but slower. You can in fact see the artifacts in your Hough transform example as vertical striations. Here's another quick example in Mathematica:

img = Import["http://i.stack.imgur.com/mODZj.gif"];
radon = Radon[img, Method -> "Radon"];
hough = Radon[img, Method -> "Hough"];
GraphicsRow[{#1, #2, ColorNegate@ImageDifference[#1, #2]} & @@ {radon,hough}]

The last image is really faint, even though I negated it to show the striations in dark color, but it is there. Tilting the monitor will help. You can click all figures for a larger image.

Part of the reason why the similarity between the two is not very well known is because different fields of science & engineering have historically used only one of these two for their needs. For example, in tomography (medical, seismic, etc.), microscopy, etc., Radon transform is perhaps used exclusively. I think the reason for this is that keeping artifacts to a minimum is of utmost importance (an artifact could be a misdiagnosed tumor). On the other hand, in image processing, computer vision, etc., it is the Hough transform that is used  because speed is primary.

You might find this article quite interesting and topical:

  M. van Ginkel, C. L. Luengo Hendriks and L. J. van Vliet, A short introduction to the Radon and Hough transforms and how they relate to each other, Quantitative Imaging Group, Imaging Science & Technology Department, TU Delft

The authors argue that although the two are very closely related (in their original definitions) and equivalent if you write the Hough transform as a continuous transform, the Radon has the advantage of being more intuitive and having a solid mathematical basis.

There is also the generalized Radon transform similar to the generalized Hough transform, which works with parametrized curves instead of lines. Here is a reference that deals with it:

  Toft, P. A., "Using the generalized Radon transform for detection of curves in noisy images", IEEE ICASSP-96, Vol. 4, 2219-2222 (1996)


CLAHE - Redistributing excess pixels among histogram bins

you have to redistributed 13 pixels(# of clipped) to the histogram.

 // clip histogram
  194
  195             if (clipLimit_ > 0)
  196             {
  197                 // how many pixels were clipped
  198                 int clipped = 0;
  199                 for (int i = 0; i < histSize; ++i)
  200                 {
  201                     if (tileHist[i] > clipLimit_)
  202                     {
  203                         clipped += tileHist[i] - clipLimit_;
  204                         tileHist[i] = clipLimit_;
  205                     }
  206                 }
  207
  208                 // redistribute clipped pixels
  209                 int redistBatch = clipped / histSize;
  210                 int residual = clipped - redistBatch * histSize;
  211
  212                 for (int i = 0; i < histSize; ++i)
  213                     tileHist[i] += redistBatch;
  214
  215                 for (int i = 0; i < residual; ++i)
  216                     tileHist[i]++;
  217             }


$$
|H(z)| = \sqrt{\Re^2 H(z) + \Im^2 H(z)} \\
\theta \big(H(z)\big) = \arctan \frac{\Im H(z)}{\Re H(z)} \\
$$

$\Re z$ is the real part of z;

$\Im z$ is imaginary part.

For a single-pole filter $y_n = b_0 x_n - a_1 y_{n-1}$,

$$
H(z) = \frac{b_0}{1 + a_1 z^{-1}}; \\
\Re H(z) = b_0 \frac{1 + a_1 \cos \omega}{1 + 2 a_1 \cos \omega + a_1^2}; \\
\Im H(z) = b_0 \frac{a_1 \sin \omega}{1 + 2 a_1 \cos \omega + a_1^2}.
$$

With $b_0 = 1$ and $a_1 = -\frac{8}{9}$, your filter may be unstable.

For a low-pass filter, try $b_0 = \frac{1}{9}, a_1 = -\frac{8}{9}$.


Instead of using the FilterDesigner Tool consider using the firls() and firpm() commands directly since they will return optimized filters (in the least squares sense or peak error /equiripple sense respectively) very easily with one MATLAB command. Type help firls or help firpm in MATLAB to read further documentation on how to use these and specifically how the parameters define the transition bandwidth desired.
This can equally be done through the FilterDesigner but I see no need for that added complexity.


What does frequency domain denote in case of images?

But what does frequency spectrum means in case of images?

The "mathematical equations" are important, so don't skip them entirely. But the 2d FFT has an intuitive interpretation, too. For illustration, I've calculated the inverse FFT of a few sample images:

As you can see, only one pixel is set in the frequency domain. The result in the image domain (I've only displayed the real part) is a "rotated cosine pattern" (the imaginary part would be the corresponding sine).

If I set a different pixel in the frequency domain (at the left border):

I get a different 2d frequency pattern.

If I set more than one pixel in the frequency domain:

you get the sum of two cosines.

So like a 1d wave, that can be represented as a sum of sines and cosines, any 2d image can be represented (loosely speaking) as a sum of "rotated sines and cosines", as shown above.

  when we take fft of a image in opencv, we get weird picture. What does this image denote?

It denotes the amplitudes and frequencies of the sines/cosines that, when added up, will give you the original image.

  And what is its application?

There are really too many to name them all. Correlation and convolution can be calculated very efficiently using an FFT, but that's more of an optimization, you don't "look" at the FFT result for that. It's used for image compression, because the high frequency components are usually just noise.


I prefer to use a more explicit notation that consists of the number of integer and fractional bits, like in your last question. I'll take them one at a time:

By Q8, I assume you mean unsigned Q0.8 (no integer bits, 8 fractional bits). In this format, there is no representation for the numbers $1.0$ or $1.5$. The largest value that an 8-bit unsigned integer can hold is $255$, which in this representation would correspond to a fixed-point value of $\frac{255}{256} = 0.9961$.
Similarly to the above, in signed Q0.7 format, there is no representation for $1.0$ or $1.5$. The largest 8-bit signed number is $127$, which would represent $\frac{127}{128} = 0.9922$.
You need to be explicit about whether your number is signed or unsigned. If you have a signed Q1.7 format, that would imply that your values are 9 bits long (1 sign bit, 1 integer bit, 7 fractional bits), which is unlikely. However, assuming that you had such a 9-bit signed number with 1 integer bit and 7 fractional bits, you can represent fixed-point numbers in the range $-2$ to $(2 - 2^{-7})$.

Recall that the conversion between the integer value and the corresponding fixed-point representation is merely a division by two raised to the number of fractional bits. A 9-bit signed integer can hold values in the range $-256$ to $255$, which leads to the fixed-point range I described above.


How does BER performance depend on gray or natural binary?

You are confusing symbol errors and bit errors. What you say is true for the symbol error rate; it only depends on the signal-to-noise ratio (SNR). However, depending on the coding, i.e. the mapping of the bits to the symbols, one symbol error can cause one or more bit errors. The most likely symbol errors are the ones where a symbol is confused with one of its closest neighbors. Gray coding makes sure that those most likely symbol errors cause only one bit error. Other types of coding will result in more bit errors per symbol error. This is why Gray coding minimizes the bit error rate for a given symbol error rate (i.e., a given SNR), if one assumes that symbols are only confused with their closest neighbors, which is very likely for moderate to high SNRs.


designing FIR narrowband pass filter using windowing technique

At 44.1 kHz sampling frequency, an FIR (and possibly also linear phase) digital filter with a pass band of 20 Hz would require very large number of taps. If this is not an issue for you then the design is straight forward by the windowing method. You can use Kaiser window design principle to find the minimum length of a Kaiser window that meets your design criteria (namely the maximum pass and stop band ripples, and the transition band width).

The following OCTAVE / MATLAB commands design a very narrow band bandpass FIR filter,specified by its parameters according to your only given cirteria.
The filter length is not unique, you can play with it to find the minimum
that meets your other criteria (if there are any).

Fs = 44100;        % Sampling Frequency
fc1 = 490;         % cutoff frequency 1
fc2 = 510;         % cutoff frequency 2

L = 20000;         % length of the FIR filter
h = fir1(L , 2*[fc1,fc2]/Fs ,'bandpass');    % bandpass FIR filter.

figure,freqz(h,1,L);
title('Frequency response of the FIR filter');

and the result is: (I've zoomed into the pass band frequency region to observe it clearly...)


Why phase information is usually ignored after Fourier transform?

Simply put:

Because in interesting cases, it remains difficult to interpret, and people tends to use other tools
Yes, potentially a lot of non-stationary features, timings, onsets, discontinuities

For a basic use, like measuring the relative energy of different frequency component, the phase is of little use. As waves tend to last over several samples, sliding/overlapping versions can help retrieving timing/localization of features.

However, phase contains a lot on information, especially on images, and perhaps plays an important  role on perception (e.g. The importance of phase information for recognizing natural images). The importance of phase is recognized in many works, see for instance The Importance of Phase in Image Processing, and it is still an active topic. Multiscale phase might detect onsets, contours, while the distribution of phase coefficients could be used for texture retrieval. The phase computed is wrapped  between $-\pi $ and $\pi $, and unwrapping it in 2D or 3D seems to remain challenging. Hence, some resort to complex transformations that are more localized than Fourier bases, like filter banks or "wavelets", to capture phase information at the correct scales.


Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples

It's your observation interval which creates the main problem.

Your reasoning based on the Nyquist sampling theorem is ok; of course with a pure sine wave at the exact Nyquist frequency you will have troubles and therefore it's wise to relax the sampling frequency (slightly) above that of Nyquist rate, such 2.2 Hz instead of a strict 2 Hz... So this is one problem you will practically have.

But your main problem that appears on the FFT plot is about the spectral resolution due to short observation interval. Mainly with two samples (or one second observation of a 1 Hz sine wave) your FFT resolution will be limited to even less than one Hz. Please search the site for PSD, spectral resolution, FFT bin frequency to get a better understanding of spectral analysis of practical windowed data.

In order to see sharp frequency peaks (ideally impulses) in the FFT output of your sine wave, you should increase the spectral resolution, which requires you to increase the observation interval.

I have modified and also expanded your code to see the result of ideal sinc based interpolator on the (near) critically sampled data. Note that I've mereley included a digital simulation of ideal sinc based interpolator (not the simulation of analog interpolator) to see that it would indeed reconstruct the pure sinudosidal from its given samples taken at close to Nyquist rate. Note that for the ideal sinc interpolator to work, the original signal should strictly be bandlimited, or at least sufficiently so, which will have a lot of consequences on the success and efficiency of interpolation.

f  = 1;                % 1 Hz. sine wave...
Fs = 4.2*f;            % sampling frequency Fs = 2.2*f  ; a bit more than the Nyquist rate.
Td = 25;               % duration of observation ultimately determines the spectral resolution.
t  = 0:1/Fs:Td;        % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end);           % get the resulting final duration
L  = length(t);        % number of samples in the sequence
M = 2^nextpow2(10*L);  % DFT / FFT length (for smoother spectral display, not better resolution! )

x = sin(2*pi*f*t);     % sinusoidal signal in [0,Td]
%x = x.*hamming(L)';   % hamming window applied for improved spectral display

% Part-II : Approximate a sinc() interpolator :
% ---------------------------------------------
K = 25;                  % expansion factor
xe = zeros(1,K*L);       % expanded signal
xe(1:K:end) = x;

D = 1024*8;
b = K*fir1(D,1/K);     % ideal lowpass filter for interpolation

y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);

subplot(3,1,1);
plot(t,x);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');

subplot(3,1,3)
plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');

Below is a plot for the result of interpolation from a set of near critical sampling.

And below is the same simulation with a more relaxed sampling, as you can see the interpolator performs much better for this set of improved samples (consequence of better bandlimitedness)


Lossless algorithm for hyperspectral imaging

I have played with SVD for 3 channel image lossy compression in the past, which prompted me to do a quick search and found the following paper which may help you out.
Randomized SVD Methods in Hyperspectral Imaging

Abstract
We present a randomized singular value decomposition (rSVD) method for the purposes of lossless compression, reconstruction, classification, and target detection with hyperspectral (HSI) data. Recent work in low-rank matrix approximations obtained from random projections suggests that these approximations are well suited for randomized dimensionality reduction. Approximation errors for the rSVD are evaluated on HSI, and comparisons are made to deterministic techniques and as well as to other randomized low-rank matrix approximation methods involving compressive principal component analysis. Numerical tests on real HSI data suggest that the method is promising and is particularly effective for HSI data interrogation.


Simulate packet loss in internet video streaming

For me your question is not clear enough. This questions come to my mind:

what exactly do you want to transmit?
Do you want to transmit raw H.264 bitstream (NAL units) or maybe MPEG2-TS, or RTP packets? Most video transmission applications would use either MPEG2-TS or RTP. I suppose that you are not interested in HTTP or other TCP based streaming solutions.
you are talking about the noise in the internet. Are you interested in simulation of some specific network topology?
what network protocols do you want to use and how does the transmission channel looks like?

You might take a look at NS-3 network simulator. With it you can simulate quite everything and it's also quite easy to use if you are familiar with C++ and network programming. But it is probably not the easiest way to solve your problem.

I'm pretty sure that you can also stream the video with FFMPEG but I don't think that FFMPEG allows you to drop packets based on some error model. But you could use NetEm to generate packet loss, delay or duplication of packets on your Linux system.


Let $\sigma^2$ denote the common variance of the random variables comprising the zero-mean stationary Gaussian process $\{X(t)\}$. We want to find
$$R_{Y}(t,s) = E[Y(t)Y(s)] = E\big[(X(t))^2(X(s))^2\big].$$
Now,
given that $X(t)=x$,
\begin{align}
E\big[(X(t))^2(X(s))^2\mid X(t)=x\big] &=
E\big[x^2(X(s))^2\mid X(t)=x\big]\\
&= x^2E\big[(X(s))^2\mid X(t)=x\big]
\end{align}
where we know (since $X(t)$ and $X(s)$ are jointly Gaussian random variables with covariance $R_X(t-s)$ and correlation coefficient $\rho = \sigma^{-2}R_X(t-s)$) that the conditional
distribution of $X(s)$ given that $X(t) = x$ is a Gaussian distribution
with mean $\rho x$ and variance $\sigma^2(1-\rho^2)$. It follows that
\begin{align}
E\big[(X(t))^2(X(s))^2\mid X(t)=x\big] &= x^2E\big[(X(s))^2\mid X(t)=x\big]\\
&= x^2 \left(\sigma^2(1-\rho^2) + \rho^2x^2\right).\end{align}
Hence, the random variable $E\big[(X(t))^2(X(s))^2\mid X(t)]$
equals $(X(t))^2 \left( \sigma^2(1-\rho^2) + \rho^2(X(t))^2\right)$ and the law of iterated expectation gives

\begin{align}
E\big[(X(t))^2(X(s))^2\big]&=E\bigr[E\big[(X(t))^2(X(s))^2\mid X(t)\big]\bigr]\\
&= E\big[(X(t))^2\left(\sigma^2(1-\rho^2) + \rho^2(X(t))^2\right)\big]\\
&= \sigma^2(1-\rho^2)E[(X(t))^2] + \rho^2 E[(X(t))^4]\\
&= \sigma^4\left((1 - \rho^2) + 3\rho^2\right)\\
&= \sigma^4\left(1 + 2\rho^2\right)\\
&= (\sigma^2)^2 + 2 (\rho\sigma^2)^2\\
&= \left(R_X(0)\right)^2 + 2 \left(R_X(t-s)\right)^2\\
&= E[(X(t))^2]E[(X(s))^2] + 2\left(E[X(t)X(s)]\right)^2.
\end{align}
This differs in subtle ways (not just the use of $t$ and $s$ in place of $t_1$ and $t_2$) from what the OP insists is the right answer as stated in several books that he claims to have read, which answer he has faithfully transcribed into his question above.


First, you need to start with in the s domain where a 2nd order transfer function (2 pole low pass filter) which can be written as follows.

H(s) = s^2 + 2zws + w^2

Lets start with omega = 1 and zeta = sqrt(2)/2, then

H(s) = s^2 + 1.414 s + 1

This filter's response and pole locations are shown here.

To peak the response, simply modify zeta. Here zeta = 0.17 and omega = 1.

H(s) = s^2 + 0.34 s  + 1

Then to move the frequncy, change Omega.
Here zeta = 0.17 and omega = 0.39

H(s) = s^2 + 0.1326s + 0.1521

Now all that is left is to determine the z domain coefficients. To do this, use the bilinear transform as shown here.

For the the last example I gave, use  A=1, B=0.136, C=0.1521, D=E=0, F=1, T=2 to calculate the a's and b's.

See this page for a derivation of the equations given here. Here is a rather simple tool that can generate any 2nd order response you want.


MATLAB: How to filter a signal using Thiran Fractional Delay Filter

try this:

h=thiran(0.3,1);
num=  cell2mat(h.Numerator);
denum=cell2mat(h.Denominator);
output=filter(num,denum,x)

Hope that helps


When you segment an image you have boundaries between segments.
Those boundary pixels (Just between 2 different segments) are the edge pixels above.


What is the practical meaning of the variance, covariance, mean value?

Mean of a signal can be practically visualized as the dc average value present in the signal (for a complete sinusoidal period), for e.g

Variance of a signal is the difference between the normalized squared sum of instantaneous values with the mean value. In other words it provides you with the deviation of the signal from its mean value. It gives you the spread of your signal's data set.

Co-variance of a signal consists of unnormalized relation with the signal correlated with itself and also with other signals. It provides you the degree of similarity between the signals of interest.

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

http://www.dspguide.com/ch2/2.htm

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

EDIT : FOR COMPLETE INSIGHT PLEASE REFER ALSO TO @Dan Boschen's ANSWER.


I will try to answer on your question with this: Least Square equalizer is minimizing the intersymbol interference of your signal (ISI), between the symbols, which does not translate directly to minimization of the bit error rate. When I am saying 'directly' that means that there is no exact linear dependency between the error reduction in the equalizer, and the bit error rate. However, they generally tend to go into the same direction.

The reason for that is the nature of equalizer, as the nonlinear estimator.

The chain of processing that you have shown consist of processing modules, and some of them are nonlinear, which interferes with the bit error rate  minimization. In other words, bit error rate minimization is not a criterion that is used by the equalizer.

Below is the mass of papers that are talking about it, and it is not a simple matter.

Please note the first one: "Approximate Minimum Bit-Error Rate Equalization for Binary Signaling" which actually argues exactly the same line of inquiry that you are advocating: "Although most linear and decision-feedback equalizers are designed to minimize a mean-squared error (MSE) performance metric, the equalizer that directly minimizes biterror rate (BER) may signiﬁcantly outperform the minimum MSE equalizer".

This means that the digital bit error based criterion would be a better choice for reducing the bit error rate than the one that is typically used.

http://www.ijesit.com/Volume%202/Issue%203/IJESIT201303_07.pdf

https://www.mathworks.com/help/comm/examples/ber-performance-of-different-equalizers.html?s_tid=gn_loc_drop

http://barry.ece.gatech.edu/pubs/conference/icc97.pdf


How to filter uniform noise in a microcontroller?

You need to check the noise bandwidth using an oscilloscope. By saying noise I mean the hardware noises before digitalization by ADC.
If the noises have high frequency component, remove it by analog filters. High frequency is defined as frequency higher than 1/2 of the sampling rate. And in practice, usually a sampling rate of 8-10 times of the pass band of the analog low pass filter is used. You may try simple RC filters.
Use FIR in a MCU to manipulate your data. Depends on the speed requirement, you can choose number of taps to use and etc. You should be able to find libraries provided by your chip manufacture, since FIR is very common in use.

Hope it helps!


How the noise in time domain affect the power in frequency domain after FFT?

You have problems in defining the noise and sampling frequency. Your noise is not zero mean and your sampling rate is not well defined. Following simple Python code computes the DFT magnitude of the sine plus noise example. (Note: I'm not a python user so please modify the code as necessary)

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
import numpy.random as npr

N = 1000                # Set signal sample length

t1 = -5;                # Simulation begins at t1
t2 =  5;                # Simulation  ends  at t2
x = np.linspace(t1,t2,N);# Define the sampling grid as x[n]

Ts = (t2-t1)/N;         # Compute resulting sampling period
Fs = 1/Ts;              # Compute resulting sampling rate

noise = (npr.random(N)-0.5); # Generate ZERO MEAN uniform noise

fsin = 0.2*Fs;            # Sine Wave Frequency in [0:Fs/2]
y = np.sin(2*np.pi*fsin*x) + noise    # Generate SINE + Noise

# Take FFT and display its MAGNITUDE
w = np.linspace(0,2,N)
waves = fft(y,N)
plt.plot(w,abs(waves))

Where the resulting plot is:


Cutoff frequency of a first order recursive filter

That formula for the cut-off frequency is a very inaccurate approximation. In this answer I derived the exact relation between the coefficient of a first order recursive averaging filter and its 3-dB cut-off frequency. Note that in the quoted answer I used the constant $\alpha=1-b$. From formula $(3)$ in that answer we get for the coefficient $b$

$$b = 2-\cos(\omega_c)-\sqrt{(2-\cos(\omega_c))^2-1)}\tag{1}$$

With $\omega_c=2\pi f_c$ we get for $f_c=0.25$ a value $b=0.26795$.

The frequency response of the filter is given by

$$H(e^{j\omega})=\frac{1-b}{1-be^{-j\omega}}\tag{2}$$

The figure below shows the magnitude response $|H(e^{j\omega})|$, and it shows that the desired cut-off frequency is achieved with the chosen value of $b$:

Note that the 'EDIT' part of the answer quoted above also explains where the inaccurate formula for the cut-off frequency comes from. The following figure shows the resulting cut-off frequency of the filter as a function of the desired cut-off frequency when the exact formula $(1)$ is used and when the approximation is used:

Clearly, the approximation is quite bad and it only kind of works for very small cut-off frequencies. I see no reason why it should be used at all.


How to detach overlapping contours but not to remove small ones. Alternative to erosion-dilation cycle

I don't know if this is a corner case or the norm in your dataset but it is a relatively easy situation to deal with. It would be much more difficult to detect trees in an urban environment, for example, where you would probably be looking at more pattern-recognition type of solutions.

The basic idea behind this solution is surface fitting, whether to directly detect the trees or detect the "floating" level of the ground and through that information detect the trees.

The following is based on Octave and using the 1D data but generalisations to other platforms (and 2 dimensions) should be straightforward.

The first thing one can do here is kernel density estimation (KDE). In this case, you would be approximating the total available surface (both the ground and the trees) via a set of Gaussian curves. In the two dimensional case, this will give you either one Gaussian per tree or higher densities of Gaussians around the locations of trees. In one shot. No further processing required. In reality, you might find that this technique approximates your disparity image via kernel density and returns the data of a set of Gaussians that you can process further to determine where there are trees.

But, having had a quick look at the data, it contains sharp discontinuities and that might not be good for KDE convergence. While you could "blur" the disparity with a Gaussian low pass and then try to approximate it with KDE, this means that your original resolution is decreased and it might reduce your ability to discriminate between small and large trees. Anyway, this is an option too (first blur with a suitable FWHM Gaussian, then approximate with KDE).

Again, this seems a relatively easy case where you could even use simple thresholding to detect the trees. The problem there of course is that the ground height now modulates your threshold. A tree on flat ground could be anything above 25 meters but certainly not a 3 meter tall shed on top of a 22 meter hill.

So, we want to threshold but we want to take into account the ground contour. Ideally, you could fit a curve that is closer to what the ground really is (an ellipsoid) but maybe you can get away with simpler fitting schemes.

The main idea here is that we want to "undo" the modulation of the ground to then apply a simple threshold. If we flatten the ground, the 22 meter hill will disappear and then anything above 22 meters is more likely to really be a tree.

If you ignore the zeros (which I am assuming are "missing values" in your case (either high uncertainty or actually missing values)), the signal looks like this:

You can threshold this at $\approx 62$ but you can see that the second tree is lower than the first but further up the slope of what looks like a hill.

Let's flatten the hill. Flattening the hill here is basically equivalent to subtracting the slowly varying component (the hill) from the signal itself and possibly even controlling the amount of removal with a coefficient. Here, I will use simple polynomial fitting but anything that could approximate the slow varying component of the height field would do.

xrange=0:(length(b)-1);
sv = polyval(polyfit(xrange,b,2),xrange); # Fit and evaluate a second order polynomial. You might want to go higher if you have ridges or depressions but not too high because we just want to approximate the slow varying ground changes.
plot(b);hold on;plot(sv,'r');grid on; # Plot everything

This now looks like this:

Obviously, we can simply subtract the red curve from the blue and get back the "flattened" version of the signal, but as you can see, the trees have "pulled the ground up" slightly and this is why if you do this blindingly, you will get this:

In terms of comparing between the trees, we are now definitely in a better condition because we have removed the ground and we can also tell that the second tree is lower than the first, but, we now have a dip from trying to remove the hill. To adjust for situations like this, you can introduce a coefficient by which to control how much subtraction you want to perform, for example:

m=0.43;plot(m.*b-(1-m).*sv);

Here, $m \in [0 . . 1.0]$ controls how much of the real signal and how much of its approximation we are taking into account. I determined 0.43 after trial and error and this results to:

Where basically, the trees start at almost the same height.

So, now you can threshold anything above -8.2 and be more certain that you are selecting trees.

There is a way by which you can flatten the ground almost completely but it requires some knowledge about the size of the tree canopy.

This works on the same premise but approximates the ground using a different method.

Basically, you can use a median filter to "chop the peaks" that the trees are causing to the height field. Anything that remains after that is basically ground variation.
The median filter is applied in a "rolling window" way and substitutes a sample with the median value of $N$ values surrounding it.

Now, to apply this, you need to determine that $N$ which should be almost as wide as your widest tree canopy. I have determined this to be approximately 125 (pixels) in this example, so by doing:

plot(b-medfilt1(b,125)); grid on;

You get something like this:

Which is almost too good. At least the trees start at exactly the same point. Again, it is now the point of applying a simple threshold (anything above 1 (?)) to detect where the trees are.

This process will give you a binary mask which you can then shrink (or even simply cluster, with simple methods) to determine a centroid that you would then associate with the location of a tree. Obviously, processing the size of each "cluster" could probably also tell you how wide is the tree canopy.

In terms of generalising this, you can either fit a 2D polynomial or a 2D NURB or even simply apply a 2D median with medfilt.

Hope this helps.


What is the advantage of using a Stepped Frequency waveform over FMCW?

If you're defining the stepped-frequency waveform as a pulse train where each pulse has a different frequency, then it indeed does have advantages. One of the main ones being bandwidth requirements, which we'll go over.
Both FMCW and stepped-frequency techniques use frequency-domain range finding via the DFT. They map beat frequencies resulting from the mixing process to a range. For a given LFM chirp of bandwidth $\beta$ and pulse width $\tau$ the mapping is given by
$$f_b = \frac{-2R\beta}{c\tau}$$
The range resolution capability of any waveform with captured bandwidth $\beta$ is
$$\Delta R = \frac{c}{2\beta}$$
In order to get fine resolution, $\beta$ must be large. In a FMCW system, the entire instantaneous bandwidth $\beta$ must be supported for a single chirp. Once you begin exceeding 100 MHz or so, this can become a big problem.
Stepped-frequency waveforms go around this issue by achieving the desired bandwidth by collecting multiple pulses, with each pulse's carrier frequency increasing by some $\Delta F$ (usually linearly). By doing this, the receiver has to only support the bandwidth of one individual pulse of length $\tau$, which is approximately $1/\tau$.
The desired range resolution can be achieved since the bandwidth necessary is now captured over $M$ pulses and is given by
$$\Delta R = \frac{c}{2\beta} = \frac{c}{2M(\Delta F)}$$
As you may already be thinking, there are disadvantages to this. The timing complexity increases and you must spend more time collecting pulses in order to achieve the desired range resolution.
Edit: Follow-up Questions

Does a step-frequency waveform ease the load on transmitter because its easy to generate?

Traditionally yes, it is more difficult to generate high chirp-rate LFM pulses when compared to simple single-frequency modulated pulses since analog devices were used. Back in the day, having oscillators that could provide the required chirp-rate were expensive to design and build. Today were we usually use a DAC, that limitation has been eliminated. You simply have to meet the sampling criteria, handle the sample throughput, some filtering and you're off to the races.

I would like to understand what kind of problems arise in LFM when the receiver has to handle a large BW?

An example of one of the many issues involves the antenna. Antenna's have defined beam patterns for a particular frequency and are designed with this in mind. When using LFM, the frequency changes intra-pulse and so does the wavelength. This means that the antenna pattern is changing as the pulse is being transmitted. For small bandwidths this can be tolerated. Once you start having hundreds of MHz of bandwidth, the antenna pattern can change drastically, which is usually undesirable.

In a stepped-frequency waveform, is the frequency constant within a chirp? If yes, then you would not get the beat frequency by mixing the transmitted and reflected signal?

In the waveform I described, there are no chirps. It is a train of simple modulated pulses, where each pulse has a higher (or lower if downchirping) frequency than the previous pulse. An example looks as follows:

We collect each individual pulse after it is reflected from a target and "stack" them to construct a pseudo-LFM sequence. Due to the mixing process, this will yield a beat frequency just like in FMCW, and we can take the DFT to process range.
There is a trade-off in the $\Delta F$ chosen and the number of pulses. You ideally want  $\Delta F$ to be small so that you approach an ideal LFM pulse, but then that requires more pulses to capture the bandwidth $\beta$.


Is sum of period samples in the DSP domain?

Ladies, Gentlemen,

Please let me say that my answer is wrong on that lower (first) frequency is represented by first least sum. By trial and error I found that lowest frequency is represented by the fourth (absolutely) least sum, indeed given that existing frequencies are two. I apologize. However, I have solved my problem and do not need any other question. My answer to my question is that this process is DSP, but I ask Messrs Administrators view so that if I have any other question about this method, ask it here.

Following C++ program outputs sums of input values (xn) and you can see that 4th least sum represents first (fundamental) existing frequency. You can change values on the condition second existing frequency be double the first. Example values are: sample rate 96000 ksamples/sec, 1st existing freq = 100 Hz, represented by 960 samples, 2nd = 200 Hz, represented by 480 samples, n = 1 to 1970 with step 11. Fourth least sum (138131) represents 1O1.373 Hz. Next frequency, 100.209 Hz, with sum -357858, is the nearest to 100 Hz. Obviously there is accuracy error. However for step n = 1, fourth absolutely least sum is -7.40625 representing 100 Hz. Eighth least sum represents half the 1st frequency, 50 Hz.  I point out that in detecting fundamental frequency of periodic signal by this method, precise is high when ratio sample rate/fundFreq is high, in opposite to method by DFT where results are bad when ratio samplFreq/existFreq is too high.

With regards and friendsip.

// Outputs sum of input values.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
float n = 0; // sample.
float xn = 0; //pure number.
float xn_cos = 0; //pure number.
float xn_sin = 0; //pure number.
float sum = 0;

//Sampling frequency and samples of segment.
int samplFreq = 96000; // sample/second.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 100; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 200; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

for (n = 1; n <= 1970; n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
   + ampl_2 * sin(n * unitAngle_2 + phase_2);
sum += xn;
std::cout << samplFreq/n << "       " << sum << "\n";
n += 11;
}
return 0;
}


Approximating $ {L}_{0} $ Norm Minimization with Non Linear Convex Inequality Constraints using Reweighted $ {L}_{1} $ Minimization

I think testing it practically rather that its theoretical inspection might help more. Therefore, I tried to solve the above mentioned problem using CVX toolbox in Matlab.

cvx_begin;
   variable x(1000);
   minimise(norm(x,1));
   subject to
      x.^2<=150; % f(x) here, is x^2
      0<=x<10;
cvx_end

I think the answer to your question depends on your function f(x), and sometimes even if there is a solution based on convexification of the stated problem, I do not see any guarantees that the answer is must be right. I mean using l1 norm as proxy for l0 norm (as what CS suggests) in this case, might be wrong and bring no fruit. Be careful with blind application of l1 norm minimization method. This great lecture from Prof. Baranuik called "Compressive nonSensing" might be help you more on this regard:
You might find it here:
http://www.norbertwiener.umd.edu/FFT/2015/


Note that if $x[n]$ and $y[n]$ are non-zero in the interval $n\in [0,2N-1]$ (and, for convenience, they are defined as being zero outside that interval), the (linear) convolution sum can be written as

$$y[n]=(x\star h)[n]=\sum_{m=0}^nx[m]h[n-m],\qquad 0\le n<4N-1\tag{1}$$

Note that the upper limit is $n$ because the argument of the term $h[n-m]$ becomes negative for $m>n$.

So if, for whatever reason, you're only interested in $y[n]$ in the interval $n\in [0,N-1]$, you can set the upper limit in the convolution sum $(1)$ to $N-1$ (assuming that $x[n]$ and $h[n]$ are defined as zero outside the interval $n\in [0,2N-1]$). However, the result wouldn't change if you chose the upper summation interval as $2N-1$, because you would just add zeros to the result.

In short, both formulas in your question give the same result, under the given assumptions on the definitions of $x[n]$ and $y[n]$.


Mysterious notation used in "BLIT" CCRMA paper

What a weird typographical error. It's a minus. Get a working copy of the PDF from http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.87.


What I believe is being done is interpolation in the frequency domain to increase the number of pixels while smoothing the transition from black to white.

Consider a single row for example, which is a rectangular function in time. The Fourier Transform would be a Sinc function, which would dominate the outer part of the spectrum in the frequency domain (go toward zero in the middle, which represents $F_s/2$ where $F_s$ is the sampling rate, and then rise back as a mirror of the first half up to $F_s$. Adding zeros in time simply interpolates more samples of this same function in frequency. Then truncating the center of the spectrum reduces the higher frequency range of the Fourier Transform, which effectively makes the pulse wider (lower frequency).

I discuss zero padding for FFT's in more detail here: What happens when N increases in N-point DFT


Assuming $P$ is the transfer function of your process and $C$ the transfer function of your controller.
The closed-loop transfer function of a standard control loop, with the controller in the normal path, is
$$ G(s) = \frac{C(s)P(s)}{1+C(s)P(s)}$$
while the transfer function of your alternative controller is
$$G_{alt}(s) = \frac{P(s)}{1+C_{alt}(s)P(s)}
$$
If you set your new alternative to controller to
$$
C_{alt}(s) = \frac{1+P(s)C(s)-C(s)}{P(s)C(s)}
$$
You will get the same closed-loop transfer function (G(s))
That being said, not sure it's a good idea for 2 reasons
1 - To get to same result as the original controller, you need to cancel poles and zeros in the alternative closed-looped transfer function. If you cancel a RHP zero with a RHP pole, your final transfer function will seem stable however your controller will not be, the output of your alternative controller could become unbounded.
2 - Even if you don't cancel RHP zeros with RHP poles, this method is way more complicated.


Output of a stable LTI system (discrete)

So I'm answering my question here; turned out that I have to solve (2) then solve (1).

To prove $y[n]$ has finite energy, we write

$$\sum_{n=-\infty}^{\infty}y^2[n] = \sum_{n=-\infty}^{\infty}(\sum_{k=-\infty}^{\infty}h[k]x[n-k])^2 = \sum_{n=-\infty}^{\infty}(\sum_{k=-\infty}^{\infty}h[k]x[n-k])(\sum_{l=-\infty}^{\infty}h[l]x[n-l])$$
$$ = \sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}h[k]h[l]\sum_{n=-\infty}^{\infty}x[n-k]x[n-l]$$
Now we know that $x[n]$ has finite energy, so
$$ \sum_{n=-\infty}^{\infty}x[n-k]x[n-l] \leq \sum_{n=-\infty}^{\infty}x^2[n] = E$$
(a generalization of $a^2 + b^2 \geq 2ab$)
Then
$$ \sum_{n=-\infty}^{\infty}y^2[n] \leq E\sum_{k=-\infty}^{\infty}|h[k]|\sum_{l=-\infty}^{\infty}|h[l]| $$
For a BIBO system, $$\sum_{l=-\infty}^{\infty}|h[l]| \leq M \leq \infty$$
So we approve $$\sum_{n=-\infty}^{\infty}y^2[n] \leq M^2E \leq \infty$$

To prove $x[n]$ tends to a constant:

$$y[n] = a\sum_{k=-\infty}^{\infty}h[k]u[n-k] + \sum_{k=-\infty}^{\infty}h[k]x_0[n-k] = a\sum_{k=-\infty}^{\infty}h[k]+y_0[n]$$
Using the conclusion of 1, we know that $$\sum_{n=-\infty}^{\infty}x^2[n] < \infty \rightarrow \sum_{n=-\infty}^{\infty}y^2[n] < \infty$$
which implies that $$\lim_{n \rightarrow \infty} |y_0[n]| = 0.$$
So
$$ \lim_{n \rightarrow \infty} y[n] = a \times \lim_{n \rightarrow \infty} \sum_{k=0}^{\infty}h[k] = ac$$
where $c$ is a constant (because of BIBO).


Why there is no paper optimize the consumption power of the circuit in energy harvest?

Everyone wants to minimize the power of their circuits, so there's nothing special about energy harvesting devices, at all.

Also, every paper about making communications more reliable while is about making the power usage lower, simply because if you want constant reliability, you can then achieve that with lower power.

Basically, all the research on transmitters  uses the energy per bit as the critical measure, so you really don't need to look for energy harvesting-specific papers.

So, this really boils down to basic textbook knowledge of transmitter design, not to speciality papers.

Now, that's on the theoretical side of things. When it comes to building silicon doing the things that you derived are good in theory, then you'll find that modern silicon's manufacturing techniques are very much limited by the power consumption of devices – so, the secrets that you're asking for a worth billions of Qualcomm's, Apple's, ZTE's, and (formerly) intel's money.


Significance of an impulse in the frequency domain

The Dirac delta is not strictly a function but a distribution. The Dirac delta is such that $\delta(x)=0 \ \forall x\neq0$ and it has to meet the following restriction:

$$\int_{-\infty}^{\infty} \delta(x) \ \mathrm{d}x = 1$$

This means that the unit impulse must integrate $1$ over all the real numbers.

Let's define this other "function" $$\tilde{\delta}(x)=2\delta(x)$$

As we can easily see, this new "function" must integrate $2$:

$$\int_{-\infty}^{\infty} \tilde{\delta}(x) \ \mathrm{d}x = 2\int_{-\infty}^{\infty} \delta(x) \ \mathrm{d}x = 2$$

And this is what we get when multiplying an impulse by a constat number. The area under the "function" changes and, due to the Dirac delta being non-zero only at the origin, this means that the impulse must change its height (if we can call it so) in order to change the area it integrates.

In the frequency domain, these constants act as weights that determine something like the relevance of a frequency in a given time-domain signal. If we have a small (i.e. it integrates a small area) impulse in a determined frequency, then there is a pure sinusoid of that frequency present in the signal whose amplitude is rather small. On the other hand, a large impulse would correspond to a pure sinusoid with a large amplitude.


At the very least, you'll need a sampling rate that's high enough to capture the full bandwidth of your Wifi transmission. In a typical WiFi receiver, that puts the sample rate at exactly the OFDM width: 5, 10, 20 or 40 MHz, depending on your Wifi.
Your dropped samples simply happen because your computer isn't fast enough to keep up with that many samples per second. There's only two solutions to that:

Get a faster computer
Make your computer do less per sample.

Also note that the sustained write rates of permanent storage devices can be lower than your sample data rate, so even just writing to disk can be a fatal bottleneck.
Sadly, we don't know the processing you're currently doing, but maybe have a look at gr-ieee802-11, which is optimal enough that a modern midrange laptop can receive 20 MHz wide Wifi with it, continuously.


Inverse Fourier transform of $\frac{\frac{2}{1+2i\pi f}-\frac{2}{4+2i\pi f}}{\frac{1}{1+2i\pi f}+\frac{1}{3+2i\pi f}}$

This looks correct; obviously you need to simplify to obtain a transfer function with polynomials as numerator and denominator. Partial fraction expansion then leads to a form that can be directly transformed to the time domain. The only thing that might be wrong is the minus sign between the two exponentials.


Why symbols have to be upsampling, if the symbol rate and sample rate is same what happen?

The baseband information-bearing signal $s(t)$ is a train of time-shifted orthogonal pulses: $$s(t) = \sum_k a_k p(t-kT),$$ where $a_k$ is the sequence of symbols and $T$ is the symbol rate. The problem is, how to generate this signal in Matlab (or any other discrete-time numerical language)?
To start, note that the $j$-th pulse $a_j p(t-jT)$ can be written as $$a_j \delta(t-jT) \ast p(t).$$ Here, the pulse $p(t)$ is interpreted as the impulse response of a filter (the "pulse-shaping filter"), which allows us to write the pulse as the result of a filtering operation (a convolution). If we define the train of impulses $\Delta(t) = \sum_k a_k\delta(t-kT)$, then the entire pulse train may then be written as the convolution $$s(t) = \Delta(t) \ast p(t).$$
This signal is easy to generate in a computer. Assume that the sampling frequency is fs, the sampling interval is Ts = 1/fs, and the symbol interval is T=LTs for integer L, which is the number of samples per symbol. Assume also that we have a vector of symbols ak, and samples from the prototype pulse p. Then, Δ = upsample(ak, L) is the discrete-time equivalent of $\Delta(t)$, and
s[n] = conv(Δ, p)

is the discrete-time equivalent of $s(t)$.


Independence of Noise at Each DFT Output

Let $m=k-k'$. Then the sum in your final equation becomes

$$\sum_{n=0}^{N-1}e^{-2\pi imn/N}=\frac{1-e^{-2\pi im}}{1-e^{-2\pi im/N}},\quad m\neq 0\tag{1}$$

where I've used the formula for the geometric series. Noting that $e^{-2\pi im}=1$ for any integer $m$ shows that the sum in $(1)$ equals zero.

[Note that your sums should only have $N$ elements, i.e., their indices should go from $0$ to $N-1$ (instead of $N$).]


How do you "split apart" the output of a linear estimator?

There is nothing wrong about subtracting the bias. For the Weiner filter, the operation is a matrix multiply, there are no non-linear operations taking place which is reason everything works out. There are two possible ways to do it:

Subtract the bias before the filter. In this case you'll have a Weiner filter, $G_1$, which will produce $\hat{x}$ given $x+n$.

\begin{align}
\hat{x}&=G_1(y-c)
\end{align}

Subtract the bias after the filter. In this case you'll have a Weiner filter, $G_2$, which will produce $\widehat{x+c}$ given $x+c+n$ (this is the case you presented).

\begin{align}
\widehat{x+c}&=G_2(y+c)\\
\hat{x}&=\widehat{x+c}-c
\end{align}
Note: Be careful with the notation of $G$. The Weiner filter is constructed using the data so the filter constructed using the $x+n$ data will be different from the filter constructed using the $x+c+n$ data.


what is the differnence between 8x8 dct coefficients and 25 x 25 dct coefficients?

There exist different variants of discrete cosine transforms. They mostly vary on a phase shift. The most commonly-used type is known as type-II. Since the works of K. R. Rao et al. in the seventies, it was shown to possess near-optimal properties.

The DCTs are invertible, possess certain symmetry properties that make them useful and practical in sound (signal) or image processing. They can technically be applied to any patch size 13, 8x8, 25x32... Power-of-two patches benefit from faster and more optimized computations. 8x8 are standard for natural image due to:

a size small enough so that image chunks are stationnary enoguh to be quite-well represented by cosine waves
a small memory footprint for digital implementation

Now, it is more common to find 4, 8, 16, and 32 sizes in image and video coding, and even non-square blocks (4x8) with close to DCT transforms.

The quality is not a real question for the DCT, since it is invertible. It is more related to a match between your data, the transform on the processing (e.g. compression). You might have data and applications where a 25x25 DCT would work better. Then you will probably hit some difficulties for optimized implementations.


Reconstructing band-limited signal's integral from discrete samples

Let me rephrase what you want to do as:

  I want to approximately calculate the integral of the product of two non-negative discrete-time signals over time. Each signal is sampled from a band-limited continuous-time signal at a sampling frequency $2N$ times its highest frequency. A method is preferred that gives about 5 % error at most, is computationally efficient and does not require a large oversampling factor $N.$

Oversampling factor

The bandwidth of the product of two signals is equal to the sum of their bandwidths. Therefore an oversampling factor of $N = 2$ will suffice to represent the product without aliasing.

Integration

Integration is a linear time-invariant operation, basically a continuous-time filter, and applying it does not increase the bandwidth of the signal. What we want to achieve is, with the product signal as "input":

$$\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

However, the impulse response of an integrator is of infinite bandwidth and sampling it to form the impulse response of a discrete-time integrator would cause aliasing, corrupting its frequency response. We can insert to the signal path an ideal lowpass filter that has its cutoff at the bandlimit of the input signal. It does nothing to the signal. The output will remain the same as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass}}\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

We can combine the lowpass filter and integrator to a single filtering operation. The impulse response of this lowpass integrator filter is the convolution of the impulse responses of the lowpass filter and the integrator filter. This shuffling of operations is enabled by the algebraic properties of convolution, namely associativity. We still get the same desired output as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass integrate}}\xrightarrow{\text{output}}$$

The lowpass integrator impulse response:

$$\begin{align}
h[k] &= \int_{-\infty}^k\frac{\sin(\pi x)}{\pi x}dx \\
\\
 &= \int_0^k\frac{\sin(\pi x)}{\pi x}dx + \frac{1}{2}, \\
\end{align}$$

is sampled from the integral of the sinc function, without aliasing. Sinc is the impulse response of the ideal continuous-time low-pass filter with a cutoff frequency at $\pi$.

Figure 1. Impulse response $h[k]$ of an integrator. The impulse response continues indefinitely beyond what is shown, approaching values 0 to the left and 1 to the right.

The impulse response almost represents the process of calculating the sum of (optionally the current sample and) the past samples. As this sum becomes larger and larger, the difference to true integration becomes arbitrarily small, proportionally. The values of $h[k]$ approaching 1 for large $k$ means that for very old samples just summing them gives virtually no error in the integral estimate.

Summary

To summarize, oversample the input signals by a factor of two and calculate with a large accumulator a running sum of obtained samples of the product. This will have an arbitrarily small error compared to true integration, as enough data gets accumulated.

Analog filtering

You should also consider what the analog filtering of the current and voltage signals does to the integral of their product. This can be a big source of error. For example let's say that both signals equal $\sin^2(x).$ If this is filtered to remove all oscillation, it becomes a constant signal $\frac{1}{2}$. The average value (power) of the square of the unfiltered signal is $\frac{3}{8}$ and for the filtered signal $\frac{1}{4}$.