Can I use any kind of transform to make spectrogram?

If your input signal is a finite length sequence and you want to analyze how the frequency domain components are varying in that input sequence, you have to use one of the techniques used for time-frequency analysis. STFT is one such technique, there are others also like wavelet transform, quadratic time-frequency distribution functions, etc.
For the $2^{nd}$ question, if you take FFT of the entire signal as 1 window, then those FFT coefficients will show all of the existing frequency components at once, you won't have the time-information. Dividing the input into smaller chunks of overlapping windows and then taking FFT of each window provides you with time information too, as in what frequency component was significant at what time, during the whole signal collection.


recovery plane from 4 point imaging using calibration camera

I found the answer in the paper: Calibrating an Overhead Video Camera by Raul Rojas in section 3 - 3.3.

for the start:
 H=K^-1*H

1)  Given four points in the image and their known coordinates in the world, the matrix H can be recovered, up to a scaling factor . We know that the first two columns of the rotation matrix R must be the first two columns of the transformation matrix. Let us denote by h1, h2, and h3 the three columns of the matrix H.

2)  Due to the scaling factor we then have that xr1 = h1 and xr2 = h2
Since
|r1| = 1,
then

x = |h1|/|r1| = |h1| and

x = |h2|/|r2| = |h2|.

3) We can thus compute the factor and eliminate it from the recovered matrix H. We just set

H'= H/x

4) In this way we recover the first two columns of the rotation matrix R. The third column of R can be found remembering that any column in a rotation matrix is the cross product of the other two columns (times the appropriate plus or minus sign). In particular

r3 = r1 × r2

5) Therefore, we can recover from H the rotation matrix R. We can also recover the translation vector (the position of the camera in field coordinates). Just remember that

h'3 = −R^t

6) Therefore the position vector of the camera pin-hole t is given by

t = −R^-1 h3


First, your supposed 'radix-4 butterfly' is a 4 point DFT, not an FFT.  It has 16 complex (ie: N squared) operations.  A typical 4 point FFT would have only Nlog(base 2)N (= 8 for N = 4).  Second, you have some supposed w[ ].r and w[ ].i 'scale' factors that don't belong.  Perhaps you obtained them from a radix-4 butterfly shown in a larger graph.  Such a butterfly would have some interstage twiddles appended to it, but they are not actually part of the butterfly.  A 4 point FFT only has an internal butterfly of -j when designed for a negative exponent FFT.

Rather than try to fix your code, it's just as easy to write my own, as shown below (DevC++ compiler; outputs appended at end of code):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft4(double* r, double* i);    // prototype declaration
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O

double r[4] = {1.5, -2.3, 4.65, -3.51}, i[4] = {-1.0, 2.6, 3.75, -2.32} ;
long n, k, j;      double  yr[4] = {0.}, yi[4] = {0.};
double ang, C, S, twopi = 6.2831853071795865;

cout<<"\n original real/imag data";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

// 4 point DFT
for (k = 0; k < 4; k++) {
    ang = twopi*k/4;
    for (j = 0; j < 4; j++) {
        C = cos(j*ang);       S = sin(j*ang);
        yr[k] = yr[k] + r[j]*C + i[j]*S;   // ( C - jS )*( r + ji )
        yi[k] = yi[k] + i[j]*C - r[j]*S;   // = ( rC + iS ) + j( iC - rS )
    }
}

cout<<"\n 4 point DFT results";
cout<<"\n n         yr[n]           yi[n]           amplitude       phase(radians)\n";
double amp, phase;
for (n = 0; n < 4; n++)  {
    yr[n] = yr[n]/4 ;      yi[n] = yi[n]/4 ;  // scale outputs
    amp = sqrt( yr[n]*yr[n] + yi[n]*yi[n] ) ;
    phase = atan2( yi[n], yr[n] ) ;
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,yr[n],yi[n],amp,phase);
} //end for loop over n

fft4(r, i) ;

cout<<"\n 4 point FFT results";
cout<<"\n n         r[n]            i[n]            amplitude       phase(radians)\n";

for (n = 0; n < 4; n++)  {
    r[n] = r[n]/4 ;      i[n] = i[n]/4 ;  // scale outputs
    amp = sqrt( r[n]*r[n] + i[n]*i[n] ) ;
    phase = atan2( i[n], r[n] ) ;
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,r[n],i[n],amp,phase);
} //end for loop over n

fft4(i, r); // this is an inverse FFT (complex in/out routine)

cout<<"\n 4 point inverse FFT results";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

system ("PAUSE");
return 0;
} // end main
//************************ fft4 **********
void fft4(double* r, double* i) {
double t;

t = r[0]; r[0] = t + r[2]; r[2] = t - r[2];
t = i[0]; i[0] = t + i[2]; i[2] = t - i[2];
t = r[1]; r[1] = t + r[3]; r[3] = t - r[3];
t = i[1]; i[1] = t + i[3]; i[3] = t - i[3];

t = r[3]; r[3] = i[3]; i[3] = -t; // (r + ji)*(-j)

t = r[0]; r[0] = t + r[1]; r[1] = t - r[1];
t = i[0]; i[0] = t + i[1]; i[1] = t - i[1];
t = r[2]; r[2] = t + r[3]; r[3] = t - r[3];
t = i[2]; i[2] = t + i[3]; i[3] = t - i[3];

t = r[1]; r[1] = r[2]; r[2] = t;  // swap 1
t = i[1]; i[1] = i[2]; i[2] = t;  //  and 2
} // end fft4

 original real/imag data
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

 4 point DFT results
 n         yr[n]           yi[n]           amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point FFT results
 n         r[n]            i[n]            amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point inverse FFT results
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

First, the input data (4 real, 4 imaginary) is printed out.  Then a 4 point DFT is taken.  The results (yr[ ] and yi[ ] plus amp/phase) are printed out.  Since the r[ ] and i[ ] original data were not overwritten when doing the DFT, those inputs are reused as inputs to the 4 point FFT.  Note that the latter has fewer +/- operations than the DFT.

The code for the FFT is not particularly elegant nor efficient – there are many ways of doing butterflies.  The code above corresponds to the four radix-2 butterflies shown in Rabiner and Gold's book “Theory and Application of Digital Signal Processing” (p. 580, Fig. 10.9), with twiddles modified to reflect a negative exponent (the ones used for the figure in the book were positive).  Note that there is only one twiddle of -j in the code, and this does not require a multiply (it's a swap/sign change).

After the FFT, the results are printed.  They're the same as the DFT

And finally, the scaled results from the FFT are used as inputs to an inverse FFT.  This is accomplished via the 'exchange' or 'reverse the list' method (ie: if FFT(r,i) is a forward FFT, then FFT(i,r) is an inverse – provided, of course, that the FFT is capable of handling complex inputs/outputs – in other words – no 'real only' routines, which usually presume that imaginary inputs are zero).  This method was described nearly 25 years ago in:

P. Duhamel, B. Piron, J. M. Etcheto, “On Computing the Inverse DFT,” IEEE Transactions
on Acoustics, Speech and Signal Processing, vol. 36, Feb. 1988, pp. 285-286.

The result of the inverse is then printed out.  It is the same as the original input data.


Butterworth-like filter with no zero at Nyquist?

Is the attenuation caused by aliasing or something different, and what is it called?

The two standard methods of mapping an analog filter to a digital one are either the bilinear transform or the impulse invariant transform. Neither does a particularly good job at preserving the transfer function at higher frequencies. The bilinear transform will map whatever happens at "infinite" analog frequency to the digital Nyquist frequency. For any analog lowpass filter (Butterworth or not), this will indeed create the steep roll off at Nyquist.

Are there digital alternatives to the Butterworth or bilinear transform that do not suffer from the attenuation near Nyquist, or is this a usual feature?

There are plenty of alternatives. Impulse invariant is one of them, although it to tends overshoot by 6dB at Nyquist. One heuristic method is to start with bilinear, and then move the zeros from z = -1 towards z = 0 until you have something you "like" by whatever criteria is important for your application.

What filter could an old resource-strained synthesizer have used to achieve a response without the attenuation?

A biquad is a biquad is a biqaud. Computational cost is not a function of the specific filter coefficients. There is really no difference between a Butterworth or any other second order filter.


How can I sample a signal at 4 samples per cycle?

You need to understand that $f=f_s/4$ and "$4$ samples per cycle" are two ways of saying the same thing:

