So,what exactly is Source entropy?

When the paper talks about "source entropy" all they mean is the entropy of the information source.  You can see that from the following passage in the paper- "Shannon showed that there must exist at least one encoding of the sequence generated by an information source which allows error-free transmission of the sequence when the channel capacity is larger than the source entropy."

  Now, problem_Set also mentions about source entropy but the formula is very different!

Yes, and no.  The formulas appear different, but they are not really.  The first formula is
$H_{source} = lim_{k→∞}\frac{1}{k}H_k$

$H_{source}$ is the total entropy of the source, and $H_k$ is the entropy of each possible discrete value that the random variables in $H$ can be.  Thus, it is almost a trivial statement that the total entropy is the sum of the individual value entropies, while the number of discrete values can go to $\infty$.

The problem set equation is-
$H(s) = \sum_i p_ilog_2\frac{1}{p_i}$

This is really the same equation, it just replaces the generic "$H_k$" with the implied $log_2\frac{1}{p_i}$.  The $\frac{1}{k}$ in the first equation is equivalent to the $p_i$ in the second equation- i.e. they are assuming that all of the discrete values of H are equally likely.  This is not generally true, but is true when the entropy is maximized.


From the audio clip posted on soundcloud, it looks like there's a discontinuity in the output every ~2060th sample. I suspect that you process your audio stream by blocks and that you do not save/restore the state variables of your filter between each call to the filtering function - so every time you process a new chunk of audio it is as if it was preceded by zeros - which is not the case obviously.

Related question.


Your proof is correct, and the result is that the system is time varying because the response to $x(t-t_0)$ does not in general equal a delayed response to the input $x(t)$. Of course, for the special case $t_0=2\pi k$ the delayed response to $x(t)$ equals the response to $x(t-t_0)$, but for time-invariance this equality must hold for any $t_0$.


Fir filtering operation? Also convolution?

Yes in general both IIR and FIR filters can be derived from the convolution summation, but since and IIR by definition has an impulse response that goes to infinity, to implement this directly would require an infinitely long convolution. The FIR can be implemented from the convolution expression directly, but a recursive subset from difference equations is used to implement IIR filters.  So convolution still takes place in both cases (the output is the convolution of the input with the impulse response), except in the case of the IIR filter the generalized convolution expression cannot be used directly for implementation.

Given the generalized convolution summation as:

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

$h[n]$ is the impulse response of the filter and $x[n]$ is our input signal.

Here is a convolution with a Finite Impulse Response (FIR), which is causal (starts at $k = 0$, and has a finite impulse response (summation is over $M$ samples) As written this is directly realizable as an FIR filter:

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

In comparison here is a convolution with an Inifinite Impulse Response (IIR), which is also causal but the impulse response extends to infinity:

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

This is impossible to create a direct realization, but a subclass of IIR systems can be implemented using recursive difference equations:

$$y[n] = x[n] -  \sum_{k=0}^{N} h[k]x[n-k]$$


Help in understanding if the maximum likelhood estimation is working properly

The maximum likelihood estimator will give you precisely the true value for the parameter for a vanishingly small number of instances --- unless the noise is zero.

Your noise instance w is not zero-mean, which will probably give more of a bias away from the true value.

If I take your numbers are just change the last entry in w from -0.65 to 0.1 then I get the "right" answer. Note that this changes the mean of w form -0.16 to -0.042 (i.e. closer to zero).

Your values:

The new values:


Extract a signal that is known to repeat every 4 samples at a known phase i.e. exactly 500 Hz signal in 2000 Hz data

To estimate the amplitude of the 500 Hz signal, a simple moving average filter will give you the best estimate in the presence of white Gaussian noise. This can be easily done by doing a moving average on the result of converting the 500 Hz signal to baseband as described below. This filter will have a Sinc response in frequency, so has high sidelobes which means it would be sensitive to the presence of other stronger signals if present. A very simple work-around is to just window the time domain signal first and then follow the same process (downcovert to baseband and sum  the result over N samples where N is the length of the moving average).   This is essentially a bandpass filter (using the windowing filter design approach) but in any filter used there will be a trade between frequency selectivity and time of convergence.

A convenience of having a "digital IF" that is exactly $Fs/4$ is you can down-convert to baseband (create the equivalent baseband analytic signal) by just multiplying by +/-1 as follows:

To get I you multiply by

I = 1 0 -1 0 1 0 -1 0

To get Q you multiply by

Q = 0 -1 0 1 0 -1 0 1

As these would be the values for the Local Oscillator $e^{-j\omega t}$ for $\omega = \pi/2$ which is the normalized radian frequency at an LO of $f_s/4$.

I detail this approach further at this post:

DAC and ADC architecture in SDRs


The approach you are alluding to involves a) first registering the images and b) using a subtraction operation to remove the background (which is the mannequin).

First on (a) (registration): Image registration deals with figuring out the displacement of corresponding pixels in two images. The problem with real-world images is that the world is 3D and because of the depth in the scene, nearer surfaces will move more than surfaces deeper in the scene. Also, a surface can self-occlude certain parts of the scene. I only state these so that you are aware of the complexities involved. Having said that, in the scenario you have posted above, you might be able to do very well with a simple translation model for the displacement. In other words, you can assume that all pixel in the second image move by $(\delta{x}, \delta{y})$ from its original position. So if a pixel was located at $(x,y)$ in $I_1$, then its location in $I_2$ will be $(x+\delta{x}, y+\delta{y})$. This can be assumed to be true for every pixel (global translation). You now need to estimate the horizontal shift $\delta_x$ and the vertical shift $\delta_y$. There are many approaches to this problem. Here is an academic reference to registration methods (may be TMI) Zitova and Flusser. You may wish to start off by trying the ITK toolkit and try using the "translation transform". If you find that after accounting for a translation transform, the image subtraction still leaves parts of the mannequins, you might want to try some more advanced transforms. For example, the camera may have rotated slightly in your hand. So in addition to a shift, you can also introduce a rotation. Also, you might have taken one shot closer to the other, so there may be a difference in scale. This necessitates the use of an affine transform which is also included in ITK. Barring these "global" methods which assumes that every pixel moved in the same way, you can also use local methods. For example, you can estimate a $(\delta_x, \delta_y)$ for every patch in the image. You can divide up the image in patches of $7 \times 7$ pixels and search for the best corresponding patch in the other image. Every pixel within this $7 \times 7$ window is assumed to have moved in the same manner as the motion that was estimated for the patch. A last resort, but the most general is to compute dense optical flow for every pixel in the image. This gives you a shift for every pixel in the image. You can find an optical flow implementation in the popular OpenCV software.

Finally, on (b) (occlusion): You will inevitably end up with some pixels without a valid shift because they simply aren't present on the other image. Or some pixels are uncovered in the second image which were covered on the first. See these paper for clues Occlusion aware optical flow and TV-L1 Optical Flow with Occlusion Detection


Which antialiasing filter before equispaced sampling

Yes, you have to apply another filter before downsampling in order to avoid aliasing. Your original signal has been acquired with a sampling rate $f_\mathrm{s} = 1.2\,\text{kHz}$. Taking every second sample effectively means that you divide the sampling rate by two, so your new sampling rate $f_\mathrm{s}'= f_\mathrm{s}/2 = 600\, \text{Hz}$. To avoid aliasing you have to make sure that the signal contains no frequencies above $f_\mathrm{s}'/2 = 300\, \text{Hz}$ before you apply the downsampling. If it is certain that your signal does not contain any frequencies above $f_\mathrm{s}'/2$ then you don't need a filter.

The anti-aliasing filter can be implemented by an FIR filter. Its properties like stopband attenuation and width of transition band depend on your signal: if you have spectral components near to 300 Hz, the filter needs to be quite steep and it will consequently have a high order. If hardware implementation is an issue, have a look at Cascaded Integrator-comb filters.


I don't understand harmonics, why do they happen?

Harmonics "happen" when your input to an FFT isn't a pure unmodulated sine wave.

Any unexpected distortion in your input waveform generation (from being exactly identical to mix of sin(wt) + cos(wt)) can be the cause of harmonics appearing in an FFT result (above the noise floor and any windowing artifacts).

Those harmonics are required to represent the energy of any differences whatsoever between a periodic signal and a perfect sinewave of the same frequency.  If harmonics aren't there, then there can't be any differences (assuming a integer periodic input), because a single result bin of an FFT can only represent a pure sinusoid.


The entropy $H(X)$ of a continuous random variable $X$ is infinite. Proof is trivial (note that we can, without loss of generality, use the natural logarithm, since any other logarithm is the same but with a factor: $\log_b(x) = \frac 1{\ln(b)}\ln(x)$):

\begin{array}
\\H(X) &= E\{I(X)\}\\
&= \int\limits_{-\infty}^{\infty}{-\ln\left(P(X=x)\right) dx}
\end{array}

The point here is that for any really continuously distributed random variable, the probability that a specific real value gets "hit" is $P(X=x)\equiv 0 \,\forall x$, and hence, the information contained in that event happening would be infinite. You hence cannot directly compare discrete alphabets with dense intervals (such as the whole $\mathbb R$).

Hence, people invented differential Entropy. For your real-valued (=1-dimensional) case, this is nothing but the obvious:

$$h(x) = \int\limits_{-\infty}^{\infty} f(x)\log_b\left(f(x)\right)dx \text{,}$$ with $f$ being the continuous probability density function of $X$.

Now, all you need to compare your discrete case with the continuous is model the discrete PDF for the alphabet $\mathbb X$

$$
f_\text{discrete}(x) = \sum\limits_{x_i\in\mathbb X} \delta(x - x_i)P_\text{discrete}(x_i) $$

$\delta$ is the dirac function.

Caveat not everything you know about the discrete entropy applies to differential entropy. Read the Wikipedia Article!


Actually the answer is already provided in the wiki page as follows:

"Since the Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing..."

Indeed those 1-D differentiation and averaging kernels are pre-multiplied to form the 2D final kernels to compute directional derivaties.

Hence, even though as you correctly expected the row and column derivates can be computed via 1D kernels, the averaging function being utilized into, renders them to be 2D.


You are right that a distributed system could be "something like a transmission line". Note that the system

$$y(t)=x(t-T)\tag{1}$$

is a simple model of a transmission line, where just a frequency-independent delay $T$ is taken into account, and the attenuation is neglected.

Note that lumped electrical systems, described by resistors, capacitors and inductors, result in ordinary differential equations, i.e., the output $y(t)$ depends on the time derivatives of the input and output signals. If the system is also time-invariant, i.e., the component properties do not change over time, you can define a transfer function, and with lumped elements that transfer function is always rational. The transfer function describing the system given in $(1)$ is

$$H(s)=e^{-sT}\tag{2}$$

which is clearly non-rational.


Introduction to fast algorithms for time–frequency analysis

If you are interested in discrete linear systems (either time-frequency or time-scale), you could invest in multirate filter banks, that provide tools for computing, optimizing, etc, such as the polyphase matrix or the lifting scheme. A  tutorial paper is

Multirate digital filters, filter banks, polyphase networks and Applications: A Tutorial by P. P. Vaidyanathan.

An earlier reference is

Frequency-domain and multirate adaptive filtering by J. Shynk.

However, additional notes could be helpful on you side whether you need:

specificities: amount of memory, real-time constraints, quantization
analysis only or analysis/synthesis
linear or nonlinear time-frequency
discrete or discrete approximation to continuous
fast in number of operations, convergence, or using hardware capabilities (finite-precision, approximate versions, etc.)

If you do consider sliding windows for Fourier analysis, this can be cast into the framework of modulated complex  oversampled filterbanks, for which fast algorithms do exist, see for instance:

Fast Implementation of Oversampled Modulated Filter Banks, 2000, S. Weiss and R.W. Stewart


Denoising approach for a combination of several ADC voltage channles

According to your processing chains (filter --> sum vs sum --> filter) the two approaches, by being LTI systems, should produce the same results under infinite precision arithmetic. You can show this mathematically by using the properties of convolution operation.

Note that a moving average filter is just an ordinary LTI system.

Note also that nonlinearities such as quantization operations may affect the above fact's validity.


Why do we need the power spectral density?

Power Spectral Density (PSD) is a theoretical construct, required to deal with random processes (signals). It's by definition the Fourier transform of the auto-correlation function (another theoretical construct) of the random processes.

Practical computation (estimation) of the PSD can be done in a number of ways and some of them utilize the way you expressed as the (scaled) square of the magnitude of the Fourier transform of an observed instance of the random process. Such an estimator of the PSD is called as the Periodogram (or its various variants).

Therefore, they are not substitudes for eachother. In practice you can only estimate the true PSD.


How to come true the convolution with a 2-D array in frequency domain with python

This is because the identity for the convolution is
[[1, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0]]

You get the correct result by replacing
img_back = np.real(img_back)

by
img_back = np.fft.ifftshift(np.round(np.real(img_back)))

It is important to np.round (to the nearest integer) before casting to np.int, otherwise you will may have wrong results due to rounding errors e.g. np.array([0.999999999999]).astype(np.int) = 0
Notes about fftshift
The the fftshift and ifftshift operations are equivalent to a convolution, we could express fftshift(X) as $S * X$.
Convolutions are commutative and associative $a * (b * c) = (a * b) * c$, and $a * b = b * a$, it is immediate in the frequency from the multiplication operator properties.
Lets say $X$ is your image $H$ is your kernel
you can check that $S^{-1} * X$, or $X = S$, is precisely the term that I mentioned in the beginning of the answer is the identity for the convolution
What you have from my solution above is $S^{-1} * (H * X)$, the solution you metioned in the comments is $H * (S^{-1} * X)$, and there would be one other variation $(S^{-1} * H) * X$ that would represents calling fftshift on the kernel instead of the image. I already saw some doing ffthisft in both inputs and then ifftshift to the output.
Alternative forms...
$(S^{-1} * X) * H = X * (S^{-1} * H) = S * (S^{-1} * X) * (S^{-1} * H)$


Why do I have frequency leakage in DFT after zero padding if frequency resolution is fine?

This phenomenon has nothing to do with spectral leakage. What you are observing is the effect of zero padding. Given a number of samples $N$, there is a maximum possible frequency resolution $\Delta f$ that can be achieved:

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

In your case $\Delta f$ is exactly $2\;\mathrm{Hz}$. If you zero-pad your signal, there is no extra information to retrieve - you will only decrease the frequency spacing.