$$\text{ # samples per cycle}=\frac{f_s\text{ samples per second }}{f \text{ cycles per second}}$$

The number of cycles per second is equal to the ratio $f_s/f$, where $f_s$ is the sampling frequency, and $f$ is the signal's frequency.

So you just need a sampled sinusoid with a frequency $f=f_s/4$ (or, equivalently, $4$ samples per cycle). The value of $f_s$ is irrelevant.

$$x[n]=\sin(2\pi f\cdot n/f_s)=\sin(n\pi/2)$$

which is exactly what you have already come up with. Multiply this signal with your low pass impulse response, and you'll get a band pass filter centered at $f_s/4$.


Have you fftshifted your results or not? The answers will be different:

If you fftshift your $256$ x $256$ 2-dimensional Fourier transform result, retain only the $64$ x $64 $ central co-efficients, inverse fftshift, and then inverse fourier transform the result, the effect will be a blurring of your original image.

This is because you have culled our your images' high frequency components, hence, the effect will be a blurring.

If however you have not fftshifted, then you are doing the opposite, and will actually be doing a high pass filtering. (Your image will certainly not be blurred).

By the way, hard thresholding fourier co-efficients like this is certainly used, but will give you ringing effects. However, those effects might be negligible and/or removable based on your application.


Why do digital modulation when the signal entering the channel is analog?

I think part of your confusion is actually between discrete-time and digital signals. Some textbooks do not clearly distinguish between them.

All signals that exist in the physical world are real. So, the signals that are produced by any transmitter, and the signal that is picked up by the receiver, are real.

In communications, we say we do "digital communications" when the analog transmitted signal conveys "symbols", which are actually numbers drawn from a finite set. This is different than "analog communications", such as AM and FM radio, where the traansmitted information is a (audio) waveform that can take arbitrary values.

In addition to this, these days communications systems are designed as much as possible using digital signal processing techniques. The nature of these techniques is that they deal with discrete-time (sampled) signals. So, a modern software-defined radio receiver will have an analog antenna, a few analog filters and amplifiers, and an analog quadrature downconverter; then, the analog signal is discretized (sampled) and handed off to a digital signal processor for further processing. The digital signal processing can happen on specialized hardware (such as an FPGA or ASIC), or as software running on a DSP or on a general-purpose CPU. This is where digital filters would be used.

Note that a receiver can be implemented using DSP regardless of whether the underlying communications technology is analog (AM, FM, old TV, etc) or digital (BPSK, QPSK, etc).


This signal is already after AM demodulation

That means this is only meaningful if the signal actually only contains information in its amplitude; luckily, most simple remotes just use some form of On-Off-Keying or Pulse Position Modulation.

  Is it possible to determine the encoding type only by looking at the signal?

If this really is just a signal that is either on or off, then this looks like on-off keying. But really, your visualization might have hidden a lot of relevant information.

This might be a lot. It might simply be that the receiver uses the clock it has, or derives its clock from the preamble. Then, it could just be on/off, or it could also be manchester encoded.


One common practice is to take the log of the squared magnitude of the DFT frames, as in log(re*re+im*im), scale that, and use the scaled log as an integer index into an RGB color lookup table ("hotter" colors for larger magnitudes).  The scale is determined so that the range in your data fits within the color index table.

A log scale is used, since it's a closer match to human perception of loudness than a linear scale.

A square root is un-needed, as that's the same as adjusting the scale after the log by 0.5.


Why is $y(n) = x(n) + x(n-1)$ a low pass filter?

What you have here is the equivalent of a moving-average filter. Specifically, it is a filter of order 1, whose impuse response is

$$h(n)=\delta(n)+\delta(n-1)$$

Taking its $Z$-transform, we get

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

There is a pole at $z=0$ and a zero at $z=-1$. Plotting the magnitude of the frequency response $H(\omega)\triangleq H(e^{-\imath \omega})=2\vert \cos(\omega/2)\vert$, you get the following curve

As you can see, this clearly is a low-pass filter. You can easily calculate the cut-off frequency from here on.


In process of image zooming, Sampling rate changes or stays constant?

Yes, it changes the sampling rate. Just like with audio, 1 kHz tone is 1 kHz tone after resampling from 48kHz to 96kHz. So the image has same content, but with more pixels.


What is the difference between pitch detection and Onset detection?

The two concepts are related to two different dimensions or aspects of music which might or might not be correlated.

Onset detection is concerned with finding the points in time at which sounds start. Doing this does not require prior knowledge of the particular pitch (or fundamental frequency) of the sound. It may indeed rely on the property that at the beginning of a sound, there is an increase of energy. Actually, you can very well perform onset detection on recordings which do not have a well-defined pitch (for example: drums, machine noises...).

You can think of pitch and onsets/rhythm as related to the "vertical" and "horizontal" organization of music. There exists types of material in which changes along these two dimensions coincide - so that a detection algorithm can take advantage of a change in one dimension to infer that a change has occurred on the other dimension. But this coincidence/correlation is contingent.


What Is the Difference Between MRF and Total Variation in Noise Removal?

These are two different concepts that you talk about. First, MRF gives you a framework to do discrete optimization of problems, which respect the Markovian property, that is a pixel is conditioned only on the neighboring ones (roughly stated). Typical applications include binary or multi-class labeling problems. Total variation on the other hand, is generally used as a regularization by adding the integral of the absolute gradient of the signal/image to the energy functional. This helps to neglect irrelevant detail and focus on important ones.

You cannot say one is better than the other, as they are not exactly contradictory things. It depends on the application and the energy function you use in the MRFs.

An Introduction to Total Variation for Image Analysis is a good total variation tutorial to begin with. Also, Chambolle provides a link between total variation and binary MRFs, leading to new algorithms in the work entitled Total Variation Minimization and a Class of Binary MRF Models:

  Total Variation Minimization and a Class of Binary MRF Models, Antonin Chambolle, 2005


Be careful with that version of the Cookbook.  I did not actually write it, although I gave permission to Doug to write it.  He has a few typographic errors.  The original has been moved.  I think this is maybe where it lives now but there appears to be a CR/LF problem with the file.

Anyway, to do trancendentals with a fixed-point processor is difficult but possible.  First you need a mathematical description of the transcendental function in terms of arithmetic operations that your fixed-point processor knows how to do.  These are usually an infinite series that is truncated to a finite length.  I have some expressions of these common transcendentals here.

The other thing you have to worry about is scaling of the coefficients to be integer values.  Do you know how to do that in fixed point?  That is fundamental.

Then finally the $\sinh(\cdot)$ is expressed in terms of the exponential function which is described in that Wikipedia article you cited.

Conceptually, doing fixed-point math is not too hard, but it is sorta a bitch.  You must worry about both precision and about overflow.


What is the meaning of Mutual Information beyond the numerical calculation?

Two additional intuitive takes on mutual information:

When two random variables are independent the joint distribution $p(x, y)$ and the product of the marginal distributions $p(x)$ and $p(y)$ are identical. One could thus assess the degree of independent between two random variables by computing a probabilistic distance between $p(x) \times p(y)$ and $p(x, y)$ - this distance being 0 when the two variables are independent. A common probabilistic distance between variable is the Kullback-Leibler divergence. If you take the Kullback-Leibler divergence between the joint distribution and the product of the marginals of two random variables, you end up with... mutual information.
From a compression / coding perspective, imagine you are given a sequence of $N$ pairs of observations $(x, y)$. You want to compress them into a file. Two strategies: storing all the (x) in one compressed file then independently all the (y) in another compressed file ; vs compressing the pairs. Using an optimal coder, the file size in the first case is $N \times H(X) + N \times H(Y)$, while in the second case the file size is $N \times H(X, Y)$. The second approach is more efficient if there's a relationship between the two observed variables! How many bits have we saved per observation? $\frac{N \times H(X) + N \times H(Y) - N \times H(X, Y)}{N} = I(X,Y)$! So mutual information tells us how many bit per observation do we save by coding two data streams jointly rather than independently.

I am not sure about your example, though... Mutual information is computed between two random variables (distributions). I can see how "book" can represent the distribution of words in a book ; but I am not sure what "word" means here. Mutual information also requires "paired" observations to be computed.


Thanks for the clarification.

If you're currently doing

plot(y)

then you should do

plot(x,y)

where

x = linspace(0,1000,length(y));


Time domain and frequency domain analysis equivalence

now, even a steady-state sinusoid can be thought of, in the limit, as a sum of weighted pulses, each with a beginning and some kinda end.  how do all of these "transient" signals add up to a steady-state sinusoid?  but they do:

$$\begin{align}
x(t) &= \lim_{T \to 0} \sum_{n=-\infty}^{+\infty} x(nT) \cdot \operatorname{rect}\left( \frac{t-nT}{T} \right) \\
 &= \lim_{T \to 0} \sum_{n=-\infty}^{+\infty} x(nT) \cdot \frac{1}{T}\operatorname{rect}\left( \frac{t-nT}{T} \right) \cdot T \\
 &= \int\limits_{-\infty}^{\infty} x(u) \cdot \delta(t-u) \cdot du
 \end{align}$$

where $ \operatorname{rect}(u) \triangleq \begin{cases} 1 \quad |u|<\tfrac12 \\ 0 \quad |u| > \tfrac12  \end{cases} $

The Fourier integral models the entire time function as an infinitely dense and infinitely large set of sinusoids.  in a similar manner, those steady-state functions can add to something with a transient.

$$ x(t) = \int\limits_{-\infty}^{+\infty} X(f) \, e^{j2\pi ft} \ df $$


Is there a way to obtain the impulse response of a discrete system by just knowing it's response to the discrete unit step function?

A simpler version of Phonon's answer is as follows.

Suppose that $y$ denotes the response of the system to the unit step function.
Then, as discussed in this answer,
in general, $y$ is the sum of scaled and time-delayed copies of the impulse
response, and in this particular case, no scaling is required; only time
delays.  Thus,
$$\begin{align}
y[0] &= h[0]\\
y[1] &= h[1] + h[0]\\
y[2] &= h[2] + h[1] + h[0]\\
y[3] &= h[3] + h[2] + h[1] + h[0]\\
\vdots~ &= ~\vdots
\end{align}$$
where each column on the right is a (unscaled and) time-delayed
impulse response.
Thus, we easily get that
$$\begin{align}
h[0] &= y[0]\\
h[1] &= y[1]-y[0]\\
h[2] &= y[2]-y[1]\\
\vdots~ &= ~\vdots\\
h[n] &~= y[n] - y[n-1]\\
\vdots~ &= ~\vdots
\end{align}$$
with nary a mention of filters, inverses, convolutions, integration, operators
and the like, just simple consequences of
the definition of linear time-invariant system.


best robust feature detectors for detecting similarities among archeological signs

There is a difference between edge detection and feature detection, which is usually a corner detection (e.g. Harris corner detector), blob detection or even region detection (e.g. MSER).

Once the features are detected in your image, they are described in some way. For point features (corners), there are many different descriptors.

I think the simplest yet powerful one is simply taking a patch around the detected feature location. The patch can be quite large and should be rotated in the direction of dominant orientation (for rotational invariance).

The feature vector can consist of just horizontal and vertical image gradients. For a 64 by 64 pixel patch, this would lead to 4096 dimensional vector. This is quite large, but can be effectively reduced by PCA (Principal Component Analysis). This is quite slow but provides good quality and robustness of the descriptor comparable to SIFT or SURF descriptors.

The SURF descriptor is based on Haar wavelet and is also quite easy to implement - it can produce smaller descriptor and hence PCA is not necessary.

You can take a look on this paper comparing various feature descriptors.


DC to phase relationship with Low pass filtered sine and square wave multiplication and it's pitfalls. IQ Demodulation

So if your sampling rate is four times the modulated signal frequency (and sampling synchronously) then $\omega= \pi/2$. So then
$$
\begin{array}
 II(t) &=& s(t) - s(t-2)\\
&=& A\sin(\pi/2 t + \theta) - A\sin(\pi/2 (t - 2) + \theta)\\
&=& 2A\sin(\pi/2 t + \theta)
\end{array}
$$
and
$$
\begin{array}
 QQ(t) &=& s(t-3) - s(t-1)\\
&=& A\sin(\pi/2 (t-3) + \theta) - A\sin(\pi/2 (t - 1) + \theta)\\
&=& A\sin(\pi/2 t + \theta - 3\pi/2) - A\sin(\pi/2 t + \theta - \pi/2)\\
&=& 2A \cos((\pi/2 t + \theta - 3\pi/2 + \pi/2 t + \theta - \pi/2)/2) \\
&&\sin((\pi/2 t + \theta - 3\pi/2 - (\pi/2 t + \theta - \pi/2))/2)\\
&=& 2A\cos(\pi/2 t + \theta)
\end{array}
$$

The trick will be to think about the square wave case, and see what this derivation tells you about that.


Estimation of Amplitude, Frequency and Phase of Linear Combination of Harmonic Signal Beyond the Leakage Resolution of DFT

Solving the Linear Combination of Real Harmonic Signal
The data model is given by:
$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \sin \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$
Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).

Estimating the parameters above requires 4 steps solution:

Estimate the model order ($ M $) / number of signals.
Estimate the frequencies accurately using direct method (Given the number of frequencies).
Estimate the amplitude and phase of the model given the frequencies (Usually using Non Linear Least Squares solver).
Estimate all amplitudes, frequencies and phase (Usually using Non Linear Least Squares solver).

The above works only in very very high SNR cases. For low SNR this problem is probably unsolvable (Certainly for large $ M $).
If the frequencies are well spread (Far from each other) then steps (1) and (2) can be achieved using DFT (fft()).
Yet in case the frequencies are close in the best case the estimation of the frequency won't be accurate (Leakage from near by frequencies will shift the peak) and in the worse the case frequencies will be blended and neither the number of signals nor the frequencies is achievable.

In the images above I used observation time of 12 [Sec] with 4 signals with frequencies of [0.95; 1.00; 1.05; 1.10] [Hz]. Since I took 1 / 10 of observation time in your question, it is equivalent of frequencies being 0.005 far apart in the 120 [Sec] observation time. I just wanted to deal with less samples.
As can be seen in the image above, in the DFT there is no way estimating the number of signals and certainly extracting the frequencies.
But using Super Resolution methods one could solve this case!
See some information about Frequency Domain Super Resolution in the section below.
Yet for this case I used different method for Super Resolution which I can't reveal (Used commercially by me as an advisor).
I also developed a method to estimate the number of signals.
In the case above indeed my method estimated 4 signals and the estimated frequencies were [0.9503, 1.0016, 1.0483, 1.0996] [Hz].
Then applying the next 2 steps I could even farther improve results:

As can be seen from above, the estimated signal is almost perfectly aligned with the ground truth (The model with known parameters with no noise).
Indeed the SNR here is very high yet still, without the Super Resolution method, it wouldn't work.
The full code is available on my StackExchange Signal Processing Q74024 GitHub Repository.
Unfortunately, I can't disclose the code of the estimation of the number of parameters and the Super Resolution.
Frequency Domain Super Resolution
Old non modern methods will use some Windowing tricks which falls short in most cases.
What you really need is Super Resolution. Super Resolution basically is equivalent to a "Window Method" with very narrow main lobe and almost no side lobes.
You can employ Compressed Sensing / Sparse Representation for Super Resolution in Frequency Domain.
One way to do so is solving the problem:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| F \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$
Super Resolution means, in that context, being able to resolve frequencies which are closer than what the observation time suggests (overcoming Leakage issues):

In the above you can see the DFT of a sum of 2 sines with the given frequencies. The Gaussian model is using $ {L}_{2} $ for regularization (Which is basically damped zero padding).
You may see that the $ {L}_{1} $ could resolve the 2 sines even when they are only 0.5 [Hz] apart with an observation windows of 1 [Sec].
The answer was taken from my answer for Super Resolution in Frequency Domain Using Compressed Sensing.


Matlab fast fourier transform (FFT) function and time in hours frequency : hertz or millihertz?

That function is not periodic. Using the Fourier transform to analyze its behaviour is questionable, to say the least. Maybe you're confusing things with the Laplace transform, which is more useful when analyzing control systems?

Why should the axis in your plot have Hertz as unit?

The discrete Fourier transform¹ does a base transform from a vector in $C^N$ to a vector in $C^N$; its physical significance is that the result is that each bin represents another frequency between $0$ and the the Nyquist rate (which, in the complex sampling case is the sample rate, and in the real sampling case is half the sample rate).

Since you only have a real function, I'd argue frequencies are relative to $\frac{f_{sample}}2$.