In the example above, when you increase $N$ to $1000$, you get a frequency spacing of $1\;\mathrm{Hz}$.  All extra observed samples are merely an interpolation, done by the window function ($\mathrm{sinc}$ in your case). You will start observing side-lobes of the window spectrum.  Since you implicitly multiplied your signal by a rectangular window, this will result in the convolution of the spectrum of your signal (two Dirac's + DC) with the $\mathrm{sinc}$ function.

Another way to look at it is to imagine that DFT is basically a filter bank, consisting of shifted $\mathrm{sinc}$ functions. Those are aligned in such a way, that peak of one is where zeroes of all the remaining ones are present. If you start looking in between those zeros, you will start taking those samples. Here is an example plot of such $\mathrm{sinc}$ filter bank.

Let's imagine that frequency corresponding to the blue filter is present. That will yield the amplitude in a corresponding bin. All the remaining frequencies are not present (orange and yellow), thus you multiply those $\mathrm{sinc}$'s by $0$ and get nothing in the bins. In the case of zero padding, that will be no longer the case. Samples of that blue $\mathrm{sinc}$ will fall in the intermediate bins and will be sinc-interpolated.

Here is what happens for $N=1000$ and $N=10000$:

And a zoomed part:

Things to notice:

For $N=500$, there is no leakage whatsoever. There are perfect spikes, representing each of your frequencies and DC offset.
We can also observe the FFT noise at the very bottom.
For $N=10000$ the shape of $\mathrm{sinc}$ function is clearly visible.

And obviously the code to reproduce the results:

Fs=1000;
Ns=500;
Ns2=1000;
Ns3=10000;
t=0:1/Fs:(Ns-1)*1/Fs;
f1=10;
f2=400;
x=5+5*sin(2*pi*f1*t)+2*sin(2*pi*f2*t);

X1 = abs(fft(x))/length(x);
X2 = abs(fft(x, Ns2))/Ns;
X3 = abs(fft(x, Ns3))/Ns;

F1 = 0:Fs/Ns:Fs-Fs/Ns;
F2 = 0:Fs/Ns2:Fs-Fs/Ns2;
F3 = 0:Fs/Ns3:Fs-Fs/Ns3;

plot(F1, 20*log10(X1))
hold on
plot(F2, 20*log10(X2))
plot(F3, 20*log10(X3))
xlim([0, Fs/2])
grid on
legend({'N=500', 'N=1000', 'N=10000'})


An image shouldn't look pixelated after filtered by a Gaussian (lowpass) kernel... instead look at the following figure:


How to append two bandlimited signals and make the result bandlimited without modifying the first signal

That's a tricky problem. First you will need to define how to design "integrity". What do you value more: time domain shape, spectrum, phase, envelope, short term Fourier spectrum envelope, etc. You will probably have to trade between those, so prioritizing would be helpful.

In general this is done bu doing "cross fades", i.e. you overlap the two signals for a small time period and then you fade out F(a) while at the same time fading in F(b). The time and shape of the fade can be controlled to satisfy the limited bandwidth requirement.

However, it seems you can't do this, since you can't modify F(a) at all. One potential work around would be to either extrapolate F(a) or use a  circular repetition of F(a) to create enough data for a cross fade. This depends on the exact properties of F(a)


Kernel Convolution in Frequency Domain - Cyclic Padding

In my StackExchange Signal Processing Q38542 GitHub Repository (See Applying Image Filtering (Circular Convolution) in Frequency Domain) you will be able to see a code which implements 2D Circular Convolution both in Spatial and Frequency Domain.

Pay attention to the function CircularExtension2D().
This function align the axis origin between the image and the kernel before working in the Frequency Domain.
Remember that for Discrete Signals the implicit assumption on signals, In frequency Domain analysis, is being periodic (Circular).
This is basically what's in FFT Based 2D Cyclic Convolution.

In the discrete case one could indeed apply Circular Convolution by element wise multiplication in the Frequency Domain.

With proper padding one could apply linear convolution using circular convolution hence Linear Convolution can also be achieved using multiplication in the Frequency Domain.

See:

Regarding your questions:

The filter is just an array of numbers. As long as you are after 2D Circular Convolution there is no constraints on the Filter. If it is valid for 2D Spatial Circular Convolution it is valid for Frequency Domain Circular Convolution.
Up to Floating Point Quantization errors both are mathematically equivalent (See Convolution Theorem).
If the Convolution Kernel is similar in its size to the image and both are large enough Frequency Domain Convolution becomes more efficient than Spatial Domain.


Hilbert transform Use as a Single Sideband Frequency Translator in Analog communication

There are many applications for the analog equivalent of the Hilbert Transform, often implemented with “90 degree power dividers” also referred to as “quadrature splitters” as well as broadband variants implemented with phase tracking networks.

One application is single sideband frequency translators done with IQ mixers; I have another posting that details this further here:

Frequency shifting of a quadrature mixed signal

Another application is phase shifters and vector modulators where the signal is split in quadrature (To I and Q components) then combined in 0 degree phase after weighting the I and Q components according to the desired phase and amplitude. This can provide frequency translation similar to the approach with IQ mixers. The architectural difference is the use of absorptive bi-phase attenuators in each quadrature leg instead of reflective mixers which is more appropriate for static phase and amplitude control.

Other common applications for the Hilbert Transform that are done digitally can also be done in the analog with the purpose of converting a double sideband IF signal to be single sided (positive or negative frequency only) or creating a baseband analytic signal (which must be complex).

The link I gave above that demonstrates using the Hilbert Transform to frequency shift a signal is converting the signal to single side band frequencies only. Since the OP further asked this question I will explain that further in this context.

First consider Double Sideband frequency translator which would use a single mixer and as shown below is demonstrated as an upconverter translating a sinusoidal tone with frequency $\omega_s$ to a carrier of frequency $\omega_c$. Here we see that both sidebands of the baseband signal with frequency $\omega_s$ have been translated up to be centered on $\omega_c$:

In comparison, the implementation for Single Sideband frequency translator uses two 90° power dividers as shown in the diagram below. I have designated one output for each of the power dividers to represent the imaginary axis by putting a label $j$ on those ports just to show similarity to the Hilbert Transform: there is no reason I cannot define one output to represent such even though many may argue that both are real signals-- whether a signal is real or imaginary rests entirely in our definition and consistency with our rules of math. Regardless the same block diagram would similarly work with real signals, and is generally considered as such. If you swap one of the quadrature legs the resulting output would be the other sideband.


Higher size of your FFT best the resolution of frequencies, but of course it has a cost resulting in a longer to process !

To apply overlap you need define one window size to analyse and one hop/step to walk around your data!

This picture illustrates how this works:

Here a simple code using Window size = 2048pts and 50% overlap step!

n_samples = length(yourdata);
winSize = 2048;

%50% overlap or 0 to not use overlap
OverlapStep = 50;

if OverlapStep > 0

    Overlap = floor((OverlapStep*winSize) / 100);
    nFrames=floor(n_samples/Overlap)-1;
else
    Overlap= winSize;
    nFrames=floor(n_samples/Overlap)-1;
end

k=1;
while ( (k+winSize-1) <= n_samples )

    FrameSignal = yourdata(k:k+winSize-1);

    %%%%%%%%%% Start process FFT here  %%%%%%%%%%%%%%%%%
    w = flattopwin(winSize);
    wdata = FrameSignal(:).*w;
    Y = fft(wdata)
    .
    .
    .

   k=k+Overlap;
end


Isolate recurrent pattern vs. evolving background sound

So, assuming the linear superposition you suggest, the sound you want to exctract lives in a periodic subspace of your signal space. This subspace is uniquely characterised by the exact spacing between the sound instances, which I assume you know in the following. If you don't know it, you can find it from calculating the autocorrelation sequence of your composed signal.

Without further knowledge, the only linear method for finding your sample is to project the composed signal onto the periodic subspace it lives in. There are two obvious ways you can do this.

The first is taking the Fourier series of the composed signal and nulling all coefficients that do not belong to an integer multiple of the frequency corresponding to the period length of the repetition. Summing the modified series will give you the desired signal.

The second is to just accurately cut the sample into segments that align perfectly with respect to the repetition period and just sum all those segments up into a single sample.

The result of these two operations will give you a representation of the composed signal the same periodic subspace. However, since the subspace is relatively large compared to the total signal, there will be components of the composed signal that are contained in this subspace by mere chance. So the separation will be far from perfect. You can improve the separation by taking more repetitions, make the full signal space much larger than the periodic subspace.

If you need a better method you won't get around characterising the nature of the signal you want to extract. That's where perception models, source models, statistical signal models, etc come in. I cannot really go into this in this answer. It's a topic that fills many books and is in general very difficult.


Find B given A and impulse response of B/A

The vector impr holds the impulse response of length nh of an all-pole filter with denominator coefficients given by a. Let's call this impulse response $g(n)$. The impulse response $h(n)$ must satisfy the equation

$$h(n)=(g*b)(n)\tag{1}$$

where * denotes convolution. With the given data, this results in an overdetermined system of linear equations, the system matrix of which is the Toeplitz matrix egqs. The last line in the code solves this overdetermined system of linear equations corresponding to (1) in a least-squares sense. You could achieve the same using the backslash operator:

b = egqs\h';


Spectral location of individual OFDM symbols

The IFFT transforms a discrete signal to the time domain. Therefore, the input signal of the IFFT should be interpreted as a frequency domain signal.

For OFDM this means that the subcarriers are defined before the IFFT. The IFFT ouput is the time domain representation that will eventually be transmitted over the channel. Before it can be transmitted it must first be converted to an analog signal. Analyzing this signal in time domain, e.g. by viewing it on an oscilloscope, won't tell you anything about the subcarrier bandwidth.

To see the bandwidth of one OFDM subcarrier you have to analyze the analog signal in the frequency domain, e.g. by using a spectrum analyzer. If you modulate all subcarriers the spectrum will look nearly rectangular, with a flat top, because the subcarriers are close together and overlap. If you modulate only one subcarrier (all inputs of the IFFT are set to zero except one), the spectrum will be much narrower, its bandwidth approximately reduced by a factor that is equal to the number of subcarriers. And it will have a sinc-like shape.


Determine whether a signal is periodic or not and get fundamental period

$\cos\left(\dfrac{n}{2}\right)*\cos\left(\pi*\dfrac{n}{4}\right)$

Now finding out period of $\cos\left(\dfrac{n}{2}\right)$
$\cos\left(\dfrac{n}{2}\right)=\cos\left(2*\pi*\dfrac{n}{N1}\right) \implies \dfrac{1}{2}=2*\dfrac{\pi}{N1} \implies N1=4*\pi$

Now finding out period of second function $\cos\left(\pi*\dfrac{n}{4}\right)$
$\cos\left(\pi*\dfrac{n}{4}\right) = \cos\left(2*\pi*\dfrac{n}{N2}\right)\implies \dfrac{\pi}{4}=2*\dfrac{\pi}{N2} \implies N2=8$.

Now for total function period is LCM of $(N1,N2)$ i.e, $\mathrm{LCM}(4*\pi,8) \implies32*\pi$.
Therefore, the period of the function is $\boxed{32*\pi}$.


So it seems that the formulas from the textbook for the bandpass are not accurate.

I used following approach to calculate second order sections given an analogue filter lowpass. As an example following code displays the coefficients for a 4th order bandpass, calculated from a second order lowpass protopye.

Paramters

fs = 1000;
fs1 = 200;
fs2 = 300;
w1 = 2*pi*fs1/fs;
w2 = 2*pi*fs2/fs;
Omega1 = 2*tan(w1/2);
Omega2 = 2*tan(w2/2);
syms s z Om2 Om1 a1 b1 A0
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Start with analogue transfer function H_s

H_s = A0 / (1 + a1*s + b1*s^2);
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Digital conversion with bilinear transformation

G_z_BP = subs(G_s_BP, s, 2*(( z - 1)/(z + 1)));
% order polynomial such that coefficients can be read out
coeffs_z = collect(G_z_BP, z)

Those coeffs_z are very similar to the ones given by the iowahills software
The band-pass filter coefficients then can be calculated with any analogue low-pass prototype, just enter numerical values on read out coefficients from coeffs_z


Consider a continuous-time signal $$x_a(t) = A \cos(\Omega_0 t + \theta) $$ where $\Omega_0$ is the frequency, and $\theta$ is the phase.
To perform spectral analysis of $x_a(t)$, sample it with a period of $T_s$ resulting in the discrete-time sequence $$x[n] = A \cos( \omega_0 n + \theta)$$ where $\omega_0 = \Omega T_s$ is the discrete-time frequency in radians.
Consider the relationship between the CTFT $X_a(\Omega)$ of $x_a(t)$ and the DTFT $X(e^{j\omega})$ of $x[n]$ for $\omega \in [-\pi,\pi]$
$$X(e^{j\omega}) = \frac{1}{T_s} X_a(\frac{\omega}{T_s})$$
Computation of $X(e^{j\omega})$ of $x[n]$ may require indefinetely long duration of $x[n]$, such as infinite duration if $x_a(t)$ is an ideal sine wave. And also $X(e^{j\omega})$ is a function with continuous domain of $\omega$. Thus it's not practical to compute the DTFT.
However, one can compute samples $X[k]$ of $X(e^{j\omega})$, known as the DFT of $x[n]$, defined over a finite length of $x[n]$ obtained by applying a window: $$v[n] = w[n]x[n]$$ where $w[n]$ is the window of length $L$.
The consequence of this windowing on the result of the computed spectrum is: $$V(e^{j\omega}) = \frac{1}{2\pi} X(e^{j\omega}) \star W(e^{j\omega}) $$
Computing the $N$-point DFT of the sequence $v[n]$ therefore provides the samples $V[k]$ of the result of the periodical convolution between the true spectrum $X(e^{j\omega})$ of $x[n]$ and the DTFT $W(e^{j\omega})$ of the window $w[n]$. Note that $N \geq L$.
For the specific example we have considered, $$x[n]=A\cos(\omega_0 n + \theta)$$ whose DTFT is $$X(e^{j\omega})= A\pi e^{j\theta} \delta(\omega-\omega_0) + A\pi e^{-j\theta} \delta(\omega+\omega_0)$$ a pair impulses weighted by $A\pi$.
And assuming a rectangular window of length $L$ for $w[n]$, its DTFT is: $$W(e^{j\omega})= \sum_{n=0}^{L-1} w[n] e^{-j\omega n} = e^{-j{\omega (L-1)/2}} \frac{\sin(\omega L / 2)}{\sin(\omega / 2)} $$ the result is $$V(e^{j\omega}) = 0.5 A e^{j\theta} W(e^{j(\omega - \omega_0)}) + 0.5 A e^{-j\theta} W(e^{j(\omega + \omega_0)})$$
Finally the DFT samples $V[k]$ are given by the following:
$$V[k] = V(e^{j\frac{2\pi k}{N} }) = 0.5 A e^{j\theta} W(e^{j(\frac{2\pi k}{N} - \omega_0)}) + 0.5 A e^{-j\theta} W(e^{j(\frac{2\pi k}{N} + \omega_0)}) $$
One peak peak of these DFT samples would occur for $\omega = \omega_0$ for some $k$ if $\frac{2\pi k}{N} = \omega_0$ is satisfied for an integer $N$ and $k$. Assuming $T_s$ ,$L$,$N$,and $\Omega_0$ are selected such that the equality is satisfied for some integer $k_0$, then the magnitude of the peak of the observed spectrum is:
$$ |V[k_0]| = | 0.5 A e^{j\theta} W(e^{j(0)}) + 0.5 A e^{-j\theta} W(e^{j(2\omega_0)}) | $$
Solving for $A$, the amplitude of the original analog signal whose spectrum analysis are after, yields:
$$A = \frac{2 |V[k_0]|}{|e^{j\theta} W(e^{j(0)}) + e^{-j\theta} W(e^{j(2\omega_0)}) |} $$
An approximation to this can be obtained by assuming, for the rectangular window, $W(e^{j 2\omega_0}) \approx 0$ and noting that $W(e^{j0}) = L$ : $$A \approx \frac{2 |V[k_0]|}{L}$$
If a window other than rectangular is used, formulation should be adjusted accordingly, specifically $W(0)$ being the sum of the window samples, which is $L$ for the rectangular window.


Modeling spring oscilliations with electric noise (50 hz)

The problem is that just taking the x[0] and x'[0] of your observed data and plugging them into the differential equation will be extremely sensitive to any kind of noise ; not only the 50Hz hum. You are trusting those two samples (rather than the entirety of your observed data) for your estimation, and you will get incorrect results from that, whether the 50Hz noise is present or not.

How to do it right? First, get rid of the 50Hz noise with a notch filter tuned to 50Hz and maybe to a few harmonics.

Then, what you want to do is fit the parameters $A$ and $\phi$ of a $x(n) = A e^{-\delta n} sin(\omega n + \phi)$ model to your data $y(n) = x(n) + noise$. You already know $\delta$ and $\omega$ since they depend on your known system parameter (damping and oscillation frequency). What you really want to know is the best estimate of $A$ and $\phi$ to reconstruct the "theoretical" part (in proper speak: signal subspace projection) of your observation.

Here is how you can estimate $A$ and $\phi$ without lame methods involving looking at peaks and zerocrossing:

You can write your signal model as $x(n) = Re (\alpha z^n) = \frac{1}{2}\alpha z^n - \frac{1}{2} \bar{\alpha} \bar{z}^n$ where $z$ is the known complex pole $z = e^{-\delta + j\omega}$, $\alpha$ is the complex amplitude $\alpha = Ae^{j \phi}$ and $\bar{z}$ denotes the conjugate.

Lay out a $N\times2$ matrix where N is your number of observation, the first column being $\frac{1}{2} [1, z, \ldots, z^N]$ and the second $-\frac{1}{2} [1, \bar{z}, \ldots, \bar{z}^N]$. Do a least square regression of that matrix (which is a basis of your signal subspace) onto the data matrix $[y(1), \ldots, y(N)]$. This will give you LS estimates of $\alpha$ and $\bar{\alpha}$ as the two solutions. Take the modulus and angle of the first to get $A$ and $\phi$. Plug this back into $A e^{-\delta n} sin(\omega n + \phi)$ to get the best fit theoretical model from your data.

# Estimate amplitude and phase of an exponentially damped model
# given observations y and system parameters delta and omega
n = numpy.arange(y.shape[0])
z = numpy.exp(-n * delta + n * omega * 1j)
z_conj = numpy.vstack((0.5 * z, -0.5 * numpy.conj(z))).T
s, _, _, _ = numpy.linalg.lstsq(z_conj, y)
print "Estimated amplitude: ", numpy.abs(s[0])
print "Estimated phase: ", numpy.angle(s[0])

If you want to also estimate $\delta$ and $\omega$ the procedure is a tad more complex and involves an EVD or SVD of the correlation matrix of your data to get a signal subspace basis, one more EVD to get the estimate of $z$, then the procedure sketched above to retrieve $A$ and $\phi$. You will have a robust estimate of all the parameters of the theoretical model from your data. Search for ESPRIT or MUSIC.

Note that this approach can be used to estimate the phase / frequency of a mixture of any number of damped sine waves. So for example, you could use this to retrieve the amplitude and phase of the 50Hz hum itself. Just put two additional columns in the matrix used for the least square procedure, with powers of $e^{j\omega_{noise}}$ and their conjugates. This will give the best fit amplitudes and phases of an exponentially damped sine wave (described by omega and delta) + 50 Hz hum model.


Improve coverage of BackgroundSubtractorMOG2 and findContours

I'm not an expert in this field, but have you tried eroding and then dilating?
I'm writing a python application to track tennis players in the court, and what I did in order to get a single contour for the players is to erode and dilate the tracked shapes.

erosion = cv2.erode(img,None,iterations = 1);
dilation = cv2.dilate(erosion,None,iterations = 1);
cv2.findContours( dilation.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE);

or you can directly call the function

opening = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel);
cv2.findContours( opening.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE);

Source: OpenCv Morphological Transformations


How to convert from Laplace Domain to Time Domain?

As per: https://math.stackexchange.com/questions/3485053/how-to-get-partial-fractions-for-this-equation-to-go-from-laplace-to-time-doma

$F(s) = \frac{-2 m k v R}{2mR\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

$F(s) = \frac{-2 m k v R}{2mR}*\frac{1}{\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

$F(s) = - k v *\frac{1}{\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

Then I can use:

$\frac{1}{(s+a)(s+b)} = \frac{e^{-at}-e^{-bt}}{b-a}$

And that should do it!


Discrete-time Biquad filter: Relation between peak frequency and pole frequency

According to the Audio EQ Cookbook the transfer function of a second-order band-pass filter with 0dB peak gain is given by

$$H(z)=\frac{\alpha-\alpha z^{-2}}{(1+\alpha)-2\cos(\omega_0)z^{-1}+(1-\alpha)z^{-2}}\\=\frac{\alpha}{1+\alpha}\cdot\frac{1-z^{-2}}{1-\frac{2}{1+\alpha}\cos(\omega_0)z^{-1}+\frac{1-\alpha}{1+\alpha}z^{-2}}\tag{1}$$

where $\omega_0$ is the angular frequency where the peak occurs ("peak frequency"), and $\alpha=\sin(\omega_0)/(2Q)$ is a parameter determined by $\omega_0$ and the desired quality factor $Q$.

If we express the same transfer function $H(z)$ in terms of its poles we get

$$H(z)=k\frac{1-z^{-2}}{\left(1-re^{j\omega_p}z^{-1}\right)\left(1-re^{-j\omega_p}z^{-1}\right)}=
k\frac{1-z^{-2}}{1-2r\cos(\omega_p)z^{-1}+r^2z^{-2}}\tag{2}$$

where $k$ is some gain constant to be determined, $0<r<1$ is the pole radius, and $\omega_p=2\pi f_p$ is the pole angle, corresponding to the pole frequency $f_p$. Comparing the denominator of $(2)$ to the filter coefficients given in this question (and shown under Eq. $(1)$ in the question above), it is clear that the frequency $f_0$ in that question equals the pole frequency $f_p$ in our formulation. And $r$ in the question above is obviously the pole radius.

Comparing $(1)$ and $(2)$ we can express $\alpha$ and $\omega_0$ of Eq. $(1)$ in terms of the pole radius $r$ and the pole angle $\omega_p$ used in Eq. $(2)$:

$$\begin{align}&\frac{1-\alpha}{1+\alpha}=r^2  \Longrightarrow \alpha=\frac{1-r^2}{1+r^2},\qquad
k=\frac{\alpha}{1+\alpha}=\frac{1-r^2}{2}\\
&\frac{2}{1+\alpha}\cos(\omega_0)=(1+r^2)\cos(\omega_0)=2r\cos(\omega_p)\end{align}\tag{3}$$

The second line of Eq. $(3)$ shows that the pole angle $\omega_p$ and the peak frequency $\omega_0$ are not equal, and that the following relationship holds:

$$\frac{1+r^2}{2r}\cos(\omega_0)=\cos(\omega_p)\tag{4}$$

The following Matlab/Octave example illustrates the above. We choose $\omega_0=0.2\pi$ and $Q=1$ and design the corresponding biquad band-pass filter according to Eq. $(1)$. The figure below shows the magnitude of the resulting frequency response, and the locations of the peak frequency $\omega_0$ (red line) and of the pole angle $\omega_p$ (green line), both normalized by $\pi$.

w0 = .2*pi;
Q = 1;
al = sin(w0)/(2*Q);
b = [al,0,-al]/(1+al);                % numerator coeffs
a = [1+al,-2*cos(w0),1-al]/(1+al);    % denominator coeffs
p = roots(a);                         % poles
r = abs(p(1));                        % pole radius
r2 = sqrt((1-al)/(1+al));             % pole radius according to formula
abs(r-r2)                             % ans = 0
wp = abs(angle(p(1)));                % pole angle: wp = 0.17877*pi
wp2 = acos((1+r^2)/(2*r)*cos(w0));    % pole angle according to formula
abs(wp-wp2)                           % ans = 0
[H,w]=freqz(b,a,1024);
plot(w/pi,abs(H))
hold on, plot([w0,w0]/pi,[0,1],'r',[wp,wp]/pi,[0,1],'g'), hold off

EDIT:

In Eq. $(2)$ I assumed that there are two complex conjugate poles, because in the original question, which triggered this question and answer, the filter coefficients were chosen in such a way that only complex conjugate poles were possible (cf. the coefficients under Eq. $(1)$ of the question above).

However, the general transfer function in Eq. $(1)$ can also have a double real-valued pole or two different real-valued poles. It can be shown that a double real-valued pole occurs when the quality factor $Q$ is chosen as $Q=\frac12$. This case is still covered by Eq. $(2)$ with $\omega_p=0$ (i.e., a positive real-valued pole) or $\omega_p=\pi$ (a negative real-valued pole). For $Q<\frac12$ we get two different real-valued poles, and Eq. $(2)$ is not valid anymore. Instead, the transfer function becomes

$$H(z)=k\frac{1-z^{-2}}{(1-\beta_1z^{-1})(1-\beta_2z^{-1})}=
k\frac{1-z^{-2}}{1-(\beta_1+\beta_2)z^{-1}+\beta_1\beta_2z^{-2}}\tag{a}$$

where $\beta_1$ and $\beta_2$ are the two real-valued poles. Comparing this transfer function to the transfer function given in Eq. $(1)$, we get the following relationships between $\alpha$, $\omega_0$, $k$ and $\beta_1$ and $\beta_2$:

$$\begin{align}\alpha&=\frac{1-\beta_1\beta_2}{1+\beta_1\beta_2}\\
\cos(\omega_0)&=\frac{\beta_1+\beta_2}{1+\beta_1\beta_2}\\
k&=\frac{1-\beta_1\beta_2}{2}\end{align}\tag{b}$$

Of course, the pole angles are now either zero or $\pi$, because both poles are real-valued. Note that the peak frequency $\omega_0$ is not  necessarily $0$ or $\pi$. The peak frequency $\omega_0$ can in fact take any value, and the two pole angles will always be either $0$ or $\pi$ as long as the quality factor $Q$ satisfies $Q<\frac12$.


What is Finite Rate of Innovation Signal?

If a signal can be exactly represented by $N$ real numbers per time interval, then its number of degrees of freedom for that time interval equals $N$. The most well-known example are band-limited signals, which can be represented by $2B$ samples per second, if their bandwidth is not greater than $B$. I.e., such signals have $2B$ degrees of freedom per second. This is equivalent to saying that these signals have a finite rate of innovation, and that rate of innovation equals $2B$.

A signal doesn't need to be band-limited in order to have a finite rate of innovation (i.e., a finite number of degrees of freedom). You can think of any signal that can be parameterized with a finite number of parameters per time unit. A simple example would be a signal of the form

$$x(t)=\sum_kw_kg(t-kT)\tag{1}$$

where $g(t)$ is a rectangular impulse with a value of $1$ in the interval $[0,T]$ and zero otherwise. You need exactly one sample per $T$ seconds to represent $x(t)$, even though $x(t)$ is clearly not band-limited.

A smoothing kernel is just an optional filter that is applied to the continuous-time signal before sampling. For sampling a (quasi-)band-limited signal, the anti-aliasing filter would be an instance of such a smoothing filter.

An annihilating filter for a specific signal is an LTI system that produces an output of zero when that specific signal is applied at its input. E.g., for a sinusoidal signal with frequency $\omega_0$, a filter with a notch at $\omega_0$ is an annihilating filter.


How to get coordinates of complex RTL-SDR signal?

You are missing several steps.

You have (complex) time-domain samples.  You need to convert that to scalar frequency domain data to get something like your example plot.  Many steps are required to do that conversion before there are any frequency peaks to be found.

You first need to:  Pick a segment of complex IQ data from the RTL-SDR of length N (you may need to convert the raw RTL-SDR IQ samples from unsigned 8-bit to signed floating point), window it (von Hann or Hamming, etc.), convert it to the frequency domain via an FFT of length N, convert the FFT result to log magnitude, and label the FFT log magnitude result bins by frequency, which will be roughly

frequency(i) = sdr.center_freq + i * sdr.sample_rate / N

for bins 0 to N/2

frequency(i) = sdr.center_freq - (N - i) * sdr.sample_rate / N

for bins N/2 to N-1

Then you can search along that log magnitude array for peaks, and apply the frequency label to the peaks found.


Which information do we get from magnitude and phase spectrum?

You get the exact same information you will receive if you analyse a 1-D signal using the Fourier analysis tools, for example. To illustrate this consider the following examples

We perform Fourier Transform on it and obtain the following spectrum

As you can see there are two symmetric dots representing the frequency that is present in the image and the central, DC, component.
Next up

Again the Fourier Transform

It is obvious that the higher frequencies are further away from the DC, central, component. To understand the logic behind this you shall remember that we are actually working with complex numbers and their properties kick in.
To answer your second question: You are usually working with sparse data structures (vectors, arrays, tensors, etc.) in the frequency domain. Now you tell me which is easier: Working on the original data set consisting of 500 different values or the transformed one with half of them 0s ?
Now, the importance of phase - drum rolls - it is the same, meaning HUGE
Again, a sample image

Fourier Magnitude Spectrum

Fourier Phase Spectrum

Now, we inverse the transformation by using just the amplitude and then just the phase information
Inverted amplitude spectrum

Inverted phase spectrum

I think you get the idea :D
Now, here are some sample images and their respective amplitude and phase spectra, so you can practice a wee bit


How to recover data/information from real complex baseband signal?

In an actual received signal you will also have to address the frequency, phase and time offsets between the transmitter and receiver since they run off of independent clocks that aren’t otherwise synchronized, as well as frequency offsets introduced through Doppler if the transmitter and receiver are in motion relative to each other. Further information on the details of this can be found under the topics of Carrier Recovery and Timing Recovery.
For CDMA this can be done with a course frequency search by doing the OP’s correlation repeatably at different frequency offsets and using an acquisition threshold on the correlator output to determine the rough offset. Once within the acquisition bandwidth of the correlator (approximately 1/T where T is the duration of the PN sequence in time), the complex I Q outputs of the correlator can be used to measure phase change per correlation interval (which as  $f=d \phi /dt$  is the fine frequency measurement (carrier offset).  This approach can be part of a frequency lock loop or the individual phase measurements can be part of a phase lock loop; either being the carrier recovery implementation.
With CDMA one approach to timing and carrier recovery (common in GPS receivers) is to use three correlators each using the matching PN sequence separated by half a chip such that the correlator outputs form “Early”, “Prompt”, and “Late”.  The prompt is used for carrier acquisition and tracking as described above by both sweeping frequency for carrier offset and the PN code for delay offset (at each frequency test). Once acquired, the early and late can maintain tracking (time offset) by subtracting the two (zero time error when early minus late equals zero).
I have other posts detailing this further as well as timing recovery that I will add links here. The above approach is sequential so slow but minimizes resource requirements.  One link I show also provides further details on a direct (one-shot) frequency and delay (PN offset) acquisition that would be relatively fast (can be completed on a single capture) when the processing resources allow.
See:
Early - late gate spread spectrum
High modulation index PSK - carrier recovery
QPSK constellation from baseband signal


importance of symmetric windows for STFT for practical applications

This example may help you to understand this.

The plot below shows the step response for two filters that have the same magnitude response, but the first has non constant group delay (non linear phase) while the second has constant group delay (linear phase).

The first filter is a minimum phase version of the second filter.


Background Noise Reduction in Audio signal

Depends a bit on the type of noise. If it's hum or buzz, you need to clean up your recording electronics and if push comes to shove you may be able to notch it out.

Real acoustic background noise is difficult to remove. You may want to read up on "blind source separation". Here are a few examples http://cnl.salk.edu/~tewon/Blind/blind_audio.html


Removing a sinusoidal artifact from a set of movie frames

Your proposed solution - calculating a sinusoid in the time domain based on the peak in the FFT, then subtracting it - should work, but there's an easier way to do essentially the same thing: modify that peak value in the FFT, then take the inverse transform.

So, for your rasterized video M[nframes, npixels] you find the frequency bin holding the artifact, then systematically flatten it (e.g., set its magnitude to the average of its neighbors) for every pixel:

import numpy as np
nframes, npixels = np.shape(M)
# Identify the bin containing the sinusoidal artifact
# Use the average intensity for each image
m = np.mean(M, axis=1)
# Calculate the FFT
Fm = np.fft.rfft(m)
# Find the largest bin away from the low-frequency region
lowfreq = 100  # or something
badbin = lowfreq + np.argmax(Fm[lowfreq:]**2)

# Now adjust the amplitude of that bin in the FFT of each pixel
for pixel in range(npixels):
   Fpix = np.fft.rfft(M[:, pixel])
   # Scale magnitude of artifact bin to be the mean of its neighbors
   Fpix[badbin] *= np.mean(np.absolute(Fpix[[badbin-1, badbin+1]]))/np.absolute(Fpix[badbin])
   # Rewrite the time sequence of that pixel
   M[:,pixel] = np.fft.irfft(Fpix)

This should work if the artifact is exactly constant amplitude and frequency, and its frequency falls right on to a submultiple of the sequence length (i.e., the sinusoids represented by the FFT).  In general, you might want to flatten out one or two bins either side of badbin to deal with a slightly broader set of narrowband corruptions, e.g.

# ...

   # Scale magnitude of artifact binS to be the mean of neighbors
   spread = 3  # flatten bins from (badbin - (spread-1)) to (badbin + (spread-1))
   # target value for new bins
   targetmag = np.mean(np.absolute(Fpix[[badbin-spread, badbin+spread]]))
   bins = range(badbin - (spread-1), badbin + spread)
   Fpix[bins] *= targetmag/np.abs(Fpix[bins])
   # ...

If you want to constrain the component removed from each pixel to have the same frequency and phase of the artifact detected in the mean intensity, you could remove just the projection of the badbin magnitude onto that phase, e.g.

badbinphase = np.angle(Fm[badbin])
# ...

   Ncomponent = np.abs(Fpix[badbin])*np.cos(np.angle(Fpix[badbin]) - badbinphase)
   Fpix[badbin] -= Ncomponent * np.exp(0+1j * badbinphase)
   # ...

Note that the resulting component at badbin will now always be 90\deg phase shifted (orthogonal) from the global badbinphase in every pixel - any signal component at exactly that frequency and phase cannot be separated from the artifact.


Time Setting of $z$ and Laplace Transforms

In typical (or all) applications, the Laplace Transform is used for continuous time systems and the z-Transform for discrete systems. However, the Laplace Transform for Discrete Time Systems certainly exists but would be more complicated than necessary to solve. The z transform exists as a mathematical simplification of the Laplace Transform that can be applied if a consistent sampling rate is used. This may be clearer with the formulas given below:

The (one-sided, for causal systems) Laplace Transform is given as:

$$X(s) = \int_{t=0}^\infty x(t) e^{-st}dt$$

Thus the Laplace Transform for a discrete time system is given simply by setting t= nT for integer n (where the integral then becomes a summation):

$$X(s) = \sum_{n=0}^\infty x(nT)e^{-snT}$$

At this point we can make a substitution to get rid of the more complicated exponential (which is one possible mapping from s to z):

$$z = e^{sT}$$

resulting in the z-Transform for a discrete-time sequence:

$$X(z) =  \sum_{n=0}^\infty x(nT)z^{-n}$$

Equivalently the z-Transform for a discrete sequence (no units of time, simply sample number) is:

$$X(z) =  \sum_{n=0}^\infty x(n)z^{-n}$$

Thus we see that the Laplace Transform of a discrete time system exists (but we have no need to solve that typically as working with the z transform once discrete is much simpler). I am not aware of a reverse case although we may be able to create one by taking the limit $T\rightarrow 0$ of the discrete time z-Transform. I haven't worked through this but we may find the same situation (we can force it to be mathematically equivalent but the math involved would be more complicated than necessary vs taking the Laplace Transform directly).


How to remove lines from captcha image which are acting as noise using java?

I worked on this problem some time ago, my solution was to compress the image with wavelet transform with an opportune coefficient that depends on captcha type that I choose with some expermients.
Then I saturate the compressed image from mean value to 0 and 1 in order to create a mask.
Next I multiply mask to the original image and I get the image without lines or noise. I post a picture example where I use also a derivative filter on the result but in reality you don't need it.

Your lines are a bit different but I think that you can start from this approach because your lines are more thin than text so wavelet compression can works but you have to improve it for your problem.

Finally your lines have different color and you can use also this information.


Histogram operation for improving visual object tracking

I don't think you have any choice other than to use the same number of bins for each observation.  Otherwise not only will you not be able to average the histograms, you will also not be able to compare them.

And you definitely need to change the histogram slowly, i. e.
$$h = (1 - \alpha)h + \alpha h_{obs}$$ where $h$ is your "moving average" histogram, $h_{obs}$ is the histogram of the current observation, and $\alpha$ is a small factor, whose value you will have to tune experimentally.

As far as different sizes of observations, the simplest thing to do is to impose a hard constraint saying that the size of the object cannot change too much from frame to frame.


Why discrete time signals are defined as sequence of numbers?

Samples of (regularly sampled) discrete time signals, like pixels of a digital images, are sequences: they are unitless on the ordinal axis (albeit time or space), and unitless in amplitude. A sequence is

an enumerated collection of objects in which repetitions
are allowed and order matters

This is close to what people call raw data.
However, to make a genuine "signal" or "image" object (meaningful in the real world or for devices) that can be better interpreted or manipulated, it is important to keep ancillary information, like: sampling rate, quantization, time-space coordinates, size, channels, amplitude units, calibration parameters, etc. Those additional parameters are found for instance  in the header of most digital sound or image raster formats, and then followed by the unitless digital sequence.
Below is an example from the WAVE PCM soundfile format (converted in png for display). The Discrete signal is typically in the bottom "data" sub-chunk, yet requires the above information to be read, listened to, analyzed.


how can I analize music as spotify/musixmath/shazam does?

The general idea for the algorithm is well presented in the linked paper. As is the case with applied machine learning and related fields, the secret sauce (which features to select, how to tweak your algorithms, etc.) is found by trial and error by people who know things about audio/audio signal processing (see for example, Julius O. Smith III's page at CCRMA @ Stanford) and by just trying a ton of things. Unfortunately, to become (and stay) the state of the art, the secret sauce has to be kept, well, secret.

That being said, the wikipedia article on acoustic fingerprinting is a decent place to look for some alternatives (and some links to some open source ones like Acoustid).


What is the value of the given periodic signal at any time, how to solve it in MATLAB?

Assuming continuous-time for the probem, for all periodic signals with a period of $T$ , for any integer $m$ , we have :
$$  x( t + mT) = x(t) \tag{1}$$
So, for example with  $T=7$ , and $m=2$ , you can see that $x(15) = x(2\cdot 7 + 1) = x(1)$. The operator that is used to find base argument for a given $t$ value is the modulus :
$$ x(t) = x( \text{mod}(t,T) ) \tag{2}$$
where $\text{mod}(t,T)$ is the modulus of $t$ wrt $T$.
For your example, your period is $T=7$, and you can see that $x(315) = x( \text{mod}(315,7) ) = x(2)$.
Hence the value of $x(t)$ for $t=2$ is $0.1 \cdot 2^3 = 0.8$.
I don't know if there's anything to solve with MATLAB here though?


The Fourier and the Laplace transform obviously have many things in common. However, there are cases where only one of them can be used, or where it's more convenient to use one or the other.

First of all, even though in the definitions you simply replace $s$ by $j\omega$ or vice versa to go from one transform to the other, this cannot generally be done when given the Laplace transform $X_L(s)$ or the Fourier transform $X_F(j\omega)$ of a function. (I use different indices because the two functions can be different for the same time domain function). There are functions for which only the Laplace transform exists, e.g., $f(t)=e^{at}u(t)$, $a>0$, where $u(t)$ is the Heaviside step function. The reason is that the integral in the definition of the Laplace transform only converges for $\Re\{s\}>a$, which implies that the corresponding integral in the definition of the Fourier transform does not converge, i.e. the Fourier transform doesn't exist in this case.

There are functions for which both transforms exist, but $X_F(j\omega)\neq X_L(j\omega)$. One example is the function $f(t)=\sin(\omega_0t)u(t)$, for which the Fourier transform contains Dirac delta impulses.

Finally, there are also functions for which only the Fourier transform exists, but not the Laplace transform. This means that the integral in the definition of the Laplace transform only converges (in a specific sense) for $s=j\omega$, but for no other values of $s$. The Laplace transform is only said to exist if the integral converges in a half-plane or in a vertical strip of finite size of the complex $s$-plane. Such functions for which only the Fourier transform exists include complex exponentials and sinusoids ($-\infty<t<\infty$), and impulse responses of ideal brick-wall filters, which are related to the sinc function. So, e.g., the functions $f(t)=\sin(\omega_0 t)$ or $f(t)=\sin(\omega_ct)/\pi t$ do not have a Laplace transform but they do have a Fourier transform.

The Laplace transform can be a convenient tool for analyzing the behavior of linear time-invariant (LTI) systems by considering their transfer function, which is the Laplace transform of their impulse response. The poles and zeros of the transfer function in the complex $s$-plane conveniently characterize many system properties and are useful for an intuitive understanding of the system's behavior. Furthermore, the unilateral Laplace transform is very useful for analyzing LTI systems with non-zero initial conditions. The Fourier transform is a useful tool for analyzing ideal (non-causal, unstable) systems, such as ideal low pass or band pass filters.

Also have a look at this answer to a related question.


Linear classification and overlapping problem reduction

It's very difficult to tell what you're asking exactly, but one thing to keep in mind is that not al data is linearly separable. You may have to obtain more features from data or try non-linear separation techniques. SVM is one of them, but since you're not considering it, you can try polynomial separation or logistic classifiers.


In my experience in the defense industry, the quality assurance engineers are mostly just “box checkers” that have limited understanding of how all of the individual facets of the algorithms and systems work. We have design engineers, integration engineers, and then the QA engineers. I have never seen an integration engineer or QA engineer get involved in any sort of design; the skill sets for each position are very, very different. It’s one thing to know what something is, but it’s something else entirely to be able to use it in a creative way to solve a problem.

Again, this is just my experience in the defense industry within the US. Other countries and industries may be a bit different though!


Is it possible to identify the "blank data" in a JPEG without knowing the image history?

The short answer is that the original encoding of a picture is lost when you re-encode it.  So a tool that takes a 70% encoded picture and re-encodes it to 100% actually does the following:

Decode from 70% to an internal, raw, uncompressed format (similar to .bmp -- basically, it'll be encoded in RGB, probably 8 bits per pixel for each color).
Encode to 100%

Note two things here: that the encoding assumes nothing about the original image, and that changing the encoding of an image (even from "perfect" RGB to "100% quality" JPEG) is going to introduce scrud -- numerical errors that are for all practical purposes impossible to predict and degrade the image.
So your 70% image isn't perfect -- compared to the original it has numerical scrud.  That -- plus a bit more because the conversion isn't perfect -- will be carried forward into the intermediate uncompressed image.
Then, during "compression" to 100% quality any numerical scrud that's in the intermediate uncompressed image will be treated as real image, not as numerical scrud.  So when the intermediate raw image -- complete with numerical scrud -- is converted to 100% image quality jpeg, that numerical scrud will be replicated as faithfully as can be managed.  Note that, because of the number of steps that the JPEG process goes through, even the best possible job you can do, that results in no compression at all, will still introduce its own scrud.
It actually gets worse than that, but not the way that your question addresses.  If you take an image and compress it, then uncompress it, then compress it again, the second compressed image will be lower quality than the first.  I honestly don't know if for certain levels of compression you'd reach some minimum of image quality at which the thing would no longer degrade, but I'm sure that if you were to compress to some slightly different image quality each time (i.e., choose a random number between 70% and 80% each time), that the image would just continue to degrade with compression -- even as the file sizes pretty much stayed the same for each compression level.


What plotting / visualization tools are useful in generating good plots?

I guess Python has great tools for this, especially IPython Notebook: http://ipython.org/notebook.html . It has good visualization opportunities and at the same time you could document your code and make it re-configurable and online. Of course it best for visualizing data low dimensions and creative visualization algorithms can be prepared either in OpenGL or high level languages such as MATLAB.

In 3D, if appropriate, time varying surfaces can be embedded in 2D and shown as a 3D volume. For generic 4D visualizations you could either use videos, or accumulated volume maps. For visualization of high dimensional data, always make sure you check methods such as Multi Dimensional Scaling.


Just to expand what I wrote in the comment. This is how you should work with the FFT in Matlab.

Compute the FFT of your time signal:

NFFT = 2^nextpow2(N);
Y = fft(y,NFFT)/N; % the division by N is to scale the amplitude

Then compute the sampling frequency: (EDIT skip this 2 lines of code if you already know the value of $Fs$)

dt = t(2)-t(1);
Fs = 1/dt;

This works if your data are sampled at constant sampling frequency. If it is not the case you need to resample them (it is very easy to do if you convert your data into timeseries and then use the resample function for timeseries. Look here for more details: What is an algorithm to re-sample from a variable rate to a fixed rate?)

Finally, plot the result:

f_plot = Fs/2*linspace(0,1,NFFT/2+1); % vector of frequencies for the plot
figure()
plot(f_plot,2*abs(Y(1:NFFT/2+1))) % multiplicated by 2 to recover the energy related
% to the negative frequencies since in this way only half spectrum is plotted
xlabel(’Frequency’)


Can error in digital signal be zero after filtering?

The discrete signal is a sum of two sines. The noise is an additive uniform noise, always positive (rand lies in open interval $(0,1)$). The latter  is the crucial point with the median. Suppose you have a frame of three signals samples: $x_0$, $x_1$, $x_2$, and three noise realizations:  $n_0$, $n_1$, $n_2$ in $[0,1]$. So you are checking if:

$$\textrm{med}\{x_0+n_0,x_1+n_1,x_2+n_2\}-x_1$$

can be zero for each window frame. You have noticed from the plots that the noisy signal is always above the clean one. So the  median of $\{x_0+n_0,x_1+n_1,x_2+n_2\}$  is always higher than the  median of $\{x_0,x_1,x_2\}$.

But could it be the same as the initial data (putting signal borders away)? Yes. Take  $\{x_0,x_1,x_2\}=\{0,0.5,0.2\}$ and add the random realization $\{0.5,0.1,0.3\}$. Then the median of the noisy signal is $0.5$, the exact value of the center sample $x_1$.

So for some well-chosen signals and noise, their are counter-examples, a median filter can in cases provides a zero error. However, and I admit I cannot give a proof at this time, it is probably very unlikely that there is a zero error for  long enough generic signals.

Note to DSP folks: I feel like I am missing a very evident proof, or such a question is not adapted at this level of signal processing.


When you compute cross correlations, you shift one signal over the other and compute the (normalized) inner product of the overlapping sections of the two signals. The point is that you can shift in both directions, i.e., your time shift (in samples) can be positive or negative.

Assume both signals have $N$ samples. Your central data point in the cross correlation is when you haven't shifted your second signal at all, i.e., all $N$ points are overlapping. Now, you can shift "to the right" by 1 sample: your overlap is $N-1$. Shift another one: your overlap is $N-2$. Keep going, until you shifted $N-1$ samples: there is only an overlap of 1 sample. This gives one half of your correlation function. You get the other half by shifting in the other direction, "to the left", you get another $N-1$ samples like that. In total, this gives $2N-1$ samples, which explains the "doubling" you observe.

Example: Correlate the sequences [1,2,3] and [4,-5,6]. Normalization left aside, we get as inner products:

[1,2,3]'*[0,0,4] = 12
[1,2,3]'*[0,4,-5] = -7
[1,2,3]'*[4,-5,6] = 12
[1,2,3]'*[-5,6,0] = 7
[1,2,3]'*[6,0,0] = 6

Therefore, from N=3 original lags, we get 2N-1=5 lags in the correlation function.

If you have a good reason to believe that your signals are periodic, you can use a cyclic convolution, which will result in N samples exactly. But his only makes sense for periodic signals.


Please clarify the mandatory condition on Vestigial Sideband Modulation filter?

For $(F_c,F+dF)$, it needs to be $2A-H(-F)$, when $H(F_c)=A$. See the picture below taken from (http://163.152.6.72/contents/vod/21/20121129093403/20121129093403.pdf). Here $A$ corrsponds to $1/2$ as per below figure, that is the value at $H(F_c)$. In order for the sum to become $1$, $A$ needs to $1/2$. Also, $\epsilon$ is value at any frequency in the range $(F-dF,F_c)$. So the value in $(F_C,F+dF)$ needs to be $2A-H(-F)$ so that sum $H(F) + 2A-H(-F) = 1$, as $H(F)=H(-F)$ along $F=0$.


Implementation of modulation/demodulation models

I would suggest you to use the GnuRadio API usualy used for Software Defined Radio (SDR) projects.
As Marcus Müller said in the comments, it is easier to develop this kind of project and test it with a software implementation.
I am also pretty sure that you will find a lot of existing examples to help you with your project.
Guided Tutorials
API documentation (not beginner-friendly)
Source code


It is quite common to combine amplitude and phase modulation. QAM is actually a good example because generally you will have symbols with different phases and different amplitudes. Just in the special case of $4$-QAM (which is actually the same as QPSK), all symbols have the same amplitude. But if you take as an example $16$-QAM then you'll see that the information is encoded in the phase as well as in the amplitude:

[Chris Watts [CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0)], from Wikimedia Commons]

Note that in general the exact position of the complex symbols is a compromise between ease of implementation and power efficiency. QAM is especially easy to implement and the penalty compared to more power efficient constellations is relatively small.


I can't interpret numpy code, but I don't think you should define an 'N' variable. For a 2-second signal you should have a total of 16000 samples. Here's how to do this in MATLAB (notice how I defined my "number of samples" variable ('Num_Samples')):
Fs = 8000;
Ts = 1/Fs; % Time between samples
T1 = 2; % Signal duration = two seconds
F0 = 0; % Start frequency
F1 = 10; % Stop frequency (at T1 seconds)
Num_Samples = 2/Ts; % Two seconds divided by 'Ts'
n = 0:Num_Samples-1; % Here's your time-domain index
Alpha = (F1-F0)/T1;
X = cos(2*pi*F0*n*Ts + Alpha*pi*(n*Ts).^2);

figure(1)
plot(n*Ts, X, '-bs'), grid on

Finding the number of zero-crossings in the 2 second duration of your chirp signal is nothing more than examining each generated 'x(nTs)' sample (starting with the 2nd sample) and seeing if its polarity is different from the previous 'x((n-1)Ts)' sample. Each time a polarity difference occurs you have a single zero-crossing.


Proof of fourier transformation of multiplication of two signals

I am very thankful to Dilip-sarwate and Gilles, who took their precious time to understand my problem and guide me.
So, Now I'm going to write the correct solution to my question. Which is as follows:
 $$
ℱ[x(t)g(t)] = \int_{t=-\infty}^{\infty}[x(t)g(t)]e^{-j\omega t}dt
$$
As we know :
$$
x(t) = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)e^{j\alpha t}d\alpha
$$
Therefore,
$$
ℱ[x(t)g(t)] = \int_{t=-\infty}^{\infty}\bigg[\frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)e^{j\alpha t}d\alpha\bigg]g(t) e^{-j\omega t}dt
$$
After rearranging, we have
$$
ℱ[x(t)g(t)] = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)\bigg[\int_{t=-\infty}^{\infty}g(t)e^{-j(\omega-\alpha)t}dt\bigg]d\alpha
$$
As we also know that :
$$
G(\omega) = \int_{t=-\infty}^{\infty}g(t)e^{-j\omega t}dt
$$
and replacing $\omega$ with $\omega-\alpha$, we get:
$$
G(\omega-\alpha) = \int_{t=-\infty}^{\infty}g(t)e^{-j(\omega-\alpha)t}dt
$$
Therefore,
$$ℱ[x(t)g(t)] = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)G(\omega-\alpha)d\alpha$$
Finally, we get
$$
ℱ[x(t)g(t)] = \frac 1{2\pi}[X(\omega)*G(\omega)]
$$

[$\because$ $x_1(t)*x_2(t) = \int_{\tau=-\infty}^{\infty}x_1(\tau)x_2(t-\tau)d\tau$]


No you can't, because the internet depends on you having a method of sending data. For example, the question you've asked here came here somehow, so it's obvious you need to be able to send data.

Now, TV is broadcasting, it's a one-way transmission, without a back channel.

There are satellite internet access services, but they are far more costly than what you've described, and much much slower.

Also, some countries had things that used a telephone line as uplink and a sattelite downlink, but they suffered from cost, the asymmetry between the fast satellite downlink and the slow telephone uplink, and the poor latency.


What is it called when you have a single-pole recursive digital lowpass filter where the coefficients don't add up to unity?

Stability only depends on $b_1$, and it needs to satisfy $|b_1|<1$ for stability of the recursion.


At complex baseband (meaning the signal is represented as I+jQ), the signal will have the same magnitude and phase versus time as the passband signal at any carrier given the system is linear.  For example, the signal 1+j0 has magnitude 1 and angle 0; we can also have a carrier at 1 GHz (or any frequency) with the envelope of magnitude 1 and phase 0 (relative to a reference at the same frequency). Likewise shifting the baseband by 45 degrees: .707 + j.707 results in magnitude 1 and angle 45 degrees is the same as shifting the 1 GHz signal by 45 degrees. Any linear operations on the two will be the same and it is significantly easier to model the baseband signal of interest versus keeping track of all those carrier cycles!  This is not very different from observing a spinning bike wheel with a strobe light, if we set the strobe light to be the nominal rpm of the spinning wheel, we can readily observe small changes in phase from nominal.

If the idea of a "complex baseband" is confusing, please also review this post here:

Frequency shifting of a quadrature mixed signal

This is also demonstrated in the following figures:

The first picture shows a complex baseband signal both in the time domain as a phasor moving randomly in magnitude and phase, but bandlimited such that it cannot move quickly. In the frequency domain we see a spectrum with bandlimited white noise (the spectral density is constant up to the bandwidth of the system).

Next we see the concept of frequency translation, where shown that multiplying by $e^{j\omega_o t}$ simply translates our spectrum to a new carrier $\omega_o$ but is otherwise unchanged.  The first plot below shows the exponential and it's Fourier transform, illuminating the frequency translation property since the transform is simply a frequency translated impulse: multiplying by the exponential (spinning phasor) in the time domain is convolution in the frequency domain, convolving with an impulse in frequency does not change any of the signal characteristics except for changing it's carrier frequency.

The next plot below shows the original signal we introduced, both at a carrier frequency and then "down-converted" to baseband by multiplying with $e^{-j\omega_o t}$, a minus sign in this case to move the signal from $\omega_o$ to baseband.  Notice specifically in this picture that the magnitude and phase shown versus time, on the IQ diagram, is identical for both: Both show the magnitude and phase relative to their respective carriers- in the plot on the left the carrier frequency is at $\omega_o$ and the plot on the right the "carrier" is at 0.  The strobe light analogy I gave earlier can also help with seeing this more intuitively.

Since these waveforms are identical under linear operations, I would usually choose to model and describe the waveform at baseband rather than at any actual carrier frequency, meaning in a simulation for example, I am not simulating an actual carrier waveform of which the signal of interest is modulating its amplitude and phase, but would do that directly at complex baseband with the same result.


LTI system output given input and frequency response

Apply the eigenfunction property of the LTI system after decomposing the sinusoidal input by Euler identity.

The eigenfunction property of the LTI system states that

$$ x(t) = e^{j \omega_0 t} \implies y(t) = H(\omega_0) e^{j \omega_0 t} $$

where $H(w)$ is the frequency response of the LTI system.

Euler identity states that :

$$ x(t) = \cos(t) = 0.5 \{ e^{jt} + e^{-jt} \} $$

Then due to system being LTI, you can define the output as:

$$ y(t) = 0.5 \{ H(1) e^{jt} + H(-1) e^{-jt} \} .$$

Since the system impulse response was a real function, its frequency response $H(w)$ will be conjugate symmetric; i.e., $H(w) = H^*(-w)$.

Let's call $H(1) = A e^{j\phi} $ where $A = |H(1)|$ is the magnitude and $\phi = \angle{H(1)}$ is the phase of the frequency response at the frequency $\omega = 1$.

Hence the output will be :
$$ y(t) = 0.5 \{ H(1) e^{jt} + H^*(1) e^{-jt} \} $$

consider the term $z = H(1)e^{jt}$  as a complex variable, then the second term will be $z^*$, hence the output is:

$$y(t) = 0.5 \{ z + z^* \} = \mathcal{Re} \{ z \} $$

$$y(t) = \mathcal{Re} \{ A e^{j\phi} e^{jt} \} = A \cos(t + \phi) = |H(1)| \cos(t + \angle{H(1)})$$


I'm a little confused by this statement given the previous one indicates you've already taken the FFT.

  Finally, I take each of the blocks (computed by the STFT) and perform PSD.

Generally the following is used; Seperate into frames (overlap if desired) -> apply window -> fft -> transform and scale into dB (if desired).

The scaling for a power spectrum or power spectral density can be a bit tricky at first, but this paper goes into some great detail about how to set it up.


Find highest frequency of a very manky signal

I think most of your questions can be solved by subtracting then mean from the signal. Namely the all sinusoidal waves, with a nonzero frequency, have a mean of zero. So when the mean of a signal is not close to zero, then it will show up at 0 Hz (you can look at this as $\cos(0\,t)=1$).

After a closer look and some testing by myself it does seem that only removing the mean would not fix your problem. Your signal seems to be contaminated by white noise, which has been integrated probably twice (so low frequencies are much more present than high frequencies). You could try to subtract a higher order polynomial fit, so linear or even quadratic. These fits can easily be obtained using some linear algebra, for example for a quadratic fit:

X = [ones(N, 1), t', t'.^2];
b = (X' * X) \ X' * y';
z = y - (X * b)';

Or you could also look at the biggest peaks above $f^{-2}$, with $f$ the frequency vector. However if you have a signal with normal white noise you should probably just stick with removing the mean, since my proposed methods try to make use of the fact that the present noise is white noise which has been integrated twice.


Discrete Fourier series for a sum of deltas

$\sum\limits_{r=-\infty}^{\infty}\delta[n-rN]$ is a periodic Kronecker delta with period $N$, but the sum to get the fourier coefficients only "uses" one period $\sum\limits_{n = 0}^{N-1}$ so you don't have to worry about the periodic extension of $\delta$ and take just $x[n] = \delta[n]$ as your function. This is going to cancel out every exponential in the sum except when $n=0$, so the result will be independent of the particular coefficient and be just $a_k = 1/N$.


You can evaluate $X[k]$ (and $x[n]$) outside the domain $[0,N-1]$. The sequences are periodic with period $N$, so $X[N]=X[0]$, and, more generally, $X[k]=X[k+N]$.


Estimating the time offset (delay) between two audio signals in real-time

There are a variety of techniques available. See this thesis for a nice write-up.

See the list of algorithms:

In particular, note that they recommend pre-whitening the signals before correlation --- to ensure that the cross-correlation has a nice, narrow peak for better "localization".


Notations to Use in Formulating of Maximum Likelihood Estimation

1) What notations to use for probability density function is it the one below: $\mathsf{P}_y(y_n|{\mathbf{u}_n})$ $\mathsf{P}_z(z_n|{\mathbf{u}_n})$ what goes in the subscript if I want to use $z$?

One way to write the density function is to subscripted upper case letters for the random variables. So $p_{Y_n}$ denotes the density of random variable $Y_n$ and $p_{Z_n}$ is the density of the random variable $Z_n$. Similarly, $p_\mathbf{Y}:=p_{Y_1,\ldots,Y_n}$ denotes the joint density of all the $Y_i$'s.

  2) Is $P(\mathbf{y}|\mathbf{\theta})$ correct or
  $P_y(\mathbf{y}|\mathbf{\theta})$ or
  $P_z(\mathbf{z}|\mathbf{\theta})$?

See my answer to (1) above. In some cases it is ok to drop the subscript if it is understood which random variable the density $p$ belongs to. However, if you have lots of different random variables it is good to explicitly include the subscripts.

  If I want to use expectation maximization, then would $\xi$ which
  denotes the complete data set consisting of the input signal and the
  received noisy signal be written as:
  $\xi = \{u_0,\ldots,u_{N-1},y_0,\ldots,y_{N-1} \}$ or
  $\xi = \{\hat u_0,\ldots,\hat u_{N-1},y_0,\ldots,y_{N-1} \}$

It is not clear what the hats mean.

  4) Would the posterior be written as