$f_{sample}$ is $\frac{4}{\text{hour}}=\frac4{60\cdot60\,\text{s}}= \frac1{900}\,\text{Hz}$ in your case. So each bin represents a frequency $$f\in [-\frac{f_{sample}}2 ; \frac{f_{sample}}2[\,=\,[-\frac{1}{1800}\,\text{Hz}; \frac{1}{1800}\,\text{Hz}[$$.

Now, frequency base vectors are uniformly spread over the nyquist rate, so the $n^\text{th}$ bin represents frequency $\frac{\frac{1}{900}\,\text{Hz}}{n_\text{bins}}$. $n_\text{bins}$ is the number of bins, which is the same as the input vector size; which is $4\cdot24=96$. So each bin is $\frac1{96\cdot900}\,\text{Hz}= 11.5741\,\mu\text{Hz}$. Not very handy.

I think you're somewhat confused what an FFT is. I think a textbook on signal theory will do you very well!

¹ you're using fft, which is a fast implementation of the DFT; not the continuous Fourier transform. That's an important difference.


What's the purpose of augmenting an image with a random "gradient"?

It's not generating a gradient in the computer vision / image processing sense but in the graphics community sense. For example, if you use Gimp's gradient tool you will fill images with values interpolated between 2 ends.

My guess regarding the augmentation is that the author is trying to learn invariance to some shading effects (eg., when the sun makes an angle with object in the scene and you can see walls or other surfaces becoming gradually darker).


Help me understanding the super equalizer source code

result[i] =
band_gain[0] * LPF(t, band_freq[0]) +
band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) +
band_gain[2] * (LPF(t, band_freq[2]) - LPF(a, band_freq[1])) +
band_gain[3] * (LPF(t, band_freq[3]) - LPF(a, band_freq[2])) +
band_gain[4] * (IMPULSE(t) - LPF(t, band_freq[3]));

is actually very simple to explain with pictures. First we assume f0 < f1 < f2 < f3.... lets look at one section

band_gain[0] * LPF(t, band_freq[0]) means create a low pass filter with cutoff frequency band_freq[0] it then says to multiply the amplitude of the filter by band_gain[0] something like this

now lets move onto the more complicated lines

band_gain[1] * (LPF(t, band_freq[1]) - LPF(a, band_freq[0])) means make a LPF with Cutoff band_freq[1] then subtract from it a LPF with cutoff band_freq[0] if you look at the image below thats the pink section. Then you multiply the result by `band_gain1'

you can find the other bands in a similar fashion . Now lets add them together, something like this, just 'superposition'

which makes our final filter look roughly like this

I hope those pictures make sense


I want to get impulse response of causal system

Several comments:

First, the usual symbol for the Z-transform is $z$ rather than $s$ (which is common for Laplace transform).

It is more common to express $H(z)$ in terms of $z^{-1}$. Note that the causal inverse transform for a single pole at $\alpha$: $H(z) = 1/(1-\alpha z^-1)$ is $h[n] = u[n] \alpha^n$.

So your transform $H(z) = \frac{3}{1 - 1/4 z^{-1} - 1/8 z^{-2}}$ can be expanded to the sum of two terms like the single pole above.


You can generate Hann window of size $N$ and use it to window your signal $x$ using the MATLAB command

w=hann(N);
x_windowed = x.*w;

Since you specified BW of 0.5Hz, you need normalize it with your sampling rate $F_s$. That is $w_s = \pi 0.5/F_s$. You need to choose $N$ such that your main lobe width of the Hann window falls at $w_s$. It is not an easy problem and you may have to do trial and error initially. You can use freqz(w) command in MATLAB to see the frequency plot of Hann window. You can hover an see at what normalized frequency your main lobe width is.


Well, without knowing the specifics of your algorithm, I assume that you basically have multiple functions (e.g. Gaussian or Lorentzian peaks, each with position and FWHM parameters, additionally maybe some polynomials for backgrounds etc.) and you sum them all up into one big "fit function" that you hand over to an optimization algorithm that jiggles the parameters around until you found a good fit in the sense of non-linear least squares.

What you can do then is computing the asymptotic standard errors. For this, you have to (most likely numerically) compute the Jacobi Matrix of your fit function (matrix of partial derivatives for each parameter (columns) at each measurement point (rows)). From this you can compute the variance-covariance matrix and then compute the errors for each parameter. This is an approach very often found in curve fitting software. In fact, there is a quite comprehensive explanation here at OriginLabs (see the section on Parameter Standard Errors). Also a more detailed overview can be found here at arxiv.

Other methods are discussed in this paper. They certainly are appropriate and maybe even mandatory in some cases, but I find them a bit too complicated for most "real world" applications.

If I were you, I would go with the asymptotic standard errors for a first start and see, if they fit the needs. It is an easy to implement strategy that is, as already said, also widely found in software packages.


Variations in MATLAB plot each time the program is executed

The plots are changing i.e. you are changing inputs to the model (be it random input sequence , random channel matrix , random gaussian noise ) .. there are lot of random combinations here.. which will definitely change the plot.. you cannot expect to have same plot everytime when you are dealing with random numbers..
You can control with using same seed in random number generator.

Use the same seed while dealing with random numbers

Quick google search on this gave me this which you can use

https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-repeatable.html


MIT exercise 6.003 HW2 - Concept of system initially at rest

The mistake is in the final integration of the term

$$\frac{1}{2i}\big(e^{it}-e^{-it}\big)u(t)=\sin(t)u(t)\tag{1}$$

because you seem to ignore the step function $u(t)$. Integrating $(1)$ gives

$$\int_{-\infty}^t\sin(\tau)u(\tau)d\tau=u(t)\int_0^t\sin(\tau)d\tau=\big[1-\cos(t)\big]u(t)\tag{2}$$


You can derive the expression for the coefficients by doing bilinear transformation of the following analog low-pass prototype filter

$$H(s) = \frac{w_0^2}{s^2 + (w_0/Q)s + w_0^2}$$

where $w_0$ is the cut-off frequency.

You can lookup the bilinear transformation on Wikipedia.

The filter used in the Android app is a Butterworth filter because the chosen value of Q is $1/\sqrt{2}$. Note that in the constructor the inverse of Q is computed and assigned to variable iQ which is used in the computation of the coefficients. Note also that the variable K holds the 'frequency-warped' value of the specified cut-off frequency. You can find more information about the frequency-warping phenomenon in the above link.

You can find many examples on digital filter design using bilinear transformation. I found this one, which is quite close to the Android example.


Wiener Filters - Obtaining / Estimating Image and Noise Spectrum

You are speaking about so called Wiener deconvolution, because Wiener filter is set by equation $ w = R^{-1}r $, where:

$ R = E[y(n)y(n)^H] $ - autocorrelation matrix of input signal
$ r = E[y(n)s(n)^*] $ - crosscorrelation vector of input signal and original signal.

So the practical solution when statistics $R$ and $r$ are unknown is the Method of Least Squares, where only data matrix $F$ and part of original signal $s$ are needed.

In your case you should have some a priori knowledge about formula entities. Try:

  https://www.umanitoba.ca/faculties/medicine/radiology/Dick_Gordon_papers/Dhawan%20et%CAal.%20%281985%29.pdf

or:

  http://www.eos.ubc.ca/research/cdsst/members/6_96.pdf

It's a pity I only concern this subject in radio applications so there is nothing more to advise.

Hope this helps!


Advantage of complex filtering in multirate applications

For one practical example, I'll point to GNURadio's frequency translating FIR filter block: https://github.com/gnuradio/gnuradio/blob/master/gr-filter/lib/freq_xlating_fir_filter_XXX_impl.cc.t

It's a channelizer  & decimator block that gains an efficiency by spinning a user provided FIR LPF up to be a complex BPF where the channel is and filtering first vs. spinning the channel down first.  The efficiency comes from performing the frequency rotation at the lower sample rate after filtering and decimation vs. before.


Eventhough you are asking for an answer about your problem with FFT usage, It seems to me, based on your definition of what you are looking for (detecting the presence of a carrier -a known signal?- inside a buffer of 4096 samples) that you would better use a "matched-filter" to detect the presence of that carrier, a known signal, inside the noisy 30ms measurement buffer.

Of course FFT would still help especially for high SNR cases, but for noisy data, matched filter is the optimal choice (under suitable conditions). And I have never heard of FFT by itself reducing noise. May be I am wrong.

To implement a matched-filter, you would simply convolve the time reversed carrier signal with the noisy observation and look for a thresholded peak.


Signal's Fundamental Frequency is different from Plotted Signal

With $k$ a continuous variable, the frequencies of the sinusoids in the first sum are $f_n=n$, i.e. they are integers. The frequencies in the second sum are integer multiples of $\frac32$. The highest fundamental frequency of which all those frequencies are integer multiples is $f_0=\frac12$. The frequencies in the first sum are given by $2mf_0$ (with integer $m$), and the frequencies in the second sum are $3mf_0$. With a fundamental frequency of $f_0=\frac12$, the fundamental period is $T=2$, as observed by you.


Indetermination in a complex fourier series

HINT:
You need to go back to the original definition of the Fourier coefficients via the integral, and treat $|n|=1$ as a special case. Your mistake happened there already (you divided by zero).


How can define an indicator which measures the degree of similarity between two signals?

The general topic of finding similarities between signals is wide ranging:

are the signals of same sampling, length, offset, shift or scale?
where do they take their values (discrete, real, complex)?
are they stationary? noisy?
what do you consider similar (whole signals, chunks, specific features)?
which are the invariances looked for?
and most important: what is your goal?

The following works provide an overview of some common similarity metrics, which one can normalize into an index:

Clustering of time series data --- a survey, section 2.2. Similarity/distance measures
Similarity Measures and Dimensionality Reduction Techniques for Time Series Data Mining, section 3. Similarity measures


It's whatever makes the most sense for the problem at hand.  Which can get confusing, and the same signal may have multiple definitions of "bandwidth" -- sometimes even in the same document, when, for example, one is reading up on or designing a communications system.
Useful Bandwidth
A signal's user is interested in the practically useful bandwidth of the signal, i.e., if it's analog can I fit comms-quality audio (300-3000Hz) on it, or hi-fidelity audio (20-20000Hz), or video (0-4.5MHz for old-style NTSC, and yes, that's debatable), and if it's digital what symbol rate could I achieve.
That 3dB point you're wondering about is often taken to be a signal's useful bandwidth.
Occupied (or assigned) Bandwidth
How much bandwidth the signal takes up in reality.  This is the amount of bandwidth that you need to assign for the signal in the real world, so that signal A doesn't interfere with signal B when you're using real-world equipment.  Typically occupied bandwidth is considered to be the bandwidth within which some large percentage (99%, 99.9%, etc.) of the total energy is confined, and in regulatory documents its usually defined in a way that makes it easy to test with reasonably-expensive instruments (i.e., a high-quality but not necessarily lab-grade spectrum analyzer).
Equivalent Noise Bandwidth
This is more a property of filters or systems, but the noise bandwidth of a filter is -- assuming the filter gain is normalized to 1 -- the ratio of the total noise power out of the filter to the total noise power in.  So, for example, the equivalent noise bandwidth of a 1$^{st}$-order lowpass filter with a 3dB point at $f_0$ is $\frac{\pi}{2}f_0$.
Whatever
Sometimes, in the context of your work, you'll come across some other bandwidth that you wish to define (such as null to null).  In that case, you state so clearly, and you use that as "bandwidth".  So if you run across some screwy bandwidth definition, accept that the author has their reasons, understand it, and roll with it.


What is the explanation for the inverse FFT of a cross-correlated signal being shifted by some amount in the time domain?

The issue is that the FFT version always assumes that the indices start numbering from $n=0$.  The result of the FFT version will be periodic, with period the same as the FFT length.

To show the same result as the correlation version, you need to use numpy.fft.fftshift.

Here is my attempt to generate your original signals:

Then the two versions of correlation:

which seems to do what you're seeing.

Whnn I apply fftshift to the FFT version, I get:

Code Only Below

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(-0.2, 0.5, 500, endpoint=False)
t2 = np.linspace(-0.2, 0.5, 250, endpoint=False)
t3 = np.linspace(-0.4, 1.0, 999)
square = (signal.square(2 * np.pi * t - np.pi/2 , duty=0.2) + 1)/3
triangle = np.zeros(500)
triangle[250:500] = (signal.sawtooth(2 * np.pi * t2 - np.pi, 0 ) + 1)/2

txcorr = signal.correlate(square, triangle, mode = 'full', method='direct')
txcorr2 = np.fft.ifft(np.fft.fft(square, 999)*np.conj(np.fft.fft(triangle, 999)))

plt.figure(0)
plt.plot(t, square)
plt.plot(t, triangle)
plt.figure(1)
plt.plot(t3, txcorr)
plt.plot(t3, txcorr2)
plt.figure(2)
plt.plot(t3, np.fft.fftshift(txcorr2), '.')
plt.plot(t3, txcorr)


Lloyd Max Quantization and Clustering : Part 1

Different.

Lloyd-Max is a special type of scalar quantizer design which is optimized (in terms of MSE) to source pdf. Hence the quantizer is generally non-uniform.
Lloyd's algorithm (and the more generalized LBG algorithm) is a scheme to design vector quantization. More information here and here.

That said, one should note that Lloyd's algorithm for VQ is essentially an extension of the scalar Lloyd-Max quantization to higher dimensions (vectorized sources).


Why do linear phase filters have symmetric impulse responses?

Actually, I think I see why.

$$X(j\Omega) = |X(j\Omega)|e^{-j\theta(\Omega)}$$

$|X(j\Omega)|$ is purely real, and therefore if we take the IFT it is even and symmetric.

$\theta(\Omega)= a\Omega$ since the phase is linear, so $e^{-ja\Omega}$ merely shifts the corresponding even and symmetric magnitude in the time domain, so the resulting impulse response is symmetric about $a$.


Notch filters: why is $H(e^{j0}) = H(e^{j\pi})$ required?

In general, there is no such requirement for notch filters that $H(e^{j0})=H(e^{j\pi})$ must be satisfied. You could definitely have a notch filter with $H(e^{j0})\neq H(e^{j\pi})$. Having the same gain at DC and at Nyquist is just a practical definition, and if you have a sufficient number of degrees of freedom (i.e., filter coefficients) you might as well make it a requirement. Ideally you want $|H(e^{j\omega})|=1$ for all frequencies except for the notch frequency $\omega_0$, for which you require $H(e^{j\omega_0})=0$.

As a simple example, take the most primitive notch filter, namely a FIR filter of length $3$. Its frequency response must be

$$H(e^{j\omega})=c\cdot \left(1-2z^{-1}\cos(\omega_0)+z^{-2}\right)\tag{1}$$

where $\omega_0$ is the notch frequency, and $c$ is a scaling constant. All degrees of freedom are taken by the requirement of the notch at $\omega_0$, so you don't get equal gains at DC and at Nyquist. Yet, it is a notch filter.


How to change the FFT phase domain of each frame by 90 degrees

It sounds like you're trying to plot the FFT such that the result is centered at 0, instead of at Fs. Matlab has a function for this fftshift.

Here's a simple example when trying to plot either the onesided or twosided FFT.

f = 100; % Signal frequency (Hz)
fs = 1024; % Sampling frequency (Hz)
N = 1024; % Number of samples
Nfft = 1024; % FFT size

xf = 20*log10(abs(fft(x, Nfft))); % Compute the PSD

fw1 = fs * (0:Nfft/2)./Nfft; % Onesided domain
fw2 = fs * (-Nfft/2:Nfft/2-1)./Nfft; % Twosided domain

figure
subplot(211), plot(fw1, xf(1:Nfft/2+1)), title('Onesided')
subplot(212), plot(fw2, fftshift(xf)), title('Twosided')


Fourier coefficients of sum of two functions with different fundamental periods?

Since this looks like a homework and solving people's homework is outside the scope of this group, here are some hints

To have a Fourier Series, a function needs to be periodic.
$$ g(t)=g(t+T_1)  \quad  \Rightarrow  \sum_k \left(a_k \sin(2\pi k \frac{t}{T_1}) + b_k \cos (2\pi k \frac{t}{T_1}) \right) $$
and
$$ h(t)=h(t+T_2)  \quad  \Rightarrow  \sum_k \left(c_k \sin(2\pi k \frac{t}{T_2}) + d_k \cos (2\pi k \frac{t}{T_2}) \right) $$
If  $z(t)=g(t)+h(t)$
Under what circumstances of $T_1$ and $T_2$ would
$$ z(t)=z(t+T_3) \quad \text{?}$$
an what circumstances of $T_1$ and $T_2$  would
$$
z(t)\ne z(t+T_3) \quad \text{for any} \quad T_3
$$

and if there were circumstance where $z(t)=z(t+T_3)$,

Could we re express $g(t)$ and $h(t)$ as Fourier Series periodic in $T_3$?


Here's just a polynomial version:

$$
\arcsin(x) = x + \frac{1}{2} \frac{x^3}{3} + \frac{1 \cdot 3}{2 \cdot 4} \frac{x^5}{5} + \frac{1\cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{x^7}{7}
$$

function y = arcsin_test3(x)
    y = x.*(1+x.*x.*(1/6+ x.*x.*(3/(2*4*5) + x.*x.*((1*3*5)/(2*4*6*7)))))
endfunction

which seems to have five multiplies (assuming you can save the result of x.*x) and three additions.

And the scilab plot is:

Top is scilab's asin vs this one, bottom is the error between the two.

Original Answer

The square root here might be a hassle, but I thought I'd write it up because it looks like fun. :-)

This page  suggests:

  from page 81 of the
  Handbook of Mathematical Functions, by Milton Abramowitz and Irene
  Stegun:

$$
   \arcsin(x) = \pi/2 - \sqrt{1 - x}(a_0 + a_1*x + a_2*x^2 + a_3*x^3),
$$
where
$$
   a_0 = 1.5707288\\
   a_1 = -0.2121144\\
   a_2 = 0.0742610\\
   a_3 = -0.0187293
$$

I've implemented this in scilab and it works OK, except around $x= -1$.  Just reflecting the $0 \le x \le 1$ over to $-1 \le x \le 0$ makes for a much better approximation.

The top plot shows scilab's asin function against the above approximation (in dashed red) against my change in green.

The bottom plot shows the error for my change (plotting that and the original on the same axes means the green looks zero everywhere).

// 25770
function y = arcsin_test(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-x).*(a0 + a1*x + a2.*x.*x + a3.*x.*x.*x)

endfunction

function y = arcsin_test2(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-xx).*(a0 + a1*xx + a2.*xx.*xx + a3.*xx.*xx.*xx)

    y = y.*sign(x);
endfunction

x = [-1: .0100001 : 1];

clf
subplot(211)
plot(x,arcsin_test2(x),'g.');
plot(x,arcsin_test(x),'r:');
plot(x,asin(x))
subplot(212)
//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')


How to Concatenate white noise signal and applying bandpass filter matlab

If I understand the question correctly this should solve your problem:

L = 5000;
P = 10000;
sigma = 2;
mu = 0;

left = sigma*randn(L,1) + mu;
right = sigma*randn(L,1) + mu;

y = [left zeros(L,1); zeros(P,2); zeros(L,1) right];    % stereo sound matrix


How to calculate arbitrary phase shift in discrete signal?

The answers to this question explain the solution for the continuous-time case. In discrete time, the solution is completely analogous. The impulse response of a phase shifter with phase shift $\theta$ is given by

$$h[n]=\cos(\theta)\delta[n]+\sin(\theta)g[n]\tag{1}$$

where $g[n]$ is the impulse response of an (ideal) discrete-time Hilbert transformer. Eq. $(1)$ shows that the solution is just a linear combination of the identity system and a Hilbert transformer (which performs a $90$ degree phase shift).


How to calculate SNR variation with period?

The FFT is equivalent to a bank of FIR band pass filters, in python you can compute the fft using numpy.fft.fft, and the central normalized (with respect to sampling rate) frequency of each filter is obtained with numpy.fft.fftfreq
For real signals you can use numpy.fft.rfft and numpy.fft.rfftfreq
If you want implementation guidance in stack overflow, give more details about the what you did, share some piece of code. You could share you slow running example and someone would help you to get it fast.


Suppose the transmitted codeword is $\mathbf C$. If no errors occurred, or one error occurred, then $\mathbf C$ is the unique codeword
that is at distance $1$ or less from the received word $\mathbf R$, and the decoding algorithm produces $\mathbf C$ as its output. When there are two errors, there is not a
unique codeword that is at distance $2$ or less from the received word. One
of the codewords at distance $2$ is $\mathbf C$, but there is at least one other
such codeword. So,
the decoding algorithm must be devised to handle this case. The possible options
are

Declare a decoding failure, that is, tell the world that an uncorrectable
error pattern has been detected. A less shameful alternative is used in ARQ
(automatic-repeat-request) systems in which the receiver sends back a
NACK signal to the transmitter which then retransmits $\mathbf C$. Thus,
the outside world is not told that the transmission is being repeated
over and over till the decoder is able to successfully decode, though the
outside world may be able to infer this because it has to wait extra long
for the next data bits to come out of the decoder.
Choose one of the codewords nearest to $\mathbf R$ and produce this as
the decoder output. In this instance, the decoding might be correct because
the decoder fortuitously chose $\mathbf C$, or it might be incorrect.

It is important to remember that decoding algorithms can output a
wrong codeword even when they are fairly certain they got it right.
If $\mathbf C$ is transmitted and the channel error pattern $\mathbf E$
is such that $\mathbf C + \mathbf E = \mathbf R$ happens to be another
codeword $\hat{\mathbf C}$, then the decoder output will be
$\hat{\mathbf C}$ and the decoder might even declare that it is
confident that it has found the right codeword because it did
not find any errors in $\mathbf R$; it matched $\hat{\mathbf C}$
in every bit! Pride, as we are told, goeth before a fall.

Finally, turning to the OP's specific question about the
$(8,4)$ extended Hamming code, the decoder can
tell that an even number of errors have occurred if the
XOR sum of the $8$ received bits is $0$ and the rest of the
syndrome (the part that is exactly the same as in the
$(7,4)$ Hamming code)  is nonzero. So, whenever this
happens, the decoder can declare a decoding failure
or send a NACK to the transmitter, or output a codeword.
An easy way out in this last case
is to assume that the overall parity bit
is in error (and thus flip it) and then use the rest of
the nonzero syndrome to correct one error in the other
$7$ bits exactly as if they were from a $(7,4)$ (non-extended)
Hamming code. This will give a codeword $\hat{\mathbf C}$
that is at distance $2$ from $\mathbf R$ and it will even
be the transmitted codeword $\mathbf C$ if there were exactly
two transmission errors, one in the $7$ bits and one in the
overall parity
bit.


Help in understanding from book expression of variance of an estimator : PRBS vs real valued input

A PRN sequence is a Pseudo-Random Noise sequence, often generated by using an Linear Feedback Shift Register (LFSR) with the feedback taps done by using a primitive irreducible polynomial in GF{2}, which is the Golois Field of 2 elements. When a primitive and irreducible polynomial in GF{2} is used, the LFSR will produce a "maximum length sequence", meaning the pattern at its output will be the longest possible for the number of registers in the LFSR.  An example generated from the polynomial $1+x^3+x^4$ is shown in the picture below. This polynomial is primitive and irreducible in GF{2}, and in this case we will use the symbols $-1$ and $+1$ as the two elements in the field), and therefore this LFSR as configured will produce a maximum length sequence; the state of the generator, being the values after each cycle at the output of each register in the generator, will cycle through every possible combination pseudo-randomly before repeating, except for the all zero state. (So in the picture shown the choices are $0001$, $0010$, $0011$, $0100$, $\ldots$, $1111$ in pseudo-random order, so the pattern of the sequence at the output will have $15$ values before repeating) Generally the maximum length is $x^N-1$ where $N$ is the order of the polynomial, the $-1$ is because the all-zero state is not allowed).

For clarity we will refer to each output of the LFSR as a "chip", and the complete sequence is referred to as a "symbol".

When you correlate to a matched PRN sequence + noise (the input), you multiply the input by each coefficient in the sequence and sum the result (correlation = multiply + accumulate).  If you are aligned to the code, the multiplication will cause the $-1$'s in the input to be $+1$'s, and the $+1$'s in the input to be $+1$'s, essentially converting the pseudo-random pattern to be all $1$'s. This output of all ones is then accumulated to achieve the correlator output.  (Below shows the multiplication that would then be fed to an accumulator after each result).

The accumulation on the signal component (which is all $1$'s when aligned) will go up in magnitude by $N$ given $N$ samples in the sequence (if we assume one sample per chip).

Assuming the input noise itself is uncorrelated to the PRN sequence, and white within our sample space (which means that each sample of noise is uncorrelated to all the other samples), then the multiplication by $+1/-1$ will not change the variance of the noise. Therefore the standard deviation of the noise will increase by $\sqrt{N}$ in the accumulator. (When you add N uncorrelated noise samples, the standard deviation goes up by the $\sqrt{N}$).

So the variance of the noise is increasing by $N$, but the variance of the signal is increasing by $N^2$ (since the signal magnitude, not power, went up by a factor of $N$). The end result is the signal-to-noise ratio (SNR) changes by a factor of N in power or $10\log(N)$ in $\textrm{dB}$, and this is referred to as the "processing gain" of the PRN correlator.

If we normalize the signal level to be the same at input and output, the noise standard deviation will have changed by a factor of $1/\sqrt{N}$, and it's variance by $1/N$, matching your equation.

Stating this mathematically, consider our samples at the output of the multiplier as
$$X_i\quad\text{for}\quad i=1, 2, \ldots, N$$

The multiplication by $+1/-1$ does not change the distribution of the noise on the input samples (they will still be independent and identically distributed, IID), with the same variance as the input given as
$$\textrm{Var}(X_i) = \sigma^2$$

Multiplying by $\pm 1$ with a synchronized PRN causes the mean to be 1 as previously described,
$$E(X_i)=1$$

We will denote the output after accumulating over all $N$ samples as $Y$, but we will also divide by $N$ to be consistent with the normalized autocorrelation $r_{uu}=1 = E(Y)$:
$$r_{uu} = \frac{1}{N}\Sigma_1^N X_i = E(Y)$$

Where $X_i$ is the product of the input waveform with the PRN.

This results in the expression for the mean and variance of $Y$ as follows:

\begin{align}
E(Y) &= E\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N}E\left[(X_1+X_2+ \ldots +X_N)\right]\\
\textrm{Var}(Y) &= \textrm{Var}\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N^2}\textrm{Var}\left(X_1+X_2+ \ldots+X_N\right)
\end{align}

\begin{align}
E\left[(X_1+X_2+ \ldots +X_N)\right] &= E[X_1]+E[X_2]+ \ldots E[X_N] = N\\
\textrm{Var}(X_1+X_2+ \ldots +X_N) &= \textrm{Var}(X_1)+\textrm{Var}(X_2) + \ldots \textrm{Var}(X_N) = N\sigma^2
\end{align}

Therefore
\begin{align}
E(Y) &= NE[(X_1+X_2+ \ldots +X_N)] = 1\quad\text{and}\\
\textrm{Var}(Y) &= \frac{1}{N^2}\textrm{Var}(X_1+X_2+ \ldots +X_N) = \frac{1}{N^2}N\sigma^2 = \frac{\sigma^2}{N}
\end{align}


A wild guess. $L$ denotes a region $R$ around pixel $g(x,y)$, potentially a symmetric one of size $(2L+1)\times (2L+1)$. $m_L$ and  $\sigma_L^2$ are the average and variance in $R$. And $\sigma_n^2$? Probably the estimated noise variance over the whole image, supposed stationary.

In a flat region with noise $\sigma_L^2 \sim \sigma_n^2$, the result could be close to the average. Around edges, one expect an higher response, so sharper. Potential problems arise when  $\sigma_L \ll \sigma_n $. Another interpretation:

$$\hat{f}(x,y)=(1-\frac{\sigma_n^2}{\sigma_L^2})g(x,y)+ \frac{\sigma_n^2}{\sigma_L^2}m_L(x,y)$$
is a sort of weighted mean (since weights sum to one) of the image and its smoothed version, whose dependency on $(x,y)$ is made more explicit.


Estimating a Low Frequency Signal Corrupted by High Frequency Noise

To answer your final question:

  Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.

That's because you're feeding them the wrong signal model.

TL;DR: Use the right tool for the job!

Gory Details

Any time you start in with the Kalman Filter, you are assuming that the signal (measurement in your equations above) was generated as:
$$
x[k+1] = \mathbf{A} x[k] + u[k]\\
\mbox{measurement}[k] = \mathbf{H} x[k] + v[k]
$$
where
$$
\mathbf{A} = \left [ \begin{array}{cc}
1 & 1 \\
0 & 1
\end{array} \right]\\
\mathbf{H} = \left [ \begin{array}{cc}
1 & 0
\end{array} \right]\\
$$
and $u$ and $v$ are zero mean, Gaussian distributed, independent random variables with covariances
$$
\mathbf{Q} = \left [ \begin{array}{cc}
1.6 & 2.4 \\
2.4 & 4.8
\end{array} \right] \times 10^{-9}\\
\mathbf{R} = 2.35 \times 10^{-6}\\
$$
At least for the equations you cite. More general application of the KF is possible, but usually not done.

So, one thing that might be done to improve your estimates is to ensure that $\mathbf{Q}$ and $\mathbf{R}$ are closer to reality.  From the look of things, your process noise ($u[k]$) has about the same variance as your measurement noise ($v[k]$).  Make them closer.

More importantly, your signal is correlated with your noise because, rather than stochastic noise, you have harmonic noise. Even worse: you have synchronized harmonic noise (meaning the noise moves in lockstep with your signal).

So: the problem is that by using the Kalman Filter your signal model is all wrong.

A better model for the signal is in the frequency domain:  the signal is low frequency. The noise is high frequency.  So, a low pass filter will do a much better job of improving things that the Kalman Filter can hope to do.

If I apply the Kalman filter to your problem (with some modifications to your equations to take account of the difference between $\mathbf{X}[k|k-1]$ and  $\mathbf{X}[k|k]$) then I get the picture below.

The black curve with dots is the measurement.  The red curve is signal.  The blue curve is the output of the Kalman filter.  The green curve is the output of a 5 point moving average filter (a simple low pass filter).

The sum of squared errors for the Kalman filter is 118.63722. The same figure for the low pass filter is 14.18046.  Clearly, the signal model of a low pass filter fits better because the error is smaller.

R code to implement this is below.

# 26489
T <- 128;
t <- 0:(T-1)/T*70

signal <- 4.4 + sin(0.06*pi*t - 0.1) + sin(0.14*pi*t + 0.07) + sin(0.2*pi*t + 0.4)
noise <- sin(0.06 + 0.5*pi*t) + sin(0.1 - 0.68*pi*t) - sin(0.04 - 0.74*pi*t) + sin(0.03 - 0.93*pi*t)

measurement <- signal + noise

xkm1km1 <- matrix(c(4.4, 0),2,1)
Pkm1km1 <- matrix(c(1,0,0,1),2,2)
H <- matrix(c(1,0),1,2)
A <- matrix(c(1,1,0,1),2,2)
Q <- 10^-6 * matrix(c(1.6,2.4,2.4,4.8),2,2)
R <- 10^-6 * matrix(c(2.35),1,1)

library("MASS") # For pseudo inverse ginv()

zhat <- t*0

for (k in 1:T)
{
  xkkm1 <- A %*% xkm1km1
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  z <- matrix(c(measurement[k]), 1, 1)
  xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
  Pkm1km1 <- (matrix(c(1,0,0,1),2,2) - K %*% H) %*% Pkkm1
  zhat[k] <- as.numeric(H %*% xkkm1)
}

plot(t, measurement, type="b", pch=19)
lines(t, signal, col="red",  lwd=10)
lines(t, zhat, col="blue", lwd=5)

lpf <-  filter(measurement, c(0.2,0.2,0.2,0.2,0.2), circular = TRUE)
lines(t,lpf, col="green", lwd=4)

errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
print(errs)


Algorithms for detecting the difference between two signals

This is a pretty common problem, and you’ll see it all the time in radar/sonar/___ar processing. Cross-correlation is really the thing you’d want to use here in my opinion.

In layman’s terms, you can use cross-correlation to do “template matching”. In radar applications, were looking for the signal we sent out to return to us. Since we know what signal we sent out, we can use that as a template and use cross-correlation as a mathematical tool to “search” for that template signal in the received signal.

So in your problem statement, you have the signal X which and the signal Y. It seems that Y = X + N, where N is some noise signal. To use cross correlation in MATLAB, you can simply just do [xc, lags]=xcorr(Y, X). Now ideally, the peak of the magnitude of xc will tell you precisely where the signal X “starts” in your signal Y. It’s important to note that the output of the MATLAB xcorr function will include “negative” time lags, so when you plot it, use plot(lags, abs(xc)).

If you plotted the magnitude of xc, that is abs(xc), you’ll likely notice a discernible peak above all of the noise. This would tell you that indeed, the signal X is present in the signal Y.

Now that’s just a simple overview of cross-correlation; naturally there are a lot of fine details, but hopefully that should get you started!

Implementation note for cross-correlation:
If you wanted to do circular cross-correlation or not use any built ins, you could do it with discrete Fourier transforms (DFT/FFT):
xc = IFFT( FFT(Y) .* conj(FFT(X)) ). If you’d like to learn more about the mathematics of this, check out Wikipedia and the convolution theorem.


You'll need to plot your signal's imaginary and real (I and Q) part at the exact symbol time – in other words, you'll need timing recovery to get the constellation from a time signal.

Timing recovery is among the more complicated things you'll need to implement for a receiver, depending on the actual kind of transmission you're dealing with.

Assuming you weren't aware of that, you might want to  look into these terms (I think the order makes sense, didactically):

pulse shaping filter
eye diagram
Nyquist ISI criteria
timing recovery

square timing recovery
Gardner timing recovery

phase error correction

Costas loop


Your signal is a superposition of fairly widely spaced impulses. Each impulse by itself has a flat magnitude spectrum and the intervals between impulses are only pseudo periodic so the phase response superpositions are effectively incoherent. The problem isn't your FFT, it is your signal model.

An actual ECG is much smoother. If you expect to see a predominant frequency at around 1Hz, the signal would look like a sine wave of around 1 Hz, not a sequence of isolated spikes.

You probably really want to measure the time intervals between spikes directly in the the time domain.


First of all, wavread has been removed in the 2015 & 2016 versions of MATLAB. So, if you plan on using your code in a newer version of MATLAB, I recommend that you use the audioread function instead.

Now, to answer your question. MATLAB already have a built-in function for plotting the power spectral density (Pxx) of a given input signal - periodogram. I've taken the liberty to copy a piece of sample code from MATLAB's documentation that should help you plot the Pxx you want (it will contain negative frequencies as well).

periodogram(x,[],NFFT,Fs,'centered')

All you need to keep from your code snippet are the input values seen here. The rest, like the axis function, aren't necessary to get the plot you want if you use this periodogram function.
The above function will plot the power spectral density of a given input signal x with no specified window ([] which just means the function defaults to a simple rectangular window). The final argument is the one that I presume you are most interested in: freqrange. It can have several different inputs, including twosided and centered. You should try both options and see which output you find more useful.

I hope this helps you!


Plotting frequency spectrum for a signal with 20MHz bandwidth and sampling frequency of 12.5GHz

Your fundamental limitation is the (comparatively) short window length of $M=400,000$ samples, at such a prohibitively high sampling frequency of $12.5$ GHz. This will roughly yield a few hundred kHz of effective FFT resolution per bin of the FFT spectral samples.

Using the formula for the spectral resolution of various windows we can deduce the following resolution, for example, for a standard Hamming window as

$$\Delta w = 8\pi / M$$ where $M$ is the window's length, that is $M=400,000$, which yields $$\Delta w = 8\pi / 400,000 = 6.2832 \times 10^{-5}$$ radians.

This resolution is determined by the window length, and its frequency equivalent under the sampling rate of $f_s = 12.5 \times 10^9$ Hz is

$$ \Delta f = \frac{\Delta w}{2\pi} f_s = \frac{6.2832 \times 10^{-5}}{2\pi} 12.5 \times 10^9 = 125,000 Hz$$ i.e., the effective frequency resolution of your window is $\Delta f = 125$ kHz. Now given your bandwdith of $B=20$ MHz, this equates to a number of $20 \times 10^6 / 125\times 10^3 = 160$ samples.

On the other hand the apparent count of the number of FFT samples (apparent resolution) inside the $20$ MHz region yields something about $k = 4\times 10^5 \times 12.5\times 10^9 / 20\times 10^6 = 640$ samples, which however does not carry more information than those inferred from $160$ subsamples.

Hence a bandwidth of about $20$ MHz would yield about less than 200 true-samples of a spectral region inside $400,000$ many apparent-FFT samples. This is quite narrow. You would either use a much longer observation interval or use a much reduced sampling frequency as @MBaz commented.


Every signal transmitted by an SDR is a quadrature signal: the sum of two amplitude-modulated signals, with carriers that are in quadrature (one is at 90 degrees from the other).
A bit of trigonometry shows that this is equivalent to an amplitude and phase modulated signal; in other words, a quadrature signal is equivalent to $A(t)\cos(2\pi f_c t + \phi(t))$.
As long as you can express your signal with one of these two (equivalent) expressions, you can do anything! You can certainly do FM on top of AM if you want to. For many of these "exotic" signals, the challenge will be in figuring out a demodulation process, but don't let that stop you from experimenting.


In signal averaging, will the noise be reduced if the noise is not Gaussian?

Consider two zero mean signals $x(n)$ and $y(n)$. Let $z(n)=x(n)+y(n)$. In the following $E()$ denotes the expectation.
$$E(z) = E(x+y)= E(x)+E(y)= 0$$
I've dropped the time index $n$ for convenience. So $z(n)$ has zero mean. Now let's look at the variance. Since $z(n)$ is zero mean - I'll just use the correlation.
$$E[z^2] =E[(x+y)(x+y)]=E[x^2+2xy+y^2]=E(x^2)+2E(xy)+E(y^2) $$

If the correlation between $x(n)$ and $y(n)$ is zero, then $E(xy)=0$, then the variances add. Usually this is done by stating the $x(n)$ and $(y)$ are independent which implies uncorrelated.

In the general case $E(xy)\neq 0$, so the resulting variance could increase or decrease depending on the sign of the correlation.

Note - stating that $x(n)$ and $y(n)$ are white can be a little unclear. Usually saying $x(n)$ is white implies it is uncorrelated with itself - it does not state how it is correlated to another noise process.

The above analysis is not dependent on the type of noise (Gaussian or not) and it does not make any assumptions on the autocorrelation functions of $x(n)$ and $y(n)$.


Determine spectrum amplitudes for half-wave rectified sine

You have to integrate by parts twice. Leaving out constants you have

$$\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt=-\frac{\cos(\omega_0t)}{\omega_0}e^{-jn\omega_0t}\bigg|_0^{T/2}-jn\int_0^{T/2}\cos(\omega_0t)e^{-jn\omega_0t}dt$$

If you again perform integration by parts on the integral on the right-hand side of the above equation you get the final result:

$$\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt=\frac{1+(-1)^n}{\omega_0}+n^2\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt$$

If $I$ is the integral on the left-hand side you have

$$I=\frac{1+(-1)^n}{\omega_0}+n^2I$$

and finally

$$I=\frac{1+(-1)^n}{\omega_0(1-n^2)}$$

Note that you have to consider the case $n=1$ separately. Also note that you could also compute the coefficients $c_n$ by direct integration using $\sin(\omega_0t)=(e^{j\omega_0t}-e^{-j\omega_0t})/2j$.


Whether Fourier transform formula be considered as Convolution or Correlation?

Correlation and convolution are basically the same operations. You can express the cross-correlation of two functions $f(t)$ and $g(t)$ by a convolution:

$$R_{fg}(\tau)=f(\tau)\star g^*(-\tau)$$

where $\star$ denotes convolution, and $*$ denotes complex conjugation.

If you evaluate the cross-correlation at $\tau=0$ you get the inner product of $f(t)$ and $g(t)$, and that's exactly what the Fourier transform is: it is the projection of $f(t)$ on the complex exponential $e^{j\omega t}$.

So, to answer your questions:

It is actually an inner product. Due to the equivalence of convolution and correlation, it can be seen as both, evaluated at $\tau=0$.
As already mentioned, the meaning of these terms is an inner product.


Relationship between frame rate and sampling frequency in speech signal processing?

The sample rate divided by the frame rate will give you the window shift in samples.  The ratio of the window (length - shift) to the window length will give you an overlap factor.  The sample rate divided by the window length in samples will give you the FFT result bin spacing.


Playing sound in left and right speaker alternatively and storing output in array in matlab

If you want to play only one channel then you must set zeros at one of your matrix columns, i.e.:

sound_L = sound_array;
sound_L(:,2) = 0;

sound_R = sound_array;
sound_R(:,1) = 0;

So your code becomes:

Fs = 44100;
N = 400
t = (1:N)*(1/Fs);
freq = 600;

sound_array = zeros(N,2);
%starts N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*freq*t);

%Left and Right channels

sound_L = sound_array;
sound_L(:,2) = 0;
sound_R = sound_array;
sound_R(:,1) = 0;

%Sound 1
sound(sound_L, Fs);
%Sound 2
sound(sound_R, Fs);
%Play Consecutive
sound(sound_array, Fs);


How to change the measured curve to be close to the target curve

It is unclear what is your application from your question.

If you aim at "gaussianizing" your data -- for example to make better use of statistical modeling, regression, or machine learning methods which work better under a gaussian assumption -- you can use the Box-Cox transform. It consist in applying the following transformation to your data:

$$y = \frac{x^\lambda - 1}{\lambda}$$

$\lambda$ is chosen as to maximize the log-likelihood of a gaussian model given the transformed data - this is a maximization problem usually solved using Levenberg–Marquardt.

I would also strongly suggest posting this to the stats SE.


Effect of zero padding an odd symmetric FIR filter in the time domain

The three options you mention all result in the same magnitude response. The ones with zeros before (to the left) of the actual impulse response will just add delay (as many samples as there are zeros). The last option has exactly the same frequency response as the original filter, without any extra delay.

I think it's pointless to add extra delay by prepending zeros. Zero padding at the end is OK, if you have a good reason for it.
If you really want a filter with $1200$ taps, why wouldn't you use the extra taps to get a better filter?

If you need a linear phase filter with an integer delay (in samples) then you could design a filter with $1199$ taps and append one zero tap at the end (if for whatever reason you need exactly $1200$ taps). This is because the delay in samples of a length $N$ linear phase filter equals $(N-1)/2$, which is only integer if $N$ is odd.


What is the DFT of $[x_1 -x_2 x_3 -x_4....x_n]$

Hint: Consider a factor 2 down-sampled sequence : $\tilde{x} = \{x_1, x_3, x_5, ...\}$ and consider another sequence which is circular shifted by $-1$ sample and downsampled by a factor of 2 : $\hat{x} = \{ x_2, x_4, x_6, ...\}$. Observe that the required sequence is difference of $\tilde{x}$ and $\hat{x}$. Each of the DFT can be computed by observation using properties.


Delay filters in Python using for loop and $\tt lfilter$

After a trial-and-error session I figured what was wrong in my code:

The coefficients I fed to lfilter had incorrect signs, because only
the signs of the coeffs from the a array have to be reversed and
those in the array b stays the same.
The backward iteration did not work properly for the IIR filter,
because when running backwards the delayed values were not yet
filtered and in that case equal to the input signal.

Here is the fixed code:

l = 1024
x = np.array(range(l))
w = x/l
y = np.sin(w*x)
d = 4

# Let's caclulate the signals filtered with IIR filter
## First using forward iteration
y11 = y.copy()
for i in range(2*d, l):
    y11[i] = y[i] - y[i-d] + 1.2*y11[i-d] - 0.8*y11[i-2*d]
## Second using backward iteration
y12 = y.copy()
for i in reversed(range(2*d, l)):
    y12[i] = y[i] - y[i-d] + 1.2*y12[i-d] - 0.8*y12[i-2*d]
## Third using lfilter
# was:
#y13 = lfilter([1, 0, 0, 0, 1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)
y13 = lfilter([1, 0, 0, 0, -1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)

# Let's caclulate the signals filtered with FIR filter
## First using forward iteration
y21 = y.copy()
for i in range(2*d, l):
    y21[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Second using backward iteration
y22 = y.copy()
for i in reversed(range(2*d, l)):
    y22[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Third using lfilter
# was:
#y23 = lfilter([1, 0, 0, 0, 2, 0, 0, 0, -1], [1], y)
y23 = lfilter([1, 0, 0, 0, -2, 0, 0, 0, 1], [1], y)

# Let's start the figure and plot the signal in the top row, then the filtered
# signals below it
fig = plt.figure()
ax1 = plt.subplot(421)
plt.title('Test signal')
plt.plot(x, y)
ax2 = plt.subplot(422)
plt.title('Test signal')
plt.plot(x, y)

# The left column is for the IIR-filtered signals
plt.subplot(423, sharex=ax1)
plt.title('Forward iter.')
plt.plot(x, y11)
plt.subplot(425, sharex=ax1)
plt.title('Backward iter.')
plt.plot(x, y12)
plt.subplot(427, sharex=ax1)
plt.title('lfilter')
plt.plot(x, y13)

# The right column is for the FIR-filtered signals
plt.subplot(424, sharex=ax2)
plt.title('Forward iter.')
plt.plot(x, y21)
plt.subplot(426, sharex=ax2)
plt.title('Backward iter.')
plt.plot(x, y22)
plt.subplot(428, sharex=ax2)
plt.title('lfilter')
plt.plot(x, y23)

plt.tight_layout()
plt.show()

Now on the plots it is clearly visible that iterations produce the same results as scipy's lfilter except for the backward iteration attempt at the IIR filter as it was explained above.


Assuming that in your question $T(r)$ denotes a map $:\mathbb{R}\to[0,1]$ and assuming that it is continuous and strictly increasing, then it is obviously an isomorphism among $r\in[0,1]$ onto its image.
As a consequence you get the following result: let $r$ be any value in $[0,1)$, and let $a=T(r)$, then for any $\rho\in[r,r+dr]$ it is $\alpha = T(\rho) \in [a,a+da]$, having denoted $a+da = T(r+dr)$. Therefore the probability of $r$ and $a$ to belong their infinitesimal intervals (which is actually the definition of PDF) is the same. hence as $dr\to 0$, also $da\to 0$ and you can write
$$
p_r(r) dr = p_a(a) da
$$
which yields the result


Least Squares Fitting to Inverse Exponential Function

I think you are going to be stuck with an iterative approach.

$$ y = a (1-e^{-bt + c}) + d $$

First rewrite your equation like this:

$$ y = g + h e^{-bt} $$

Where:

$$ g = a + d $$

$$ h = -a e^c $$

You have a set of values $(t_n, y_n)$.  Consider the $y$ values as a single vector $\vec Y $.  Then pick an initial value for $b$ and construct a vector $\vec B$ using your $t$ values in $ e^{-bt} $.  You now have a vector form of your problem:

$$ \vec Y  = g \vec U  + h \vec B $$

Where $\vec U$ is a vector of all ones.

You can solve for $g$ and $h$ by dotting $\vec Y$ with $\vec U$ and $\vec B$.  This gives you a linear system of two equations two unknowns which is easily solvable.

$$ (\vec Y \cdot \vec U)  = g (\vec U \cdot \vec U) + h (\vec B \cdot \vec U) $$

$$ (\vec Y \cdot \vec B)  = g (\vec U \cdot \vec B) + h (\vec B \cdot \vec B) $$

Note that $ \vec U \cdot \vec U $ is your sample count and $ \vec Y \cdot \vec U $ is the sum of your $y$ values.  Calculating $ \vec U \cdot \vec B $, $ \vec B \cdot \vec B $, and $ \vec Y \cdot \vec B $  can be done in single loop without actually having to store the vector values if you are memory constrained. Storing $ \vec B $ will save recalculation when you are measuring the fit.

To measure your fit, define the difference vector:

$$ \vec D  = \vec Y  - \left( g \vec U  + h \vec B \right) $$

Your square fit is now:

$$ f = \vec D \cdot \vec D $$

Your problem now is to find the value of $b$ which minimizes $f$.  There are several ways to do this.  The fastest converging one might be to pick three values of b, calculate an f for each, fit a parabola to your $(b,f)$ points, find the minimum $b_*$.  To iterate, choose three new tighter $b$'s around $b_*$.  Repeat until good enough.

The value of $g$ is the answer you are looking for.

Hope this helps,

Ced

Followup

As I mentioned in my comment to Ben, this equation is overparameterized.  My recommendation would be to get rid of $d$ by setting it to zero.  This leaves

$$ a = g $$

$$ c = \ln( -h/g ) $$

As a differential equation, per r b-j:

$$ y = g + h e^{-bt+c} $$

$$ y' = h e^{-bt+c} (-b) $$

$$ y' = b ( g - y ) $$

Getting a good initial value of b:

There's nothing like sleeping on a problem.  If your samples are evenly spaced in time then you can find a good initial value of b.

The simplest involves using three evenly spaced points at time $t_A$, $t_B$, and $t_C$.

$$ t_B = t_A + s $$
$$ t_C = t_B + s $$

Calculate the following ratio:

$$ r = \frac{ y_B - y_A }{ y_C - y_B } $$

The value of b can be now be solved for:

$$ r = \frac{ ( g + h e^{-bt_B} ) - ( g + h e^{-bt_A} ) }{ ( g + h e^{-bt_C} ) - ( g + h e^{-bt_B} ) } $$
$$ r = \frac{ h e^{-b( t_A + s ) } - h e^{-bt_A} }{ h e^{-b( t_B + s ) } - h e^{-bt_B}  } $$
$$ r = \frac{ e^{-b t_A }( e^{-bs} + 1 ) }{ e^{-b t_B }( e^{-bs} + 1 )  } $$
$$ r = e^{-b ( t_A - t_B ) } $$
$$ b = \frac{ \ln( r ) } { t_B - t_A } $$

In order to mitigate the effects of noise, more points can be used.  Another formula would be:

$$ r = \frac{ y_1 - y_2 + y_3 - y_4 + y_5 - y_6 }{ y_7 - y_8 + y_9 - y_{10} + y_{11} - y_{12} } $$

$$ b = \frac{ \ln( r ) } { t_7 - t_1 } $$

Any number of formulas can be created.  The rules are that the coefficients have to add up to zero and the time spacing has to be the same in the numerator and the denominator.  Also, you want to stay on the steeper part of the curve so the differences in values are larger compared to your noise level.

With enough sample values used, the value of $b$ might be good enough that the first calculation of $g$ will be accurate enough that no iterations are needed.


Using FFT to demodulated a sine wave envelope with a cosine carrier

As you can see from your "Original FFT" plot, the real part of the FFT has even symmetry and the imaginary part has odd symmetry. When you zero out the lower sidebands, shift, and re-create the negative frequencies, you're forcing the FFT to have even symmetry (look at your "shifted FFT") plot. The signal you re-created corresponds to a cosine; this is what you get in the "Original vs Demodulated" plot, where the recovered signal is 90 degrees out of phase with the original.

The more common approach to demodulating in the frequency domain is as follows:

Zero out the negative frequencies. This is usually accomplished with the Hilbert transform. A signal $s(t)+j\hat s(t)$, where $\hat s(t)$ is the Hilbert transform of $s(t)$, is called an "analytic signal" and its spectrum only exists for positive frequencies.
Shift the spectrum back to baseband by multiplying the analytic signal with a complex exponential of frequency equal to $-f_c$.

The result of this process is the "complex envelope" which, despite its name, is real in the case of DSB-SC modulation of a single sinusoid.


Phase-coded modulation is another name for direct-sequence spread-spectrum (DS/SS) modulation. In a typical DS/SS communication scheme using binary phase modulation, a single bit is transmitted not as a sinusoidal pulse of duration $T$ and phase either $0$ or $\pi$ depending on whether the bit being transmitted is $0$ or $1$, but rather as a sequence of successive $N$ shorter sinusoidal pulses of duration $T/N$ with phases $0$ or $\pi$ as defined by a specific binary sequence (called a spreading sequence) $\mathbf b$ of length $N$ or its complement $\mathbf{\bar{b}}$ according as the data bit $d$ is $0$ or $1$. For example, with $N=7$, we would have a phase-coded transmission describable as
$$\begin{array}{ccccc}
d=0 &\leftrightarrow &\mathbf b = (1,1,1,0,1,0,0) &\leftrightarrow &(\pi,\pi,\pi,0,\pi,0,0)\\
d=1 &\leftrightarrow &\mathbf{\bar{b}} = (0,0,0,1,0,1,1) &\leftrightarrow &(0,0,0,\pi,0,\pi,\pi)
\end{array}$$
In the communication context, the interference between two different phase-coded transmissions using two different spreading sequences is determined by the cross-correlation between the spreading sequences. In the radar context, $\mathbf{\bar{b}}$ is not needed, but interference between two radar signals with different spreading sequences is still determined by the cross-correlation between the spreading sequences. The term spreading refers to the fact that the phase-coded signal has a bandwidth $N$ times larger than the non-phase-coded signal. Because the energy is spread over a larger bandwidth, the signal doesn't interfere as much with narrowband signals, and is not easy to detect even by broadband scanners. Thus, covert communication or covert scanning is easier, and signals are not that easy to jam, either by narrowband jamming or by broadband jamming, because the correlation receivers reject both kinds of jamming more easily.


How to calculate particles sizes for spheric particles with overlapping and superposition? (Example image included)

Here is what I experimented with:

Use ELSD to generate elliptic contours. You could basically use any edge detector, but since in the following stages I will benefit from circle detectors, it is good to already have some geometrical edges. Here is what the output looks like:

Collect all circular (also non-circular) edges via collecting the channels into a single image. This can also be done via simple thresholding. Again, the output looks like:

I will denote this image to be $C$.

We will now compute the circular hough transform on this image. However, there are many spheres (which appear to be circular) in varying radii. Thus, we select a range $\mathbf{r}=\{9, 12, 15..150\}$ to be the expected radii per hough transform. The pseudocode looks like:

   $c_x$ $\leftarrow$ $\emptyset$, $c_y$ $\leftarrow$ $\emptyset$

   for all ($r_i \in \mathbf{r}$) $\{$

      {x, y} = hough_circles($C$, $r_i$)

      $c_x \leftarrow \{c_x \cup x \}$

      $c_y \leftarrow \{c_y \cup y \}$

   $\}$

Depending on the peaks of the accumulator, Hough might report multiple center locations. A simple clustering might help avoiding this (Just average the close-by centers).

The final image with the resulting centers overlayed (in different colors), looks like:

You could customize the radius to be robust to different circles. Moreover, one can use elliptic Hough transform to have a better coverage of the particles (as some of them are not circular). Finally, the much-occluded large circles can be obtained after masking out the small ones - even with the help of a distance transform. If the false circle hypothesis are desired to be eliminated, you can always reduce the Hough sensitivity. Or you could always render the contours of the circle and compare to the detected edges.


How to determine if the system is invertible

It depends on what exactly you mean by "invertible". In system theory, what is often meant is if there is a causal and stable system that can invert a given system, because otherwise there might be an inverse system but you can't implement it.

For linear time-invariant systems there is a straightforward method, as mentioned in the comments by Robert Bristow-Johnson. First, note that the input-output relationship of an LTI system is given by the convolution integral

$$y(t)=\int_{-\infty}^{\infty}h(t-\tau)x(\tau)d\tau\tag{1}$$

Comparing $(1)$ to the system given in your question shows that the system is LTI with impulse response

$$h(t)=e^{-t}u(t)\tag{2}$$

where $u(t)$ is the unit step function. The corresponding transfer function is

$$H(s)=\frac{1}{1+s}\tag{3}$$

This system is causal and stable. However, its inverse system

$$G(s)=\frac{1}{H(s)}=1+s\tag{4}$$

is not stable. It has a pole at infinity. The input-output relation of the inverse system is

$$y(t)=x(t)+x'(t)\tag{5}$$

where $x'(t)$ is the derivative of $x(t)$. So in order to implement the inverse system you need a differentiator, which is not stable in the bounded-input-bounded-output (BIBO) sense. I.e., bounded input signals can lead to unbounded output signals (just use a rectangular input function as an example).


How to recognise partial text on mechanical indicators

If you are using an out-of-the box OCR engine, then there's not much you can do with significantly occluded letters.

However, if you are rolling your own OCR, say with some template-matching, then, there is nothing preventing you from learning/teaching partial letters and their valid combinations (like the 8-9 and 9-0 above). Just treat it as another letter in the set to be recognized.


what is the effect of variation of duty cycle of sampling frequency on a reconstructed signal

First establish the following convention:

$$DutyCycle = \frac{SampleClockOnTime}{SampleClockPeriod}$$

In the following discussion,

$$SampleClockOnTime = T_p$$
$$SampleClockPeriod = T_s$$

You can model the effects of $T_p$ on your recovered signal as a pulse amplitude modulation signal recovery problem.  The sample clock signal is a pulse with finite duty cycle.  The reconstruction signal starts by modulating the sample clock signal with the data values of your sampled signal.  Literally you can think of multiplying your samples values by a unit amplitude pulse (your sample clock) to create a physical time sequence.

In the frequency domain, you end up with a reconstructed signal that has altered frequency characteristics (roll off at the high frequencies) that are a direct consequence of the sample clock pulse width (aka duty cycle).

To complete the recovery of your signal, you will pass the modulated pulse signal through a low pass filter to obtain only the baseband components of the signal.  The resulting signal with have frequency magnitude characteristis of the form:

$$|F'(f)| = |F(f)||\frac{T_psin(pifT_p)}{T_s*f*T_p}$$

Where:
F'(f) is the recovered signal and
F(f) is the spectrum of the origianl sampled signal

The effects of the sample clock pulse width is represented by:
$$\frac{T_psin(pifT_p)}{T_s*f*T_p}$$

If you plot this, you will see the amplitude rolls off as frequency approaches the sample clock frequency and the roll off is becomes more pronounced as $T_p$ increases.

You can look for more information on this by searching on the topic of "Sample and Hold effect".  Most of the work I've seen on this topic is in the context of switched capacitor filters (which is how I was exposed to it).


Using oversampling to increase resolution of a DC-signal as input

My input signal is a dc signal (sensor output) and im getting kind of a headache to understand why oversampling can increase the resolution of dc signals.

The ADC puts out integers.  So, let $x$ be the integer that would come out of the ADC (100.3, say, or -333.3).  Now let $y = \lfloor r \rfloor$ be the quantization operation, where you get straight integers, like 100 or -334.
Now, let $r$ be the sum of some actual value $x$ and some nice well-behaved zero-mean noise -- $r = x + n$.  For the right sort of noise (zero-mean Gaussian that's big enough is one kind) the ADC output is a random variable with a mean equal to $x$.  So treat $y = x + n_{ADC}$, where $n_{ADC}$ is a zero-mean noise process.
Now when you average a bunch of samples of $y$, the effect of $n_{ADC}$ on the signal diminishes (in fact, its deviation goes down as $\sqrt n$ for an $n$-sample average).  That's why oversampling and averaging reduces the noise and increases resolution.

How is the low pass filtering achieved?

In the method you outline, you're averaging the input samples.  That's effectively low-pass filtering.
Note that the app note gets it slightly wrong: it's having you add up $2^n$ samples, then shift by $n$, effectively dividing by $2^n$.  However, by virtue of the averaging, you're effectively adding bits to the ADC resolution*.  Usually when I do this I just take the sum of the $2^n$ samples, knowing that there's a lot of noise in the LSBs.  If I do shift down, I shift down by $n/2$, not by $n$.

Is something like SNR a good specification to describe the quality of a dc signal?

Not really.  SNR is good when you have a well-defined signal strength.
For DC, it's much better to use $\pm x_e$, where $x_e$ is the anticipated error bars, or effective number of bits, or something.  There's plenty of examples out there -- I'd take my guidance from voltmeter specs, or from industrial instrument specs.
* You are adding to precision, not accuracy.  Accuracy is a different thing**.  The accuracy of an ADC is limited by things other than its noise characteristics.
** The best short description I have is throwing darts: accuracy is how well-centered the darts are around the bullseye -- precision is how tightly clumped they are.  If all your darts land on a point 3" over and 1" down from the bullseye, you're very precise -- but not terribly accurate.


This question seems to highlight some of the perils of translating concepts between the continuous domain and a discrete one.  My understanding is that bandwidth is defined in the continuous domain as the difference between the upper and lower band limits.  Since it is the interval length, whether the interval is open or closed doesn't matter.  Open or closed means if the endpoints are considered part of the interval.  Since we are always saying "up to the Nyquist", considering the endpoints to be open seems logical.

The extra twist in this is your stipulation that the DC component is zero.  This would seem to narrow the bandwidth at the lower end, but to what?  In the continuous case, it would be an infinitesimal amount.

Can this properly be translated to the discrete case you describe?  I don't think so.  There may be a common convention for doing so that I am unaware of.  Perhaps your approach, but then you should strike off the gap at both ends, since you have excluded  DC.

Upon some reflection, I am going to disagree with MBaz again, and say that $f_s/2$ is the theoretical accurate answer based on the continuous definition of bandwidth.

I will agree that these are purely theoretical considerations with negligible practical value, but then again, so are a lot of exam questions.

You should probably strike the word "radio", it stopped RB-J from giving a good answer.  Probably a better answer than what I just gave.

$$ H_k(\omega) = \sum_{n=0}^{N-1}W_N^{-kn}e^{-jn\omega}  $$

$$ W_N^{-kn} = \left( e^{-j \frac{2\pi}{N}} \right) ^{-kn} = e^{j kn \frac{2\pi}{N}} $$

$$
\begin{aligned}
H_k(\omega) &= \sum_{n=0}^{N-1}e^{j kn \frac{2\pi}{N}}e^{-jn\omega} \\
 &= \sum_{n=0}^{N-1} e^{-jn ( -k \frac{2\pi}{N} + \omega) } \\
 &= \sum_{n=0}^{N-1} \left( e^{-j ( -k \frac{2\pi}{N} + \omega )} \right)^n \\
 &= \frac{1 - \left( e^{-j ( -k \frac{2\pi}{N} + \omega )} \right)^N }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{1 - e^{-j ( -k 2\pi + \omega N )}  }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{1 - e^{-j \omega N }  }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{e^{-j \frac{\omega N}{2} }}{e^{-j ( \frac{-k\pi}{N} + \frac{\omega}{2} )}}
 \cdot \frac{ e^{j \frac{\omega N}{2} } - e^{-j \frac{\omega N}{2} }  }{ e^{j ( \frac{-k\pi}{N} + \frac{\omega}{2} )} - e^{-j ( \frac{-k\pi}{N} + \frac{\omega}{2} )} } \\
 &= e^{j \left( \frac{-k\pi}{N} - \frac{\omega (N-1)}{2} \right) }
 \cdot \frac{ \sin \left( \frac{\omega N}{2} \right) }{ \sin \left(  \frac{-k\pi}{N} + \frac{\omega}{2}  \right) } \\
\end{aligned}
$$

For similar math, see my blog article DFT Bin Value Formulas for Pure Complex Tones.  The section called "As N Gets Large" covers the relation between "alias" sinc, aka Dirichlet Kernel, and the normalized sinc function.


How to create a distortion compensation filter

The following addresses the CIC compensation and is not a general "distortion" technique.  But it is a straight-forward method to "compensate" non-varying "distortion".  If the frequency response is known the inverse of the frequency can be used to compensate.  Examples like the CIC filter, where a poor filter might be used because of reduced complexity, are compensated later the signal chain.  In this example the frequency response is know and the inverse can be used.  Note, with multi-rate filters you only want to use the "usable" spectrum after decimation.

Generally, to compensate a CIC filter the inverse of the CIC filters response
can be used to generate the compensation filter.  The CIC has
a response of (see reference [r2] or [r3])

$$
  H(\omega) = \left| \frac{sin(\omega D/2)}{sin(\omega M/2)} \right|^N
$$

Where D is the differentiate delay, M is the decimation rate, and N is the
filter order (number of cascaded filters).  The inverse can be specified as

$$
  H(\omega) = \left| \frac{sin(\omega M/2)}{sin(\omega D/2)} \right|^N
$$

Once we have the
frequency response of the compensation filter, we can simply choose the length
of FIR filter that we desire.  The length of the FIR is application
specific.  Obviously the longer the FIR filter the better compensation.

The following are plots of this straight forward compensation.

The following is the Python code to create the frequency responses and plots.

import numpy as np
from numpy import sin, abs, pi
import pylab

D = 1; M = 7; N = 3

Hfunc = lambda w : abs( (sin((w*M)/2)) / (sin((w*D)/2.)) )**N
HfuncC = lambda w : abs( (sin((w*D)/2.)) / (sin((w*M)/2.)) )**N

w = np.arange(1024) * pi/1024

G = (M*D)**N
H = np.array(map(Hfunc, w))
Hc = np.array(map(HfuncC, w))
# only use the inverse (compensation) roughly to the first null.
Hc[int(1024*pi/M/2):] = 1e-8
plot(w, 20*log10(H/G))
plot(w, 20*log10(Hc*G))
grid('on')

See [r1] for other approaches and $sinc^{-1}$ approximation.

[r1] Altera, "Understanding CIC compensation filters"

[r2] R. Lyons, "Understanding Digital Signal Processing",
    2nd ed., Prentice Hall, Upper Saddle River, New Jersey, 2004

[r3] R. Lyons, "Understanding Cascaded Integrator Comb Filters"