You need to clarify the posterior of what conditioned on what. I assume you are interested in estimating the $u$'s so the posterior you are referring to is

\begin{eqnarray}
p_{\mathbf{U}|\mathbf{Y}}(\mathbf u|\mathbf y) &\propto& p_{\mathbf{Y}|\mathbf{U}}(\mathbf y|\mathbf u) p_\mathbf{U}(\mathbf u) \\
&=& p_{Y_1,\ldots,Y_N|U_1,\ldots,U_{N}}(y_1,\ldots, y_N|u_1,\ldots u_N) \; p_{U_1,\ldots,U_{N}}(u_1,\ldots u_N)\\
&=& \left(\prod_{i=1}^N p_{Y_i | U_{i},\ldots U_{1 \wedge (i-p+1)}} (y_i|u_{i},\ldots u_{1\wedge (i-p+1)}) \right) \; p_{U_1,\ldots,U_{N}}(u_1,\ldots u_N).
\end{eqnarray}

That $\wedge$ symbol denotes a min. operation and makes sure that the subscripts don't go below 1. If you have a probabilistic model for the $U_i$'s you can further factor out $p_{U_1,\ldots,U_{N}}$ into a product of conditional densities. For example, if it's a 1st order Markov model with a stationary initial distribution $p_{U_1}(\cdot)$ you can write,
$$
p_{U_1,\ldots,U_{N}}(u_1,\ldots,u_N) = p_{U_1}(u_1)\;\prod_{i=2}^N p_{U_i|U_{i-1}}(u_i|u_{i-1}).
$$

  pdf, likelihood, log likelihood

Finally, I think, you may be interested in writing the "complete data" log-likelihood function i.e. the density of $\xi$ parameterized by $\theta$. This can be written as follows:

\begin{eqnarray}
p_{\xi|\theta}(\xi|\theta) &=& p_{U_1,\ldots ,U_N,Y_1,\ldots ,Y_N|\theta}(u_1,\ldots u_N, y_1,\ldots y_N | \theta) \\
&=& p_{Y_1,\ldots ,Y_N | U_1,\ldots ,U_N,\theta}(y_1,\ldots y_N | u_1,\ldots u_N,  \theta) \; p_{U_1,\ldots ,U_N|\theta}(u_1,\ldots u_N| \theta)
\end{eqnarray}

and simplify this further using your observation model for $Y_i$'s and Markov transition model for $U_i$'s.

It seems that your confusion stems from being unfamiliar with how random variable notation works. I would recommend reviewing a textbook on probability and random processes for electrical engineers.


You want a time-varying low pass filter, automated so that the cutoff frequency follows your curve. For such slow modulations most IIR low pass filters are stable, so you can just find a steep enough low pass filter and update its coefficients while filtering while keeping the filter state.


Creating an Eye Diagram for QPSK Simulation

Your eye diagrams for the inphase and quadrature baseband signals are perfectly correct. For QPSK with rectangular signals and perfect matched filtering, the inphase and quadrature signals are indeed just BPSK with rectangular signals, and the matched filter outputs can be deduced from, for example, the figure (reproduced below) at the end of this answer of mine.

If you sketch the matched filter output between $T$ and $5T$, overlaying the signal in chunks of length $2T$ centered at $2T$ and $4T$, you will see that the signal is flat at the lower level when two $1$s are transmitted in succession. Similarly, the signal will be flat at the upper level if two $0$'s are transmitted in succession. For transitions between two data bits, the matched filter output crosses over.

As to where the smooth transitions instead of straight ramps come from, if the modulating signal is not a rectangular pulse, or if the matched filter is not quite matched perfectly, the filter output can well be smoothly transitioning from one level to the other. For example, if the BPSK pulse were one cycle of a sinusoid of period $T$, specifically
$$s(t) = \begin{cases}\sin\left(\frac{2\pi t}{T}\right), &0 \leq t < T,\\
0, &\text{otherwise,}\end{cases}$$
then the matched filter output would be a perfect sinusoid of period $T$ between two successive sampling instants if the two data bits were the same, and make a smooth transition from peaking at $nT$ to peaking with opposite polarity at time $(n+1)T$.
See, for example, the figure below.


Could there be any reason to prefer convolution-based calculation of autocorrelation?

And as it is well known, the computational complexity is $\mathcal{O}(N^2)$ and $\mathcal{O}(N \log N)$ for them respectively.

$\mathcal{O}$ notation ignores any constants that determine exactly how fast those functions run. Depending on the constants out front, it may be faster to compute the convolution (or correlation) directly.

For example, the Fourier transform method to convolve may take $10^5 N \log N$ seconds while the direct method may take $10^{-3} N^2$ seconds. If $N = 100$, the direct method is faster.

In general for large $N$, the Fourier method tends to be faster with speedups of up. For more detail on computing which method is faster, see scipy PR #5608.

Below is a graph that times the Fourier transform method of convolution as fftconvolve and the direct method as convolve. We see that the direct method is faster for many arrays of practical size.

Note: this graph is when convolving 1D signals. For more higher dimensional signals (2D, 3D, etc) fftconvolve tends to be faster more often

To time these graphs, I used scipy.signal.fftconvolve and np.convolve with the default parameters. As of scipy version 0.19, scipy.signal.convolve will choose between the direct and Fourier transform methods to choose the faster method.


If you look at Gary Bishop, Greg Welch - An Introduction to the Kalman Filter (SIGGRAPH 2001, Course 8) then, on page 17, there are equations 3.3 and 3.4:

This is probably not the KF formulation you normally see because the system equation, 3.3, is kind of lagged one behind the observation equation, 3.4, time-wise ( with respect to their right hand sides ). This is less common than the one where the RHS's of both equations has the same time subscript. So, without going into all the gory details, if you derive the KF based on 3.3 and 3.4, you get totally different recursions but they are equivalent to those in the KF formulation where the RHS's of the observation and system equations are contemporaneous rather than lagged by one. Hopefully this gives the general idea. I would have put it in comments since it doesn't provide the details but it wouldn't fit.


The issue is possibly that the input signal you have chosen is not persistently exciting.  This means that the signal doesn't "excite" enough modes of the filter in order to be able to accurately estimate its parameters.  Another way to think about it is that it doesn't have enough energy in enough places in the spectrum: just at the frequency of the cosine, whereas white noise has energy everywhere in the spectrum.

From the paper:

  Recent work (e.g. [1–5]) has shown the convergence and robustness properties of adaptive estimation
  and control algorithms to be closely related to persistence of excitation of regression vectors (equivalently
  that they have a sufficient number of spectral lines [5])


Convert a 1D array into a 'rotated' 2D matrix

I assume that the length of the array is odd and that the array is symmetrical. Also note that from your requirements not all elements of the resulting matrix are specified. Those elements will be zero in the code shown below.

v=[1,2,3,2,1];      % odd-length array
A=diag(v);
A=A+rot90(A);
n=(length(v)+1)/2;
A(n,:)=v;
A(:,n)=v;
A =

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


Difference between real fft and complex fft with imaginary part of zero in fftw?

I haven't tried this myself but I think it's just a problem of ordering the output data. In the real-to-complex case, your output matrix is not a $3\times 3$ matrix, but - due to the implied symmetry - a $3 \times 2$ matrix. So re-arranging your output gives

 36            -4.5+2.59808i
-13.5+7.79423i  0
-13.5-7.79423i  0

where the "missing" third column would just be the complex conjugate of the second column. This is also explained here.


Wiener filter to estimate constant signal

You can use this code

 function w=WienerFIRFilter(input,desired,M)
     auto_corr=xcorr(input,M,'unbiased');
     r=auto_corr(M+1:end); % positive lags only
     R=toeplitz(r); %correlation matrix
     cross_corr=xcorr(input,desired,M,'unbiased');
     p=(cross_corr(M+1:end));
     w=inv(R)*p'
 end


BPSK modulation under bandlimited channel

1. Matlab can't handle symbolic calculations by itself. You need to use a symbolics-capable package, such as Matlab's own symbolic toolbox, Sage, SymPy, Maple, Mathematica, etc.

2. No, it is not correct. The sequence of -1 and 1 is a sequence of symbols, but what is actually transmitted over the channel is an analog signal. Between the symbols and the channel you need a pulse-shaping filter, which determines the transmitted signal's bandwidth and spectral characteristics.

3. The receive process will depend on the channel's frequency response. Say your transmitted signal has bandwidth $B$. If the channel response is (mostly) flat over that bandwidth, then the channel doesn't distort the signal and you can recover the information present in the signal by matched filtering. On the other hand, if the channel's response is not flat, then the signal will be distorted. In this case, you can perform channel equalization followed by matched filtering. A channel equalizer is a filter that is the "inverse" of the channel; the equalizer's output is an undistorted signal.

I've used a few concepts without explaining them, like pulse shaping and matched filtering. There are already many answers in this site regarding those concepts. If you need further help with these, please leave a comment.


Before I address your questions, you should understand:

a. the integral branch of the loop filter maintains a average phase increment in units of radians/sample.  It is not a frequency value, though in the right context it can be converted to a frequency value, using your sample rate as a conversion factor.

b. the total output of the loop filter is an instantaneous phase increment which is used to advance the numerically controlled oscillator's phase.

So to your questions:

  Why is the input phase error of the Costas loop limited in the [1,-1] region?

It looks like an arbitrary limit to keep the instantaneous phase increment output of the loop filter to some sane number of radians/sample.  (Correcting more than $\pm\pi$ radians/sample is pointless.)

The average phase increment maintained by the integral branch is limited to $\pm 1.0$ radians/sample as well.  In a no-signal, noise processing situation, the integrator output can grow without bound, so the limiters are needed.  The $\pm 1.0$ radians/sample is just an arbitrary limit to keep the loop near sane operating values for when a signal does come in.  A Costas loop that is applying a correction of more than 1 radian/sample is probably not anywhere near converged anyway.

  How do you choose the frequency region [max_freq and min_freq] of the integrator part of the loop? why do we need to limit the frequency range?

In your application context, which you really haven't provided, you have to decide what are reasonable limits for the average phase increment coming out of the integrator.  You must provide limits, because the integrator is unstable and its output can grow without bound, which would make the PLL useless.

For a Costas loop, the spacing in radians of your constellation points and the number of samples/symbols will dictate what are sensible limits for the maximum & minimum increment of the NCO in radians for each sample processed.  Anything outside of the interval $(-\pi, \pi)$ radians/sample would be nonsensical due to phase wrapping.

For a PLL that is trying to lock on to a tone in a particular frequency range $[f_{low}, f_{high}]$ in Hz, you can use your sample rate, $F_s$, to convert those to limits in units of radians/sample, i.e. $\left[\dfrac{\pi f_{low}}{F_s/2},\dfrac{\pi f_{high}}{F_s/2}\right]$ radians/sample.

  Why is the output phase Φ_o[n] (in both Costas and FLL) wrapped between -2PI and 2PI?

Because computation of $e^{j\phi_o}$, used for setting the NCO output, is generally numerically better with smaller angle arguments in most libraries.  Phase is cyclic so it's OK to just wrap down to small magnitude angles.


Do $|s(t)|$ and $|S(f)|$ uniquely determine $s(t)$?

As pointed out in a comment by Jazzmaniac, the proof below is flawed. The mistake I made was the assumption leading to Eq. $(3)$. Of course, $|S_1(f)|=|S_2(f)|$ implies that the two spectra are related by an all-pass transformation, but the phase of that all-pass can be dependent on the input signal, and hence need not be LTI.

One example of two signals satisfying $|s_1(t)|=|s_2(t)|$ and $|S_1(f)|=|S_2(f)|$ without satisfying $(1)$ is $s_1(t)=e^{j\phi(t)}$ and $s_2(t)=e^{-j\phi(-t)}$. Clearly, $|s_1(t)|=|s_2(t)|=1$. Since $s_2(t)=s_1^*(-t)$, the corresponding spectra satisfy $S_2(f)=S_1^*(f)$, and, consequently, $|S_1(f)|=|S_2(f)|$.

The dual case would be $S_1(f)=e^{j\Phi(f)}$ and $S_2(f)=S_1^*(-f)=e^{-j\Phi(-f)}$. The signals $s_1(t)$ and $s_2(t)$ then satisfy $s_2(t)=s^*_1(t)$ and $|s_1(t)|=|s_2(t)|$.

The second example above requires the signals to be complex-valued (otherwise $s_1(t)=s_2(t)$ would hold). However, the first example also includes the real-valued case, i.e., a signal $s_1(t)$ that only attains two values: either $1$ or $-1$ (like a telegraph process). With $s_2(t)=s_1(-t)$ we guarantee $|s_1(t)|=|s_2(t)|$, and since $S_2(f)=S_1(-f)=S_1^*(f)$ we also have $|S_1(f)|=|S_2(f)|$.

OLD answer:

I think I can show that the trivial case

$$s_2(t)=e^{j\alpha}s_1(t),\qquad\alpha\in\mathbb{R}\tag{1}$$

is the only case for which $|s_1(t)|=|s_2(t)|$ and $|S_1(f)|=|S_2(f)|$ holds.

From the requirement $|s_1(t)|=|s_2(t)|$ it follows that

$$s_2(t)=e^{j\phi(t)}s_1(t)\tag{2}$$

with a real-valued function $\phi(t)$. From the requirement $|S_1(f)|=|S_2(f)|$ it follows that $s_2(t)$ can be obtained from $s_1(t)$ by all-pass filtering:

$$s_2(t)=\int_{-\infty}^{\infty}s_1(\tau)h_{AP}(t-\tau)d\tau\tag{3}$$

with $H_{AP}(f)=\mathcal{F}\{h_{AP}(t)\}$ satisfying $|H_{AP}(f)|=1$.

Since the multiplication in $(2)$ is memoryless, Equations $(2)$ and $(3)$ can only give the same result if $h_{AP}(t)=e^{j\alpha}\delta(t)$, i.e., for a trivial all-pass which just scales the input by a complex constant with magnitude equal to $1$. Consequently, the function $\phi(t)$ in $(2)$ must be constant (and equal to $\alpha$). This means that $(2)$ and $(3)$ are reduced to Eq. $(1)$.

Consequently, given $|s_1(t)|$ and $|S_1(f)|$, the signal $s_1(t)$ is determined up to a complex factor with magnitude $1$. I just do not know of any method to actually compute $s_1(t)$ from these two magnitude functions.


Question about zeroes of simple 2nd order FIR filter

The reason for this is that the maths centres on a false assumption:

  We know that the zeroes need to be complex conjugates

Do the zeros need to be complex conjugates? Can they not be 2 real zeros?

I sympathise here because I also read this book and the maths portion ended up utterly confusing me because of very odd methods.

A more general method would be at this point:

$0=1+\alpha_1z^{-1}+\alpha_2z^{-2}$

Rearrange thus:

$0=\frac{z^{2}+\alpha_1z^{1}+\alpha_2}{z^{2}}$

Here on the top we have a quadratic equation, which you can use any quadratic solver for.

For your first case, I arrive at 0.635 $\pm$ 0.637i, which is a pair of complex numbers denoting the conjugate zeros - correct.

For your second, I arrive at two real zeros, one at 3.36 and one at 0.24 - so it is most definitely able to be computed!

And finally, you will have noted the $z^{2}$ term lying on the bottom, which indicates that there is a square pole at 0+0i. If these zeros were poles then there would similarly be a square zero at 0+0i.

All of these results can be verified with the zplane function in MATLAB.

It took me some time with this book and checking results with the zplane function in matlab to realise that his methods of working really don't work on a general basis, and frankly overcomplicate things. It has however been awhile since I've looked at it so maybe I similarly misinterpreted his maths