Why we need fourier transform of periodic signal although we have fourier series for periodic signal?

If you take the Fourier transform of a periodic signal, you'll arrive at the same result that you would if you had computed the signal's Fourier series. Its spectrum will consist of discrete lines that correspond directly with the terms in the Fourier series. As an example, note the following Fourier transform pair:

$$
\cos(\omega_0 t) \Leftrightarrow \pi \left(\delta(\omega - \omega_0) + \delta(\omega+\omega_0)\right)
$$

If you computed the complex Fourier series of $\cos(\omega_0 t)$, you would find exactly two nonzero coefficients: those that correspond to frequencies $\omega_0$ and $-\omega_0$.


If you had supersonic hearing, would you be able to hear the carrier wave of an AM radio station?

The AM or medium wave band extends from 526.5 kHz to 1606.5 kHz in Europe or 535 kHz to 1705 kHz in the US.

Ultrasound (not supersound! :-) frequencies range from 20kHz and up:

However, a bigger problem is that AM stations broadcast using electromagnetic radiation, whereas sound travels by acoustic waves in the air.


What is a good way to distinguish aircraft noise from other sounds?

You need to train a classifier. The way I would approach it is to allow users to submit audio recordings through your app, and indicate which sections of the recordings contain the plane. This would enable you to quickly crowd source a lot of data with which to build a model. Then you would find the STFT and volume of the signal over time. Finally, train the classifier (such as an SVM or neural network) using the amplitude spectrum and volume features. Once you have trained the classifier, you can use it to detect planes by extracting the same features on live audio (STFT isn't computationally expensive, so it should be feasible to do on the phone).

Here's a backgrounder on Audio signal classification.


On the Use of OMP Algorithm to Estimate Sparse Vector

Well, in your example, the channel isn't exactly sparse.

It has been shown that $\ell_0$ minimization can recover any $K$-sparse vector $x$ from observations $\Phi x$ as long as $2K < {\rm spark}(\Phi) \leq M+1$, when $\Phi$ is $M \times N$ (so that $x$ is $N\times 1$), i.e., $K<M/2$ more or less. This is a necessary condition which means that if $K$ exceeds $M/2$ it can be shown that there is no unique solution anymore (sidenote: for an MMV problem, i.e., if you have multiple snapshots, this condition can be relaxed to $K<M-1$).

The example you picked has $M=N=8$ and $K=4$. This means that you are operating right at the limit of $K=M/2$. To have an $\ell_0$ algorithm perform well, you would want to have $K$ a bit below this limit, even more so if you choose a heuristic one like OMP. Try $K=2$ or $3$ or increase $M$ and $N$ a bit, it should work better.

Further notes:

Results  like these are typically asymptotic in $M$ and $N$ and tend to become tight only as $M$ and $N$ are somewhat larger. Performance for small-scale systems is much harder to analyze.
OMP is a sub-optimal technique. It is quite sensitive to correlations in the measurement matrix $\Phi$. For $M=N$, your could use orth(randn(n,n)) to make $\Phi$ orthonormal. OMP would appreciate that.
As you asked whether there are conditions on the matrix $\Phi$. Yes, there are. It should satisfy the restricted isometry condition (RIC). If you draw it randomly, chances are high that it does, even more so the bigger it is (see first point).

*edit: Regarding your comment: yes, OMP is just not such a great algorithm. And your randn design is also not the greatest. Success rates for M=8, N=8 with your design (empirical): K=1: 100%, K=2: ~88%, K=3: ~65%, K=4: ~45%. If you use orth(randn(n,n)) instead, you should see 100% (though in this case you wouldn't need OMP). Going from 8 to 16 also improves things. The truly interesting case for sparse signal recovery is $M<N$. Here you can make $\Phi$ at least row-orthogonal to help OMP a little. Normalizing columns would be another good idea.

*edit2: To emphasize more clearly what I said above, here is a numerical example varying $K$ and showing the empirical success rate. Comparing pure random Gaussian design (red line) to one where the rows are made orthogonal to each other (blue line). It helps. Note that here $M=8 < N = 16$.


Help with Audio Signal noise removing using FFT and IFFT in Matlab

The short answer is do not null out the "mirror frequencies" that are located above $f_s/2$ that match the frequencies you want to keep. If your FFT was generated from a real signal, then when you do the IFFT you will get the real signal back as long as you did not zero out those upper frequencies (as you did).

The DFT (which the FFT computes) returns complex values but for real signals the “positive” frequency values will be the complex conjugate of the negative frequency values which together represent real signals. Due to how the FFT is computed, instead of returning the positive and negative frequencies from $-f_s/2$ to $+f_s/2$, instead those same frequencies are given from $0$ to $f_s$.

Consider Euler’s Identity:

$$2cos(\omega t) = e^{j\omega t}+ e^{-j\omega t}$$

If it wasn’t already clear, the general form of  $Ae^{j\theta}$ represents the magnitude and phase of a complex number as $A\angle{\theta}$ so above we see how the cosine is represented by two phasors each rotating in time in opposite directions (a positive and negative frequency) such that their sum is always on the real axis.

The FFT returns the frequency bins from 0 to one sample less than the sampling frequency: $n =0$ to $N-1$ bins where each bin is spaced by $f_s/N$ with $f_s$ as the sampling rate. Due to the cyclical nature of the FFT this represents the positive and negative frequencies by mapping everything above $f_s/2$ to the negative frequencies, specifically samples $N/2 + n$ for $n\ge N/2$ map to $-N/2+n$. The command fftshift in MATLAB or Octave does this specifically.

For this reason you should NOT null out all the "mirror" frequencies that are located above $N/2$ as without them your signal will be complex. Consider the example above with the cosine.

$$2cos(2\omega_n t) = e^{j2\omega_n t}+ e^{-j2\omega_n t}$$

If sampled such that there were 100 total bins, and $\omega_n$ was the spacing on one DFT bin, then the DFT would have non zero values at the 3rd bin (bin 0 is DC) and the 2nd to last DFT bin which after shifting would be at the expected +/-2 bins. If you removed the higher bin then the resulting waveform would simply be:

$$e^{j2\omega_n t}$$


n real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams?

You can't "agree" on levels, since the received amplitude depends on the attenuation on the channel, and that depends on distance and a lot of other factors.
So, your receiver depends on detecting something that it knows the right amplitude (or power) of – for example, a preamble, or an average squared symbol, or something like that.
How that is done depends on the system you look at – generally, this is part of the synchronization and equalization efforts every system has to make. And it really differs! For example, some continuously transmitting system will not insert any specific symbols for that, but depend on the receiver doing statistics. Others have a "silence" period followed by known symbols; other systems are packetized and depend on preambles with specific autocorrelation properties (but not necessarily 100% defined content). Other systems depend on a pilot tone, and there's many more options.


I will try to avoid math, because math and "how to do it" tutorials can be easily found.

So, I start by pointing out one VERY important thing: One does not compute Harris for a single pixel, but for a vicinity (a patch of image) around that pixel! Let $I(i)_{xx}, I(i)_{xy} ...$ be your derivatives for a point $i_0$, then,

$H = \left[
\begin{array}{cc}
\sum_{i\in V}I(i)_{xx} w (i-i_0) & \sum_{i\in V}I(i)_{xy}w (i-i_0) \\
\sum_{i\in V}I(i)_{xy} w (i-i_0)& \sum_{i\in V}I(i)_{yy} w (i-i_0)\\
\end{array}
  \right] $

The $w(t)$ is a Gaussian kernel.
The previous eq tells you to integrate the derivative values over vicinity $V$ around current pixel. Each value of the neighbors is multiplied with a value that shrinks as the distance increases. The law of decreasing follows a Gaussian, because $w(t)$ is Gaussian centred at $i_0$. And that's it with math.

Now, back to the empirical observations. If you use solely the derivatives, and that pixel is part of a linear structure (edge), then, you get a strong response for the derivatives.
On the other hand, if the pixel is at a corner (an intersection of two edges) then, the derivative responses will cancel themselves off.

Saying that, the the Hessian is able to capture the local structure in that vicinity without "cancelling" effect. BUT very important, you have to integrate in order to get a proper Hessian.

Having a Hessian, obtained using Harris method or by other means, one might want to extract information about the vicinity. There are methods to get numerical values on how likely is to have an edge at current pixel, a corner, etc. Check the corner detection theory.

Now, about "stable points" or salient points. Picture that you are in a foreign town with no GPS and only with a good map. If you are "teleported" in a middle of a street, you might locate the street on the map, but you cannot tell where exactly are you on that street or in what direction you should go to move left or right (wrt to the map).
Imagine now that you are at a intersection. Then, you precisely can point your position on the map!. (Of course, assume that two streets doesn't intersect more than once).

Imagine now that you must match two images. One acts as a map, and the other as the city. You must find pixels that can be uniquely described, so you can do the matching. Check images on this post for example of matching. These points are called salient points. Moreover, the corner points tend not to change their 'cornerness' properties when the image is scaled, translated, rotated, skewed, etc. (affine transforms) This is why they are called "stable".

Some points in the image allows you to uniquely identify them. These pixels are located at corners or at intersection of lines. Imagine that your vicinity $V$ is on a line. Except for the orientation of the line, you cannot find anything else from that vicinity. But if $V$ is on a corner, than, you can find out the directions of the lines that intersect, maybe the angle, etc.

Not all corner points are salient, but only corner points have great chances of being salient.

Hope it helps!

p.s. How to find if a point is corner or not, take a look at Harris paper.

p.p.s. More on matching, search for SIFT or SURF.

p.p.p.s. There is a "generalization" of the Harris method, called Structure Tensor. Check Knutsson seminal work!


Wavelet is a small wave, what is problem with using sin(x) of one cycle?

A 1-cycle sinewave wavelet corresponds to imaginary component basis vector for bin 1 of a rectangularly windowed DFT result, and thus has the same problems of the lowest bins near DC of a rectangular windowed FFT.

The frequency response of a rectangular windowed FFT is a Sinc function with an expanse of large ripples in the response. The main lobe of the frequency response would be very wide relative to the frequency of that wavelet's basis. And the complex conjugate image, or negative frequency response, is nearby to the lowest bin(s); and thus their side lobes can potentially cause interference with the positive frequency response, depending on phase.

e.g. the wavelet would be highly subject to interference  from noise and any other spectral content that wasn’t purely sinusoidal and exactly integer periodic relative to the wavelets length.


what is difference between structure image and texture image

There is a growing awareness that images are not homogeneous and that it could be useful to model them as a combinaison or superposition of morphological components that have distinct properties. One can cite four main components: piecewise-smooth parts, contours, geometrical/periodic textures and a last one, barely modelled, sometimes termed noise.

The edge-cartoon/texture is a common name given to a more macroscopic version with two components: structure and texture. Structure often refers to the piecewise-smooth parts and large-scale contours, texture to geometrical and almost periodic textures, noise. Sometimes, they are called geometry and oscillatory components. Related works are mentioned in Structure-Texture Image Decomposition: See for instance: Image decomposition toolbox

When the contours of an object are not distinct, its textural aspect can help to segment it. The color aspect is often put on the "texture" side.


How does a low-pass filter programmatically work?

You asked how a low pass filter works and mentioned that the filter uses past values of you're data.  This is a non-technical discussion of what happens in a low pass filter.

The low pass filter takes differing views (shifted in time) of your signal, scales them and adds them together.  You can imagine drawing your signal 3 times, one being current, the second being shifted by one sample time, the third being shifted by 2 sample times.

At low frequencies, all the views look very similar (shifting by a single sample barely changes where you are on the signal at any instant in time).  In this case, the three versions will add together in a constructive (or at least non destructive) way, so the signal passes through the filter.

Now moving to higher frequencies, each shifted version of the signal becomes more distinct at any given instant (sample point) and in fact may even reverse in sign.  At these higher frequencies, the three versions of your signal tend to cancel out (added destructively) so the signal becomes attenuated.

Different types of filters arrange for this kind of constructive/destructive interference to take place over appropriate frequency bands to create low-pass, band-pass or high-pass filters.


What you have there is "wrap around" and a bad case at that.
It's easy enough to fix though. Let's call your samples $x[n]$ and your max amplitude $x_{max}$, so your signal is currently bounded to $[-x_{max},+x_{max}]$
Go through the waveform sample by sample. Whenever you see a jump that's larger than your max amplitude, add or subtract twice the max amplitude to bring it back to where it was.
Something like
for (i = 1; i<n; i++)
{
   sampleDiff = x[i]-x[i-1];
   if (sampleDiff > xmax) { x[i] = x[i] - 2*xmax; }
   if (sampleDiff < -xmax) { x[i] = x[i] + 2*xmax;}
}

Make sure you are using a data type that accommodate a data range of at least  $[-2x_{max},+2x_{max}]$, which is quite likely NOT the original data type. Also make sure you don't start in a "wrapped" region, i.e. x[0] needs to be good


Is ICA appropriate for separating mixed signals when all source signals are NOT detectable at all sensors?

You should be fine, zeros in the mixing matrix are not a problem....and theoretically it should converge even quicker than if all sources existed in all sensors.


How to simulate analog VU meter response to audio?

if you have input signal as vector of floating point numbers, try following code. It won't behave as real VU meter, but by modifying two constants you should get something very similar. Visualize the value of C in logarithmical scale

// memory element (can be imagined as capacitor)
c = 0
kCharge = 0.1
kDischarge = 0.001

for each sample x:
{
  // rectify the signal
  x1 = abs(x)

  if x1 > c then
  {
    // charge
    c = c * (1-kCharge) + x1*kCharge
  } else
  {
    // discharge
    c = c * (1-kDischarge)
  }
}


How to design a filter with custom group delay -- why is it an IIR?

As pointed out by Peter K., it is true that many well-known techniques for designing FIR filters actually only design linear phase filters. However, FIR filters are very well suited for delay equalization, simply because the design process is much simpler than for IIR filters. The reason for this is the fact that the design problem can be formulated in such a way that it is still linear in the unknown filter coefficients, which is usually not the case for IIR filters. This type of design problem is known as complex approximation problem, because in the case of a desired non-linear phase response a complex frequency response must be approximated by the filter, unlike a real one in the case of a linear phase filter. You can find a lot about this topic, including Matlab programs, in this thesis. Have a look at the design example #2 on page 119 ff. where an FIR magnitude and phase equalizer is designed compensating the effects of an analog anti-aliasing filter. The following figure is taken from the thesis and it shows the group delay of the analog low pass filter, the group delay of the equalizer(s) (there are two with different magnitude responses), and the total group delay of the cascade of analog and digital filter, which is supposed to be constant in the passband:


Electret microphone RC filter design contradiction among information sources

First, the youtube source shows an 1st order active filter, maxim source shows 2nd order active filter and last source shows first order passive filter. That's three different designs.

A 1st order active filter can be implemented either in the negative feedback branch of the op amp (as in youtube source), this does not need to be grounded.

Or it can be done as a classical RC circuit followed by impedance converter or amplifier. Here, either the R goes to ground (highpass) or the C goes to ground (lowpass)

If you combine these two approaches you get a 2nd order filter. So, wether "the filter needs to be grounded", as you ask, depends strongly on where its implemented.


Kalman filter restricted by areal constraints

The simplest solution in terms of computational resources as much as simplicity, I guess, is the Unscented Kalman Filter (UKF). There, the clipping might easily be added to the calculation of the sigma-points, in a way that they remain in certain boundary. I recommend that you take a look at this article, for example (working in 2016.04.22, Constrained State Estimation Using the Unscented Kalman Filter, Rambabu Kandepu, Lars Imsland and Bjarne A. Foss)


How do i set the band filter such that it has a pass band of 100Hz to 300Mhz?

That's much easier in the frequency domain. You simply get
$$V_o(\omega) = V_i(\omega)\cdot \frac{1}{1+j\cdot \omega\cdot R_1 \cdot C_1}\cdot \frac{j\cdot \omega\cdot R_2 \cdot C_2}{1+j\cdot \omega\cdot R_2 \cdot C_2}$$

The first term is the lowpass and the second term is the highpass. To set the frequencies you need
$$R_1 \cdot C_1 = \frac{1}{2\cdot \pi \cdot 300MHz}$$
and
$$R_2 \cdot C_2 = \frac{1}{2\cdot \pi \cdot 100MHz}$$


Training a NN classifier on a single\double channel out of a surround sound dataset

Generally the order of a 5.1 channel signal is as follows: L (Left), R (Right), C (Center), LFE (Low frequency enhancement), Ls (Left surround), Rs (Right surround).
What you need is a downmix and there are many downmixing methods. Here's the recommendation of ITU:
$$
L' = L+0.7071C+0.7071Ls
$$$$
R' = R+0.7071C+0.7071Rs
$$
However, multichannel signal may contain some hidden information for your classification task. To get a better performance, I don't think downmixing to single channel is a good idea. Stereo signals may be a compromise between performance and training time.


Filtering Frequency Bands Out of a Signal

What would be the most simple approach to do this ?

I'm a GNU Radio person. So the simplest approach was this GNU Radio companion-designed flowgraph:

Generating this flowgraph file generates a python program that uses GNU Radio to do the signal processing, on as many CPU cores as there are.
You can then run it as
$ ./bandpass_filters.py -f {input file containing float32 samples, one after the other}

To decrease the complexity of all these bandpasses, I first used a low pass to reduce the sampling rate of the signal by half – that's Ok, since you are only interested in frequencies below half of half of your input sampling rate.
I tried this with a 200MB "dummy" file, creating eight 100MB output files – this, to and from a temp directory, took a whole 13 seconds.
I'm pretty hopeful it's fast enough! I benchmarked this a bit, and it seems the source is able to push through about 20 million samples per second; the slow part is the writing to the 8 files.
Note that this is a quick and easy solution – a proper, speed-optimized solution would probably use

decimation, because you're output sampling rate is still 110 Hz, which contains a maximum of 9 Hz (that's a waste of processing power, and hence, speed)
another pre-decimation by 2 through low-pass filtering to 27.5 Hz for the four filters below that rate

However, really, at your modest amount of samples (16GB = 4 Billion 32bit floats?) this isn't really necessary, if you wanted to grab a coffee, anyways.
That flowgraph was really nice and easy to design. Here I am, shamelessly plugging GNU Radio's Guided tutorials if you want to learn to design such signal processing flow graphs yourself; in fact, it's even more fun to use GNU Radio in live flowgraphs, i.e. in systems where the in- or output (or both) are physical devices (microphones, ecg sensors, sonars, radio frontends). Exactly the same flow graph can work with a sound card, if you replace the "file source" by an "audio source" (both come with GNU Radio) and specify a sound card-compatible sampling rate via the -r flag.
Now, since you probably know better than me how to work with Mathematica, here's what I'd recommend you implement:

(first, one time) You use the GNU Radio companion to generate a python program that does the signal processing for you.
You use Mathematica's BinaryWrite() to write your raw samples to a file as 32 bit floating point numbers, e.g. BinaryWrite("/tmp/samples.dat.f32", your_sample_vector, "Real32").
You use Mathematica's abilities to call executables on your machine to execute that python script. Now, I'm no Mathematica expert, so this is all Google-based-guessing: RunProcess({"/path/of/python/executable", "/path/of/generated/python/program", "-f", "/tmp/samples.dat.f32"})
You read (BinaryRead) the resulting /tmp/samples.dat.f32_1_3.dat (and …_4_7.dat etc) in Mathematica.

Optionally, you really consider doing less workflow in Mathematica and use the sheer mass of signal processing things that come with GNU Radio or its module ecosystem.


The samples at the input to the IFFT block $X_0$, $X_1$, ... represent the magnitude and phase for each output frequency represented by the index. The inverse FFT translates this block to the equivalent time domain signal which has real and imaginary components (complex time domain representing your baseband signal). The dual DAC combined with a quadrature mixer converts this to a real signal centered at $f_c$ whose spectrum is identical to the baseband signal.

By declaring our data at the input to represent frequency with magnitude and phase represented by the constellation is central to the OFDM concept. You could imagine (properly) that each individual input $X_n$ to be a QPSK modulated waveform (such as your picture suggests), which is then up-converted onto its own sub-carrier, and banded together with the other sub-carriers with their own modulation, side by side in frequency, which is exactly what is occurring. The use of the IFFT just does this efficiently by putting each QPSK modulation onto it's appropriate sub-carrier by assigning it to an input of the IFFT (which must have frequency inputs), and allowing the IFFT to compute what the composite time domain signal of all of those sub-carriers would represent.


With such short (known) stimuli, does it make sense to use generic analysis methods, or is it «better» (less work, results that are more usable) to tailor your processing to the transmitted processing (matched filter?) or apply a parametric model and some numerical optimization?


Entropy : do we prefer higher or lower entropy?

You seem to have a number of misunderstandings, which I'll try to clarify while also trying to help with your questions.

  The entropy of a source $H(S)$ gives the average codeword length to encode a given source alphabet. i.e. it is the average number of bits per symbol required to encode the information in the source.

While this is true, I think it's not the more enlightening way to look at the entropy. Think of the entropy as an inherent property of a source. The entropy measures the average information of each message produced by the source.

As an analogy, think of the source as a black box with a button; every time you press the button, the source tells you something. A high-entropy source will consistently (not necessarily always) produce messages that leave you saying "Wow! That's amazing". A low-entropy source will be boring most of the time.

Now, Shannon proved that it's possible to find an enconding from source messages to bits, so that the average number of bits required per message is equal to the entropy. This is why your statement above is technically correct. Note that, typical of Shannon, he didn't actually tell us how to design such an encoder (but he provided a few hints).

  Let $S_1=101101...$ be a message of length $N$ [...] The entropy $H(S_1)=\log(2)=0.69$ bits/symbol.

Now you can see this statement makes no sense. $S_1$ is a message, so it cannot have an entropy. What it does have is information. To find the entropy of the corresponding source, you need to look at its message set (alphabet) and its probability distribution. If the source has memory, you also have to take that into account.

  Is there a way to determine $n$, the minimum word length for a message sequence of $N$? For example, $S_1=(10)(11)(01)$ here $n = 2$.

I'm not sure what you're asking here, but I'll make a guess. One of the hints Shannon gave us about source encoding is that efficient encoding requires us to gather many source messages and encode them together. This is equivalent to defining a new source with a different alphabet (whose elements are concatenations of the original symbols) and a different probability distribution.

You seem to be asking about moving in the other direction: breaking up the source alphabet into simpler symbols. In the case of a source whose messages are strings of 0s and 1s, you can always define it as a source with alphabet $\lbrace 0,1 \rbrace$, but you have to be careful when doing so, as you may need to introduce memory to the source.

For example, say you have a source with alphabet $\lbrace 00, 01, 11 \rbrace$. To redefine it a source with alphabet $\lbrace 0,1 \rbrace$, you'll need to introduce a mechanism to avoid messages such as $0,0,1,0$, which the original source cannot produce.

  It is desired that the entropy of the system be maximized, so $H(S1)<H(S2)$, so should I prefer $k=8$ symbols then? Which entropy do we prefer - message with higher entropy and hence is this the reason why we prefer higher modulation?

This question does not make much sense to me, since a given source has a given entropy (remember that entropy is an inherent property of a source). So, what is there to maximize?

Now, if you're given a choice between two sources, one of which has higher entropy, what you should do is not clear either. What do you want to achieve?

From the point of view of a source coder (compressor), which is trying to compress the messages as much as possible, having a low-entropy source is preferable. For example, a PNG compressor which is fed single-color images will be able to compress them a lot.

From the point of view of the communications engineer, what you want to do is to match the information content of the source to the capacity of the channel.

Maybe what you mean is this: remember that Shannon says that we need to group source messages together to achieve larger encoding efficiency. This process will also result in a source of larger entropy -- but this means more bits per message, so how can this be better?!

To clarify this point, consider a binary source with alphabet $\lbrace 0,1 \rbrace$ and probability distribution $\lbrace 0.9,0.1 \rbrace$. Even though this source is low-entropy (most of its messages are boring), it's very hard to compress: you still need one bit per message!

Let's say that you redefine this source to have 1024 messages of 10 bits each, and calculate the new probability distribution. You'll find that a few of the new 10-bit messages are much more likely than the others. You can then use (for instance) the Huffman algorithm to find an encoding that approaches the entropy.

So, from this perspective, redefining your source to have larger entropy is better, because you'll be able to find better encoders.

About higher-order modulation -- that doesn't play a role in this problem at all.


Using superposition and time invariance to find RL circuit response

The answer to your last question is no. Nobody can explain how that answer was obtained because it's wrong. The fact that it is wrong is obvious, because due the term $e^{(t-2)}$ (with a positive power of $t$!), the output will keep increasing indefinitely, which is impossible for a stable system (such as the given RL-circuit).

Your answer is the correct one, and you must have made a mistake while checking the book's answer.


How do impulse response guitar amp simulators work?

When talking about modeling, there are two things that usually get modeled: 1. the guitar amp, and 2. the speaker cabinet. Only the latter is modeled by an impulse response, which means that the cabinet is simply represented by an LTI system and implemented by convolution. This is of course an approximation but it works fairly well. You can find a lot of such measured impulse responses on the internet.

As for the guitar amps, there is no such thing as an "impulse response guitar amp simulator". As mentioned in the comments and in Marcus Müller's answer, guitar amps are much too complicated to be modeled by a simple impulse response. They are non-linear and dynamical, and the corresponding digital models are quite complex. There's a good reason why the big players in amp modeling (Kemper, Fractal Audio Systems, Line6) do not publish their modeling algorithms.

There are two basic approaches to modeling of nonlinear analog circuits: 1. black-box models, which only try to emulate a measured input-output relationship using a set of test signals, and 2. white-box models, which try to simulate the actual analog circuit in all details.
You can read a good summary of the basics in A Review of Digital Techniques for Modeling Vacuum-Tube Guitar Amplifiers.


What is the need of IIR filter? When FIR have generally good/desirable characteristics?

The biggest advantage of IIR filter is that you can achieve the filter specifications with much lesser order compared to FIR.
In a recent question (Lowpass Hann filter in Python) a transition band $\omega_{cutoff-passband}=0.05\pi$ to $\omega_{stopband}=0.136\pi$ was achieved with just $N_{IIR}=4$, but using FIR filter with Kaiser window, the order $N_{FIR}=55$. This results in much less complexity in implementation (multipliers, adders) reducing area required for implementation on a chip. In the above exmaple, for a given output, IIR requires 4 multipliers and 5 addition operations. For FIR, 55 multiplications and 54 addition operations are needed. Of course, there are provisions for optimizing this, but you get the idea..
Also, for the given example, in the pass band the phase is nearly linear for the IIR filter, so it results in much less distorted output for the pass band frequencies (stop band frequencies will get distorted but they are attenuated away by the filter). In IIR, you have the additional flexibility of placing poles as per your requirement to fine tune the filter attenuation/ripple. So it may seems FIR has advantage in terms of pure linear phase but IIR isn't far behind.
Code used to compare the order
clc
clear all

fc1=1200;
fc2=3000;
fs=44100;
rp=1;
att=40;
[n,wp]=cheb1ord(fc1/(fs/2), fc2/(fs/2), rp, att);
[b,a]=cheby1(n,rp,wp,'low');
freqz(b,a)
[n1,Wn,beta,ftype] = kaiserord([fc1 fc2],[1 0],[0.9 0.01],fs);
hh = fir1(n1,Wn,ftype,kaiser(n1+1,beta),'noscale');
freqz(hh,1,1024,fs)

Zoomed view near transition band

As you can see Chebyshev Type-1 provides a much higher attenuation of frequencies beyond pass band cutoff ($\text{1200Hz}$).


Wifi uses OFDM, which uses symbols that can transport many bytes at once. That means you can't have arbitrary long packets, but always need to use the next multiple of a symbol duration.


Octave: How to attenuate all frequencies except one?

Think of the morse code as an OOK (on-off keying) modulated signal corrupted by additive noise. The optimal linear receiver that can maximize the SNR is matched filter. So "the best" filter you can choose is a filter matched to the frequency of $1$s. After doing that, you have a hypothesis test problem. So you need to choose a threshold (again, the optimal threshold is something known).

If you decided to go that way, I can give you the details.

Before any details, your code appears to be just a pseudo code. It is really not clear what specific waveform you assume. So I try to give you a general answer.

How to do it:

1 - Find the frequency of the $1$s, using FFT or any other scheme. Let's say it is $f_0$.

2 - You should know the bit duration (i.e. the time each $1$ or $0$ takes). Assume it is $T$.

3 - If the waveform of $1$'s is $s(t)$, then your matched filter is $h(t)=s(T-t)$.

4 - Calculate $y(t)=\int_{0}^{t}r(\tau)h(t-\tau)d\tau$, where $r(t)$ is your noisy signal at each interval $T$. This will give you the output of matched filter. We will look at this output at $t=T$ in each timeslot. In other words, you will have discrete values $y_k$ corresponding to each timeslot.

5 - Now you need a threshold $\gamma$ to decide which part of $y_k$ corresponds to a $0$ and which one is $1$. There are a bunch of options, depending on what you want to achieve, the statistics of noise, and the prior probabilities of $0$s and $1$s. If you are not so strict about optimality, assume equiprobable bits. In AWGN you can assume you will expect equally likely $0$s and $1$s which simplifies your choice. So an option would be just $\gamma$ equal to half of the magnitude difference between $0$s and $1$s.


How to determine sampling frequency for an x(t) signal avoiding aliasing?

The function you show has a frequency range of $-\infty$ to $\infty$ so it isn’t strictly band limited. To completely avoid aliasing, it needs infinite sampling.

Of course, you can use tables or Fourier identities to obtain the exact inverse.

Usually some acceptable criteria is specified for tolerable aliasing when specifying a sample rate.


Do filters work on a "per sample basis"?

Obviously, filters need more than one sample of input – I mean, a single sample is a single number, and how would a single number have something like a frequency? Filtering is something you apply to a digital signal, and signal is defined by being a changing entity – i.e. different samples.

Hence, digital filters always need sequences of samples. The typical (non-rate-changing) filter works in a matter that each time you "push in" a new sample, an output sample "pushes out" of the filter.

You should definitely brush up your knowledge on what a digital filter is, and especially how FIR filters work. The important point here is the understanding that application of a filter to a signal is a discrete convolution.

You're absolutely free to change the system with which you convolve any time – but you'd break the time-invariance that makes LTI (linear, time invariant) systems such as FIRs so "easy" to handle.

EDIT: with your clarification:

Yes, any FIR will need "history" of samples. But also, yes, you can recompute the coefficients anytime, and "switch over" to a different FIR. However, you must make sure your new FIR's "delay"/memory elements were already filled with the previous input samples, otherwise you'll get inconsistencies.


Step response given impulse response $h = e^{-|t|}$

HINT: In your work you implicitly assume that $t\ge 0$, so your solution is (probably) correct for $t\ge 0$, but you must also consider the case $t<0$.


Noise removal using adaptive noise cancellation algorithms in real time systems

In real time application adaptive noise cancellation uses single sensor for signal reception.The single sensor consist of a  delay  z-Δ to produce a delayed version of incoming signal d(n), denoted by x(n), which de-correlates the noise while leaving the target signal component correlated. So here not needed reference signal separately but delay creates uncorrelated noise from the incoming noise corrupted signal.


MUSIC algorithm for frequency and DOA estimation

The estimation of frequency for one dimensional signal (dimension is time) is same as one dimensional space data for narrow band signal. To make it clear, consider you have n time samples (Nyquist criteria is met) data. The problem of frequency estimation from this n sample is same as the following case. Consider you have a n element linear array (equal spacing) and you have received a time stamp data from a narrow band signal  i.e you have n samples collected at the same time from the n element array. Now once you demodulate or remove the carrier effect from the n samples then the direction finding problem is same as frequency estimation problem mentioned in case of time domain data.

The above is true when there are less then 'n' frequencies present in the signals. Why this is true: the phase difference of the received signal on the array directly related to the angle of arrival. Kindly recall that phase of a sinusoidal is constant.  So time domain data with "fixed sampling frequency" is same as samples on Linear array with equal spacing element(fixed time) after demodulation .

Please understand it is true for narrow-band signal only not broad-band chirp like signal.

Hope it helps.


Single-sided Z transform with difference equations and the system function

You misinterpreted the result of residuez. Note that there is also a fourth output argument, namely the multiplicity of the poles (which is $2$ in your example). In your case, the residuez function doesn't help much, because you know the poles already, and the result is already in a partial fractions form:
$$Y(z)=X(z)H(z)=\frac{5}{\left(1-\frac14 z^{-1}\right)^2}\tag{1}$$
You can either look up the inverse $\mathcal{Z}$-transform of $(1)$ in a table, or, maybe even more straightforward, compute the result by convolution in the time domain:
$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\tag{2}$$
Solving $(2)$ for the given sequences is very straightforward.


Since this question has multiple sub-questions in edits, comments on answers, etc., and these have not been addressed, here goes.

  Matched filters

Consider a finite-energy signal $s(t)$ that is the input to a (linear
time-invariant BIBO-stable) filter with impulse response $h(t)$, transfer function $H(f)$,
and produces the output
signal
$$y(\tau) = \int_{-\infty}^\infty s(\tau-t)h(t)\,\mathrm dt.\tag{1}$$
What choice of $h(t)$ will produce a maximum response at a given time
$t_0$? That is, we are looking for a filter such that the global maximum
of $y(\tau)$ occurs at $t_0$. This really is a very loosely phrased
(and really unanswerable) question because clearly the filter
with impulse response $2h(t)$ will have larger response than
the filter with impulse response $h(t)$, and so there is
no such thing as the filter that maximizes the response.
So, rather than compare apples and oranges, let us include the
constraint that we seek the filter that maximizes $y(t_0)$ subject
to the impulse response having a fixed energy, for example, subject to
$$\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt = \mathbb E
= \int_{-\infty}^\infty |s(t)|^2 \,\mathrm dt.\tag{2}$$

  Here onwards, "filter" shall mean a linear time-invariant filter
  whose impulse response satisfies (2).

The Cauchy-Schwarz inequality provides an answer to this question. We have
$$y(t_0) = \int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt
\leq \sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}
= \mathbb E$$
with equality occurring if $h(t) = \lambda s(t_0-t)$ with $\lambda > 0$
where from (2) we get that $\lambda = 1$, that
is, the filter with impulse response $h(t) = s(t_0-t)$ produces
the maximal response $y(t_0) = \mathbb E$ at the specified time $t_0$.
In the (non-stochastic) sense described above, this filter is
said to be

  the filter matched to $s(t)$ at time $t_0$ or
  the matched filter for $s(t)$ at time $t_0.$

There are several points worth noting about this result.

The output of the matched filter has a
unique global maximum value of $\mathbb E$ at $t_0$; for any other
$t$, we have $y(t) < y(t_0) = \mathbb E$.
The impulse response $s(t_0-t) = s(-(t-t_0))$
of the matched filter for time $t_0$ is just $s(t)$ "reversed in time"
and moved to the right by $t_0$.

a. If $s(t)$ has finite support, say, $[0,T]$, then the matched filter is
noncausal if $t_0 < T$.

b. The filter matched to $s(t)$ at time $t_1 > t_0$ is just the filter
matched at time $t_0$ with an additional delay of $t_1-t_0$. For this
reason, some people call the filter with impulse response $s(-t)$,
(that is, the filter matched to $s(t)$ at $t=0$) the matched filter for $s(t)$ with the
understanding that the exact time of match can be incorporated into
the discussion as and when needed. If $s(t) = 0$ for $t < 0$, then
the matched filter is noncausal. With this, we can rephrase 1. as
The matched filter for $s(t)$ produces a unique global maximum
value $y(0) = \mathbb E$ at time $t=0$. Furthermore,
$$y(t) = \int_{-\infty}^\infty s(t-\tau)s(-\tau)\,\mathrm d\tau
= \int_{-\infty}^\infty s(\tau-t)s(\tau)\,\mathrm d\tau = R_s(t)$$
is the autocorrelation function of the signal $s(t)$. It is
well-known, of course, that $R_s(t)$ is an even function of $t$
with a unique peak at the origin.  Note that the output of the
filter matched at time $t_0$ is just $R_s(t-t_0)$, the autocorrelation
function delayed to peak at time $t_0$.
No filter other than the
matched filter for time $t_0$ can produce an output as large
as $\mathbb E$ at $t_0$. However, for any $t_0$,
it is possible to find filters that
have outputs that exceed $R_s(t_0)$ at $t_0$. Note that $R_s(t_0) < \mathbb E$.
The transfer function of the matched filter is $H(f)=S^*(f)$, the
complex conjugate of the spectrum of $S(f)$.
Thus, $Y(f) = \mathfrak F[y(t)]= |S(f)|^2$.
Think of this result as follows. Since $x^2 > x$ for $x > 1$ and $x^2< x$ for
$0 < x < 1$, the matched filter has low gain at those frequencies where
$S(f)$ is small, and high gain at those frequencies where $S(f)$ is large.
Thus, the matched filter is reducing the weak spectral components
and enhancing the strong spectral components in $S(f)$. (It is also
doing phase compensation to adjust all the "sinusoids" so that
they all peak at $t=0$).

-------

But what about noise and SNR and stuff like that which is what the OP
was asking about?

If the signal $s(t)$ plus additive white Gaussian noise with
two-sided power spectral density $\frac{N_0}{2}$ is processed
through a filter with impulse response $h(t)$, then the output
noise process is a  zero-mean stationary Gaussian process with
autocorrelation function $\frac{N_0}{2}R_s(t)$. Thus, the
variance is
$$\sigma^2 = \frac{N_0}{2} R_s(0) = \frac{N_0}{2}\int_{-\infty}^{\infty} |h(t)|^2\,\mathrm dt.$$
It is important to note that the variance is the same regardless
of when we sample the filter output.  So, what choice of $h(t)$
will maximize the SNR $y(t_0)/\sigma$ at time $t_0$?  Well, from the
Cauchy-Schwarz inequality, we have
$$\text{SNR} = \frac{y(t_0)}{\sigma}
= \frac{\int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}
\leq \frac{\sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}} = \sqrt{\frac{2\mathbb E}{N_0}}$$
with equality exactly when $h(t) = s(t_0-t)$, the filter that is matched
to $s(t)$ at time $t_0$!! Note that $\sigma^2 = \mathbb EN_0/2$.
If we use this matched filter for our desired sample time, then at other
times $t_1$, the SNR will be
$y(t_1)/\sigma < y(t_0)/\sigma = \sqrt{\frac{2\mathbb E}{N_0}}$.  Could
another filter give a larger SNR at time $t_1$? Sure, because $\sigma$
is the same for all filters under consideration, and we have noted above that
it is possible to have a signal output larger than $y(t_1)$ at time
$t_1$ by use of a different non-matched filter.

In short,

"does the matched filter maximize the SNR only at the sampling
instant, or everywhere?" has the answer that the SNR is maximized only
at the sampling instant $t_0$. At other times, other filters could give
a larger SNR than what the matched filter is providing at time $t_1$,
but this still smaller than the SNR $\sqrt{\frac{2\mathbb E}{N_0}}$
that the matched filter is giving you at $t_0$, and if desired,
the matched filter could be redesigned to produce its peak at time
$t_1$ instead of $t_0$.
"why not make a filter that makes a really tall skinny spike at the point of decision. Wouldn't that make the SNR even better?"
The matched filter does produce a spike of sorts at the sampling time
but it is constrained by the shape of the autocorrelation function. Any
other filter that you can devise to produce a tall skinny (time-domain)
spike is not a matched filter and so will not give you the largest possible
SNR. Note that increasing the amplitude of the filter impulse response
(or using a time-varying filter that boosts the gain at the time
of sampling) does not change the SNR since both the signal and the noise standard deviation increase proportionately.
"The I&D will basically ramp up until it gets to the sampling time and the idea is that one samples at the peak of the I&D because at that point, the SNR is a maximum."
For NRZ data and rectangular pulses, the matched filter impulse response is
also a rectangular pulse. The integrate-and-dump circuit is a correlator
whose output equals the matched filter output only at the sampling instants, and not in-between.  See the figure below.

If you sample the correlator output at other times,
you get noise with smaller variance but you can't simply add up the samples
of I&D output taken at different times because the noise variables are highly correlated, and
the net variance works out to be much larger.  Nor should you expect to be able
to take multiple samples from the matched filter output and combine them
in any way to get a better SNR.  It doesn't work. What you have in effect
is a different filter, and you cannot do better than the (linear)
matched filter in Gaussian noise; no nonlinear processing  will give
a smaller error probability than the matched fiter.


Actually, if it's got 160 samples every 20ms, you get a sampling rate of 8000 (you're missing a zero), which is the standard for over-the-phone voice communication. from a Google search, I couldn't find any AMR-NB specification that uses a different sampling rate.

Also, I found this line here pretty useful:

  For AMR-NB, only an “8000” sample rate and one of the following
  bitrates are valid:
  4.75k, 5.15k, 5.9k, 6.7k, 7.4k, 7.95k, 10.2k or 12.2k

Wikipedia also mentions that AMR-NB uses a 8000 Hz sampling rate, you can be pretty sure that this is the case all across the board.


filter and resample or resample and smooth?

Could you describe what you want to find in your data -- peaks, trends, model ?
and what kind of noise there is -- $\infty$ possibilities ?
A plot would be good.
Without some more info, I can only say, with others, "try it and see".

As you may know,
a good way to smooth and resample (interpolate) in Python is with
scipy.ndimage.map_coordinates,
along the lines

import numpy as np
from scipy.ndimage import map_coordinates

def smooth_resample( y, newlen ):
    """ resample y to newlen, with B-spline smoothing """
    n = len(y)
    newgrid = np.linspace( 0, n - 1, newlen )
        # e.g. n = 3, newlen = 5 -> [0 .5 1 1.5 2]
    return map_coordinates( y, newgrid, mode="nearest", prefilter=False )

map_coordinates also works for 2d images, 3d voxels ... hence the name "ndimage".
For more, stackoverflow has quite a few Q+A s that use it.

(Cubic-spline interpolation uses the nearest 4 points to each x in newgrid,
whereas dsp smoothing filters are usually much longer.
For smoothing that preserves peaks,
see Savitzky-Golay filters, e.g.
Numerical Recipes p. 771 ff.)


As per the exact solution of a discrete LTI system:

$$y(t)=C A^k x_0+C\sum_{k=0}^{t-1}A^{t-k-1} B u(k)+D u(t)$$

and considering that $u(t)=\delta(t)$ and $x_0=0$ for the impulse response, we have:

$$h(t)=C A^{t-1} B+D \delta(t)$$

where $A$,$B$,$C$ and $D$ are the state space discrete matrices:

$$x(t+1)=Ax(t)+Bu(t)\\y(t)=Cx(t)+Du(t)$$

Which in this case are obtained as:

b=[1/6 -1/3 0];
a=[1 1/6 -1/3];
[A,B,C,D]=tf2ss(b,a);

Hence, aside numerical issues, both methods gave the same results:

[h,t] = impz([1 -2],[6 1 -2]);
stem(t,h);

t=(0:23)';
h(1)=D*1;
for k=2:length(t)
    h(k,1)=C*A^(k-2)*B;
end
stem(t,h);

So the $t$=23th term on the $h(t)$ signal is:
$$h(23)=C A^{t-1} B|_{t=23}=C A^{22} B=-3.39702260269318e^{-05} $$

Though small, the $t$th term never vanishes, hence, transfer functions with denominator different than the identity vector, will always be an IIR.


What happens if I choose the length of signal L > NFFT? and what's about choosing L different form NFFT?

Did you read the documentation? http://www.mathworks.com/help/techdoc/ref/fft.html

  Y = fft(X,n) returns the n-point DFT. fft(X) is equivalent to fft(X, n) where n is the size of X in the first nonsingleton dimension. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.

So it just ignores the signal after n samples.

  2 In the following code, I'm not sure if I used window correctly.

Yep.  You just multiply the window by the signal:

  But when I use window (hanning in the following code), I can't get the exact values of amplitudes?

Yes, It attenuates anything at the ends of the window, so the amplitude varies depending on when the window hits the signal.  If there's a huge peak in the signal, but it happens at the end of the window, it gets multiplied by 0.  If it happens in the middle of the window, it gets multiplied by 1.

For signals that are stationary across the window, then it doesn't matter when the window is applied to the signal, right? So you can just multiply by a constant (2 for a Hann window, I think, because the area under the window is 0.5) to get the true amplitude.

Also, if you use overlapped windows and certain windows (constant-overlap-add (COLA) constraint), I believe you can combine them to get the true amplitude for any waveform, even non-stationary ones.  https://ccrma.stanford.edu/~jos/sasp/Overlap_Add_Decomposition.html

  3 When L and NFFT get different values (I mean when I change L and NFFT,) then the values of amplitudes were different too. How can I get the exact value of amplitude of input signal?

The signal is stationary, so it shouldn't really matter how much of it you measure.  To get the true amplitude you probably need to divide by N, depending on matlab's implementation. https://gist.github.com/236567


How does sampling rate affect discrete filters?

IF we take the sampling period as T.First I will convert the filter to time domain impulse response as
$$
h(s)=\frac{1}{\tau s+1}\\h(s)=\frac{1/\tau}{s+\frac{1}{\tau}}
$$
this in time domain is
$$
h(t)=\frac{e^{-t/\tau}}{\tau}
$$
IF the above signal is sampled at time T the it can be written as
$$
h(Tn)=\frac{e^{\frac{-Tn}{\tau}}}{\tau}\\h(n)=\frac{e^{\frac{-Tn}{\tau}}}{\tau}
$$
if Z transform taken for the impulse response.Then the equation is
$$
H(z)=\frac{1/\tau}{1-e^{\frac{T}{\tau}}z^{-1}}
$$
Here we can see the sampling time coming in denominator which could put the pole inside or outside the unit circle


Deriving the Fourier transform of cosine and sine

Your work is OK except for the problem that the Fourier transform of
$\cos(2\pi f_0 t)$ does not exist in the usual sense of a function of $f$,
and we have to extend the notion to include what are called distributions,
or impulses, or Dirac deltas, or (as we engineers are wont to do, much
to the disgust of mathematicians) delta functions.  Read
about the conditions that must be satisfied in order for the Fourier
transform $X(f)$ of the signal $x(t)$ to exist (in the usual sense)
and you will see that $\cos(2\pi f_0 t)$ does not have a Fourier transform
in the usual sense.

Turning to your specific question, once you understand that impulses
are defined only in terms of how they behave as integrands in an
integral, that is, for $a < x_0 < b$,
$$\int_{a}^{b} \delta(x-x_0)g(x)\,\mathrm dx = g(x_0)$$
provided that $g(x)$ is continuous at $x_0$, then
it is easier to deduce the Fourier transform of
$$\cos(2\pi f_0 t)
= \left.\left.\frac{1}{2}\right[e^{j2\pi f_0 t} + e^{-j2\pi f_0 t}\right]$$
by musing on the fact that
$$\int_{-\infty}^\infty \delta(f-f_0)e^{j2\pi ft}\,\mathrm df = e^{j2\pi f_0t}$$
and so it must be that $\cos(2\pi f_0 t)$ is the inverse
Fourier transform of $\displaystyle \left.\left.\frac{1}{2}\right[\delta(f-f_0) + \delta(f+f_0)\right]$.


Calculating the Norm of Each Sample of the DFT Trasform

The simple answer: yes this can be done in various different ways:

Calculate the $a_{k}$ and $b_{k}$ coefficients of the regular
Fourier Series (http://en.wikipedia.org/wiki/Fourier_series) and then $a_{k}^{2}+b_{k}^{2}$ is the norm at frequency k. This is equivalent of splitting the original DFT equation into it's real and imaginary parts
Calculate the auto-correlation of $c_{xx}[n]$ of $x[n]$. Then calculate only the cosine coefficients $a_{k}$ of the Fourier Series of $c_{xx}[n]$. $a_{k}$ is the norm at frequency k.

However, each of this methods is computationally much less efficient than simply doing the DFT and then taking the magnitude squared. Using the DFT (with complex math) is by far the most efficient way of doing this. If you worried about computationally efficiency, you can shave about 1/3 to 1/2 of the computation time by using the fact that your input signal is real. The DFT of a real function can be calculated by using DFT of a complex signal of half the length. That requires a bit if extra math, which I'll skip here.


Possible to get transfer function coefficients from window?

What you are describing is an FIR filter, such that all the denominator coefficients are zero, save the basis, a[0]=1.  So you could do something like:
data_smoothed = scipy.signal.filtfilt(window, 1, data_noisy)

There is a notable point.  The DC gain of the filter is equal to the sum of the coefficients for FIR filters.  Your window is likely normalized to 1, so the sum is probably higher than 1, which means your filter will have gain at low frequencies.  You would want to divide all the coefficients by the sum of the window to keep the gain to unity.


My first suggestion is that you understand FFT in 1 dimension before trying to interpret results in 2D.

The Discrete Fourier Transform (FFT is an implementation of DFT) is a complex transform: it transforms between 2 vectors complex vectors of size N.

So in the 1D case, you will get not only negative values, but complex values in general.

The same applies in 2D. The imshow function is probably taking the real part of the complex matrix (it is not clear in the imshow documentation).


Your window is centered around the middle of the time window. This corresponds (roughly speaking) to a delay of half the window size and corresponding additional phase of $e^{-j \cdot \omega \cdot t}$

The "wiggles" that you see is a steep linear phase that wraps around 2*pi occasionally. There a few ways to fix this.

circularly shift the time domain signal so that the max of the window is at 0. This will be better but still isn't great since the windowing will significantly impact the phase.
Make sure that the fundamental fits into your time window and integer amount of times. Then do an FFT without windowing. This should give you an exact result.


Instead of calculating the STFT, you could just implement a 3 band filter bank as shown in the image below.

Your input signal x[n] could then be broken into 3 subbands: x0[n], x1[n], x2[n] which contain the energy associated with each frequency range. You could then feed each individual signal into an accumulator block in which the output would be the sum of the last x number of points. Or in you case, the x number of points that occurred in 50 ms. This would give you a single number in which you could map to some visualization output. You could implement normalization factors however you would like.

All you would have to do is define the filters H0, H1, and H2 as bandpass filters corresponding to your desired cutoff frequencies.


Determine the Signal Curve from Parameters of a Power Curve by Noisy Measurement

Since your domain for the parameters is limited (Only 9 options) the best way for White Noise would be going through them and pick then one with the least Mean Squared Error (MSE) which is the parameter to minimize for AWGN.
in MATLAB it will be something like:
function [ paramAlpha, paramBeta ] = EstimateModelParameters( vT, vY )

vParamAlpha = [1, 2, 3];
vParamBeta  = [1, 2, 3];

bestMse = 1e50;

vX = zeros(size(vY, 1), 1, class(vY));

for ii = 1:length(vParamAlpha)
    currParamAlpha = vParamAlpha(ii);
    for jj = 1:length(vParamBeta)
        currParamBeta = vParamBeta(jj);
        vX(:) = currParamAlpha * (vT .^ currParamBeta);
        currMse = mean((vX - vY) .^ 2);
        if(currMse < bestMse)
            bestMse     = currMse;
            paramAlpha  = currParamAlpha;
            paramBeta   = currParamBeta;
        end
    end
end

end

In a simple performance simulation I got the following:

Which means that for STD in reasonable value for this case (Since your maximum value is around ~3) your performance will be great.
The code is easily adaptable to any values of the parameters you'd like.
Though I think such brute force methods are reasonable for up to ~100 combinations.
Another thing you might try is use Affine Model on the log of the values.
Then apply some "Rounding" to the values. Though this will change the properties of the noise, in real world it seems to be effective way.
The full code is available on my StackExchange Signal Processing Q70753 GitHub Repository.


What is the phase and magnitude response of white noise?

Isn't white noise supposed to have a flat magnitude response?
      (equal amounts for all frequencies)

The expected magnitude response of white noise is flat (this is what JasonR calls the power spectral density).  Any particular instance of a white noise sequence will not have precisely flat response (this is what JasonR's comment refers to as the power spectrum).

In fact, the Fourier transform of white noise is... white noise!

  What is the relationship between the standard deviation
      (1 in my example) and the magnitude and phase?

There will be no relationship between the standard deviation and the phase.  As for the magnitude, suppose $n(t)$ is stationary white noise with zero mean and standard deviation $\sigma$. Then the autocorrelation (covariance) is:

$$
R_{nn}(\tau) = E[n(t)n(t+\tau)] = \sigma^2\delta(\tau)
$$

So the power spectral density is just $\sigma^2$ (though for discrete-time, there will be a scaling based on the duration of the signal).

Questions from comment:

  When you say that the Fourier transform is also white noise,
  how can I measure the std-dev when the transform is complex? Real,
  imaginary part or some combination?

Assume that our noise is discrete-time and is $n[m]$ (zero mean, Gaussian, white noise with variance $\sigma^2$).  Then the transform is:

$$
\begin{eqnarray}
N[k] &=& \sum_{m=0}^{M-1} n[m] e^{-j2\pi mk/M} \\
 &=& \sum_{m=0}^{M-1} n[m] \cos(2\pi mk/M) + j n[m] \sin(2\pi mk/M)
\end{eqnarray}
$$

and the expected value is:

$$
\begin{eqnarray}
E\left[N[k]\right] &=& E\left[\sum_{m=0}^{M-1} n[m] e^{-j2\pi mk/M} \right]\\
&=& \sum_{m=0}^{M-1} E\left[n[m]\right] e^{-j2\pi mk/M} \\
&=& 0
\end{eqnarray}
$$

The variance of the real part is given by:

$$
\begin{eqnarray}
E\left[ (\Re N[k])^2 \right] &=& E \left [ \sum_{m=0}^{M-1} n[m] \cos(2\pi mk/M)  \cdot \sum_{p=0}^{M-1} n[p] \cos(2\pi pk/M) \right ]\\
&=& E\left[ \sum_{m=0}^{M-1} \sum_{p=0}^{M-1} n[m]n[p] \delta[n-p] \cos(2\pi mk/M)\cos(2\pi pk/M) \right]\\
&=& \sum_{m=0}^{M-1} E\left[ n[m]^2 \right]\cos^2(2\pi mk/M)\\
&=& \sigma^2\sum_{m=0}^{M-1} \cos^2(2\pi mk/M)\\
&=& \sigma^2 \left( \frac{M}{2} +  \frac{\cos(M+1) 2\pi k/M \sin(2\pi Mk/M) }{ 2 \sin (2\pi k/M) } \ \ \ \right)\\
&=& \frac{\sigma^2 M}{2}
\end{eqnarray}
$$

I believe the imaginary part will behave the same way.

  Could you please enlighten me how the duration of the signal relates
  to the power spectral density (for discrete time situations)

I believe that (based on the above derivation), the power spectral density (the expected value of the square of the DFT) will scale linearly as the duration.

  If the phase is not affected by the std-dev, what determines the 3 degree amplitude, and the type of distribution (seems to be uniform rather than normal)

Check out the table on page 2 of this PDF file. it says that the argument (phase) of the coefficients will be uniformly distributed, as you state. Screenshot of the table included below.


Understanding the lowpass filtering of a digital signal

Any signal can be constructed from the sum of sinusoids of different frequencies (to within an arbitrarily small error).  This reconstruction is unique and can be calculated with the Fourier transform.

When we talk about the various frequencies that comprise a signal, we are talking about these sinusoids of various frequencies.  In fact, in signal processing, we are used to thinking about the same signal in the time domain -- as a function that gives a magnitude for every instant in time, or in the frequency domain -- as a function that gives a magnitude and phase for every frequency.

When you apply a low-pass filter, you are modifying the frequency domain view of the signal so that all frequencies higherthan the cutoff have their magnitudes set to pretty much zero.

Another way of thinking about it is, you take all the sinusoids that make up the original signal, and throw away the ones with frequency higher than your cutoff.

The secret to how this actually works is called the "convolution theorem": https://en.wikipedia.org/wiki/Convolution_theorem

The convolution theorem shows that convolving two signals together is the same as multiplying their frequency-domain representations together.  A "low pass" filter is just a signal that has a constant magnitude for frequencies below the cutoff, and near zero magnitude for frequencies above the cutoff.

When the filter is convolved with the input signal, the input signal's high-frequency components are multiplied by zero, effectively removing them.


Generating a noisy sine wave in Python, efficiently

The are some ways to improve your solution.

First, use numpy and vectorize your code:

def genSine(f0,fs,dur):
    sinusoid = []
    for t in range(dur):
        sinusoid.append(math.sin(2*math.pi*t(*freq/fs)))
    sinusoid = normalise(sinusoid)
    return sinusoid

becomes

import numpy as np
def genSine(f0, fs, dur):
    t = np.arange(dur)
    sinusoid = np.sin(2*np.pi*t*(f0/fs))
    sinusoid = normalise(sinusoid)
    return sinusoid

You will see the biggest improvment in efficiency when vectorizing all your code like this.

Second, concerning the normalization, you should consider the mathematics of your problem instead of doing it 3 times. You have 3 pieces of signal, the sinusoid, the noise, the mix. The first is already normalized (between -1 and 1), the second also (hence the name normal distribution for describing the noise). So you should only normalize the third.

Hope this helps!


Proper way to detect and decode transmission in gnuradio

It sounds like you are talking about a "bursty" transmission. This is common in real-world communications systems (for example, WiFi and Bluetooth operate this way).

What you say is exactly right - relying on power levels alone is not a reliable way of detecting a transmission. A real-world bursty communications system solves this problem by inserting some amount of completely "known" signal into every transmission. This carries no payload information whatsoever - it is specified by the communications standard (so both the transmitter and the receiver know it) purely to aid reception.

In many cases, this known redundant signal is inserted somewhere at the start of the transmission, in which case you might see it referred to as a "preamble". (However, in unusual cases, like TETRA which is what the police radios use in much of Europe, a "midamble" is used instead). In other contexts, this kind of known signal might be referred to as a "pilot" or "training sequence".

Although it might seem wasteful to transmit this redundant signal, it turns out to be very useful - not only for detecting a transmission, but also for performing other tasks that improve the reception (such as timing synchronization, frequency correction, gain control and channel estimation).

Therefore, if you have control over the data that is being modulated, then you should try putting some known data at the front of every transmission. Then, at the receiver, (in the simplest case) you can correlate using this known sequence in order to reliably detect the transmission. There is quite a lot of study behind how to choose a sequence that will produce really reliable detections in different scenarios. You could read about M-sequences, Gold sequences, Hadamard codes, Barker codes and many more.

If you don't have control over the transmissions, then you need to find some information about them. If you really know nothing about their structure, and can't estimate anything, then unfortunately power level detection may be the best you can do.


Does Function to Function Multiplication and Convolution, Follow Laws of Convolution

In a nutshell, you are asking whether the following statement is true:

$$[x(n)\cdot y(n)] *h(n) = [x(n)*h(n)] \cdot [y(n)*h(n)]$$

In words, if the convolution operator is distributive with multiplication.

The short answer is no.

An easy example would be to think of two sequences $h(n)$ and $x(n)$ such that their convolution is $0 \ \forall n$ (if you can't think of any, check this out). We can see that in this case the right side of the equality given above would be $0$ no matter what $y(n)$ is, but the left one would not (at least not for all $y(n)$).


How is noise added in an image during image acquisition or image transmission?

Imaging sensors are just electronic devices that receive photons (light) and turn them in to electrons (electric signal). As such, they are affected by ambient conditions such as temperature.
For example, you will have some noise coming from the electronic agitation (false currents due to freely moving electrons, without incoming light) that grows with the sensor temperature.

You do also have non-linear responses in the sensor that cause other noise patterns. A pixel needs to receive enough light before providing a reliable measurement. Before that threshold, it will be affected by Poisson noise (a noise whose distribution depends on the value of the pixel). On the other side of the per-pixel energy, it can get saturated (too much light), and its energy will "spread" onto neighboring pixels, thus modifying their value.

Finally, since you have a digital device, then you have an analog-to-digital conversion step that is likely to introduce further quantization noise, i.e., the recorded values will not be the "real" ones but something close up to the measurement precision.

Outside your imaging device, there is light transmission in the atmosphere (turbulence is very likely to be a disturbance cause here, but it's usually not considered as noise), and image transmission over communication channels where packet losses (and thus pixel loss!) can happen.


Why are there no poles for finite length signals

Hand-wavy? A rational fraction $$\frac{N(X)}{D(X)}$$
where $N(X)$ and $D(X)$ are polynomials can be written as some polynomial $P(X)$ plus a finite amount of terms of the shape:

$$\frac{a_i}{X-X_i}\,.$$

Let us say that a $0$ pole ($X_i=0$) is harmless, because it is
 just some $X^{-1}$. So  on the one hand, you have the monomials $X^d$, $d\ge0$, from $P$ and possibly a $X^{-1}$. If you replace the $X$ by a $z^{-n}$, you will still have a finite number of $z$ to the power of some integer, which are just delays.

On the other hand, if you have one non-zero pole:

$$\frac{1}{X-X_j} = \frac{1}{X}\left(1+ \left(\frac{X_j}{X}\right)^1+\left(\frac{X_j}{X}\right)^2+\cdots\right)$$
and replace the $X$ by some $z^{-n}$, you will get an infinite quantity of $z$ to the power of some integer, and thus an infinite impulse response (IIR). The complicated part is to show that several non-zero poles do not "cancel" each other. A lot of people confuse IIR and recursive.

Be cautious though that an apparently "recursive implementation" does not always imply non-zero poles: a FIR system with $z$-transform $1-z^{-2}+z^{-4}-z^{-6}$ can be written as:

$$ \frac{1-z^{-8}}{1+z^{-2}}$$

that is, you can rewrite: $$y ( k ) = x ( k ) - x ( k - 2 ) + x ( k - 4 ) - x ( k - 6 )$$ as $$y ( k ) = x ( k ) - x ( k - 8 ) - y ( k - 2 )\,,$$

which is a recursive difference equation, yet with a FIR.


Psudocode implementation of true envelope

As far as I can say, your interpretation is valid, but given the pseudo-code you provide, somehow, there seems to be some confusion with the definition of the quantities.

In my understanding, the $A_i(k)$ are the (real-valued) $\log$-amplitudes. For this reason, you do not need to compute the $\log$ of these values.

Then, the $V(k)$ are actually also $\log$-amplitudes. While I can understand to some extend that you use the same name both for $V$ and its inverse Fourier transform, this may be confusing when trying to understand the algorithm (and debugging it...). In addition, this could be an issue, depending whether the result of the IFT is real or not - this, in turn, depends on whether you "hermitianize" your vector $A$ before the IFT.

A tentative correction for your pseudo-code would be:

v = array of -9999999
a = log(mag(fft(x))
repeat a couple of times:
    for k from 0 to a.size():
        a[k] = max(a[k], v[k])

    vc = ifft(a) // let's assume vc is complex Hermitian, so we just keep the first half of the vector
    vc[from 50 to v.size()] = 0
    v = fft(vc) // here, be careful that the result is real

I did not try to implement it. Nevertheless, it seems to make sense that way. Let me know if that solves your problem!


Recommendation for book - Writing DSP code in C

I can recommend you two books about DSP for C language.

  Embree P. M. - C Language Algorithms for Digital Signal Processing

It is old and you can easily get it second-hand for a decent price. It covers pretty much all 4 topics that you described.

The other one I recommend is:

  Malepati H. - Digital Media Processing: DSP Algorithms Using C

It covers way more topics (including error correction algorithms and image/video processing algorithms). The nice thing about that one is that test signals are provided to test your implementation.

I just recalled that there is another book on my desk which has some C code for DSP in it.

  Press W. H. - Numerical Recipes: The Art of Scientific Computing

There are two chapters about the Fourier Transform and its applications.

One last suggestion - if you can, please use an existing C library instead of writing your own tools. There is no point in re-inventing the wheel.


Use of pulse shaping in digital modulation

QPSK in its raw form would appear just as the OP has shown, where the information is encoded into four symbols given by the four phase states alone. The problem is this has a very wide transmit spectrum given the pulse shapes are rectangular and the Fourier Transform of a rectangular pulse is a Sinc function (whose peaks roll-off in frequency only at the rate of 1/f). Wireless spectrum is an expensive commodity, so we go through great lengths in terrestrial wireless communications to limit spectral occupancy and in the process improve spectral efficiency. This is the only reason for pulse shaping, but it is significant. By more slowly transitioning between symbols rather than the abrupt transitions shown, we greatly limit the amount of spectrum needed (the spectrum will be given by the Fourier Transform of the pulse shape). Effective pulse shaping does this while controlling the inter-symbol interference that is added in the process.
The graphic below demonstrates this for a 16QAM waveform where the real portion of the time domain waveform is shown in the upper part of the graphic with and without pulse shaping. The lower part of the graphic shows comparative spectrum where we see the significant reduction in spectral occupancy due to the pulse shaping provided.

I see this graphic is the same one as what the OP linked, so to further explain: This is how the envelope of the signal would look like in the air (we don't see the actual carrier frequency in this graphic). This graphic was for QAM so does have information in the amplitude as well as phase, but for QPSK case the amplitude would carry no information. But even with the QAM case, the pulse shaping modifies the amplitude ONLY to reduce the spectrum, and modifies it in such a way to control the trajectory between our symbol samples of interest, while still passing through those samples defined by the modulation exactly at the proper symbol timing locations.
Below is an eye diagram showing the real and complex waveform for the QPSK case (The eye-diagram similar to the time domain waveform shown above, but repeating synchronous to symbol boundaries)

Show this shows how the over the air waveform would look (in actuality we typically do half the pulse shaping in the transmitter with a root-raised cosine filter while this plot shows the result after the second root-raised cosine filter in the receiver), but in general we are seeing the amplitude and phase of the time domain waveform for all samples, and the job of the receiver is to determine the sampling locations where the red dots are shown (the center of each symbol), and we then determine the magnitude and phase of those to demodulate (or demap) the QPSK data.
To see what is transmitted "over the air", the sketch below begins to demonstrate how the envelope of the carrier would be modulated by the pulse shaping, but doesn't provide the whole picture as the phase would also be transitioning from symbol to symbol so would not align with what would occur if we continued to transmit one symbol repeatedly. I added more details further below showing the exact waveforms expected. These are the real waveforms that are transmitted at the antenna rather than the baseband equivalent waveforms showing the complex envelope typically with two plots as I (real) and Q (imaginary). The sinusoidal carrier could be any frequency and the envelope would not change, including frequency = 0 which would then represent the eye diagram and I/Q plots that I have shown above (and must then be a complex waveform to represent).

A further note is the OP used symbol mapping that is not optimally gray coded (the 180 degree transition depicted in the diagram going from symbol 00 to 10 should be reserved for transitions where both bits change such that only one bit changes for all symbols that are close to each other).
Also note how this is similar to the operation of windowing prior to computing the FFT to reduce spectral leakage.

To avoid any confusion based on subsequent questions in the comments, here is the actual results on what the OP's QPSK waveform would look like at "RF" using properly gray-coded symbols using the following mapping of symbol to constellation:
0: -1-1j
1: 1-1j
2: -1+1j
3: 1+1j
We will start with the final waveform in the receiver as the desired QPSK constellation and work back to what is transmitted over the air by adding in the carrier offset and seeing the result:
First the I and Q waveform at baseband for a QPSK waveform pulse-shaped with a raised cosine filter with alpha = 0.3 would appear as the following in the receiver once all the carrier offsets are removed:

The constellation of the above is formed by plotting I vs Q on the complex plane. These trajectories represent the phase and amplitude of our RF carrier at every moment in time, and here I show the one instance that occurs at every integer symbol number in red, which is the only point over the duration of each symbol where the phase and magnitudes of the carrier will actually represent what we expect for QPSK (every other sample in between will be in transition at other phases and amplitudes!).

With an actual carrier frequency added, this same sequence would become the waveform shown in the plot below, where the amplitude of the envelope will match the amplitude of the constellation in the plot above:

Zooming in on the first 10 samples, this is what the OP would see for the sequence given (this is for the sequence 0,2,1,3,0,2,1,3,1,1,2,2.... specifically). The symbol sample points are aligned in this plot with the integer symbol locations on the horizontal axis.  This was done with a carrier frequency that is 10 times the symbol rate, but if we increased the carrier frequency further this plot showing the envelope would just get solid but appear identically the same otherwise.

And the same thing as we reduce the carrier to 4 times the symbol rate, showing that we will see the same basic envelope shape regardless of which carrier frequency is used (and as noted if the carrier gets too small compared to the symbol rate we need to use a complex signal to properly represent the waveform, exactly as we do in transceiver implementations):

Note that these are all done with a raised cosine filtered pulse-shape, but what is typically transmitted over the air is root-raised cosine filtered.  To be complete, below is the plot showing the same waveform with a root-raised cosine pulse shape:

With close observation of the above plot against the unmodulated carrier aligned with the first symbol we see how the modulated carrier is continuously changing in phase during the course of each symbol and, as also made clear by the constellation diagram of the modulated waveform, is only at the correct amplitude and phase according to the QPSK constellation only over very short durations, and in this specific case after being passed through a second RRC filter. (If the pulse-shape is altogether eliminated then the waveform will stay at the correct amplitude and phase for the entire symbol, and as the pulse shape roll-off factor is increased, this will approach that).

The above plots and descriptions used Raised Cosine (and Root-Raised Cosine) pulse shaping but an alternate harris-Moerder filter is described in this paper http://eon.sdsu.edu/~seshagir/SDR05.pdf. Using this as the pulse shaping filter results in an order of magnitude lower EVM for the same number of taps. This isn't commonly used even today (15 years later) to my knowledge since the root-raised cosine filter is baked into our specifications but certainly is very compelling to be aware of for consideration in future radios (I added a question to see if there is further experience with its use: Use of the harris-Moerder Nyquist Pulse Shaping Filter which may have additional information on this.)


Suppose in time domain a Gaussian kernel is expressed as:
f(x,y) = exp (-r.^2/r0.^2);

where r.^2 = x.^2 + y.^2;, r0 is the radius of the exp(-1) point.
Then the 2D FFT should be:
F(u,v) = pi/r0.^2 * exp(-w.^2/w0.^2);

where w.^2 = u.^2 + v.^2; , w0 = 1/ (pi * r0);
As a result, I would change your function to:
function H = gaussian_kernel(size, sigma)

    half_size = floor(size/2);
    [u,v]=meshgrid(-half_size:half_size,-half_size:half_size);
    H = exp(-(u.^2+v.^2) ./ (2 * sigma^2));
    r0 = 2 * sigma^2 * pi^2;
    H = H * r0 * pi;
end

And,
H = gaussian_kernel(16, 2);
subplot(2,1,1),imagesc(H)  % frequency domain
subplot(2,1,2),imagesc(real(fftshift((ifft2(fftshift(H))))))  % time domain

result:

Suppose it is (-N/2+1 : N/2) /N * fs in the frequency axis (N is the sampling point number, and fs is the sampling rate), then it is supposed to be (0:N-1)/(N * fs) in spatial axis.


Audio distortion while implementing a low-pass filter

I suspect that most likely the audio coming in is integer, and not floating point. Depending on your device, operating system and type of file input, you can get one or the other. Read documentation for audio input API you're using to see which format the audio is in. Also, if it's integer, it's likely to be 16 bit rather than 32 bit. In any case, again, read the documentation.

EDIT:
BASS library has a flag you can set to get floating point values everywhere:

BASS_SetConfig(
    BASS_CONFIG_FLOATDSP,
    BOOL floatdsp
);


How to apply baseline methods to MNE preprocessed raw EEG signals?

MNE has a fantastic website with a lot of documentation and many tutorial examples. Those examples include how to access machine learning classifiers in scikit-learn (to play with SVMs) as well as the methods for getting / plotting the power, extracting frequency components with FFT. You should definitely check out the website.

Alex


Why do I see ringing in the output of a digital filter with a narrow transition band?

Your observations are an example of the Gibbs phenomenon. When you apply a filter with a very sharp transition band, you will observe oscillations in the filter output (or "ringing") near any sharp transitions in the input signal (e.g. boundaries of pulsed waveforms). The apparent "frequency" of the oscillations is dependent upon the bandwidth of the filter; as you increase the filter's cutoff frequency, the oscillations will become more localized in time (i.e. "higher in frequency"), but the peak overshoot does not change. The Wikipedia article linked above has a good explanation halfway through or so.


What is the name for a constant-heading Kalman filter model for vehicle tracking?

I think the magic acronym is CHCV, "constant heading constant velocity". This returns at least a few results on Google.


Count the vehicles flow that pass over one bridge

It is relatively easy to create the "sensing element" for the cars but difficult to get a short term valid count for the number of vehicles.

As far as the sensing element is concerned, this is a simple pulse counter. The signal that creates the "pulse", comes from the optical flow of the vehicles as they pass by a region of interest in front of the camera.

You basically need to differentiate between successive frames, rectify the differentiated signal and integrate it with respect to time with the possible addition of a threshold.

When you differentiate between two successive frames (for instance $I_n, I_{n-1}$, where $I$ is a two-dimensional array of pixel brightness values at some color depth and $n$ is the frame number (integer))  you are basically removing everything that is not moving and highlighting parts of the image that appear to have changed between time instances $n$ and $n-1$.

At this point, you don't really care about the direction of the change and this is why you need to take the absolute difference: $Q_n = |I_n-I_{n-1}|$

If you now focus on a small region of interest of the $Q$ series of images, you will notice that when nothing happens, that region is pitch black and when an object MOVES in front of it the region "lights up" because of all the differences. Therefore, what you can do now is obtain the mean value of all the pixels of this small region of interest (over the $Q$ sequence) and observe it in time.

This will be producing a "pulse" every time that something disturbs the region of interest. The duration of the pulse (in time) will be equal to the duration of the disturbance.

Now, if we were to think in terms of "counting cars over a bridge", you can position your region of interest around a light pole. Past the first frame, the actual scene will not matter anymore because the light frame doesn't move (actually it does move because the whole bridge vibrates slightly and the light poles are tuning in but let's assume that this is negligible in your case). If you do this then you can count disturbances under the light pole during the day and also during the night in a similar way (rather than having to rely on the irregular shapes of the car lights).

Fine, so now we have a way to transform passing vehicles into "pulses". What is the problem? Count the pulses and you get the number of vehicles.

Not exactly, so here is where things become difficult:

There are more than one lanes per direction. Therefore, it is quite possible that two cars fit through the region of interest at the same time thus creating just one pulse. Even worse, a car can be masked by a bus (be behind it) and never be counted. The only way to counteract this is to repeat the same ultra simple "counting" process at multiple points on the bridge (multiple light poles) and obtain an average count. (Of course, a car can travel at the same speed as a bus and remain hidden from the perspective of the camera for as long as it remains on the bridge...but there is no remedy for this....from your point of view).
As you can see, this pulse counting is actually counting disturbances, not necessarily vehicles. Here are a few disturbances that have nothing to do with the vehicles: Atmospheric disturbance (smoke, sea evaporation, anything that causes a slight flickering), movement of the camera, birds and light signals. Imagine for example what happens when the bridge lights will come on. Your pulse counter will register a pulse (because before, the area under the light pole was darker and now, it is lighter) but this pulse is going to be extremely brief. The only way to counteract spurious movement somewhere across the image is by the solution to #1 (multiple counting points + average) and the only way to counteract invalid pulses is by setting thresholds. So, very short pulses (in time) are not counted and very weak pulses (in terms of magnitude) are not counted too.

This is the principle. You can do this with openCV, python, (openCV through Python), MATLAB, even off-line video processing.

Hope this helps.


Compressed / Compressive Sensing - Sensor Placement

Suppose your beam has ten vertical degree of freedom. You can only apply CS to compress measurement vector of each DOF. You can not reconstruct all of the responses using limited number of them.
The ten-dof structure would have ten pick- related to ten natural frequency- in frequency domain, so its sparsity has been already guaranteed. Hope this help.


NDA timing estimation for mutli-level QAM(RC pulse)

Consider using a Gardner Timing Error Detector which in the following form is usable for higher order QAM:

$$TED = I_n(I_{n+1}-I_{n-1}) + Q_n(Q_{n+1}-Q_{n-1})$$

Where $I_{n-1}$,$Q_{n-1}$,  $I_{n}$, $Q_{n}$ and $I_{n+1}$, $Q_{n+1}$ are the early, prompt and late QAM samples at 2 samples per symbol. The TED will drive the prompt sample to zero error at the zero crossings of the waveform.

This is from deriving a timing error by taking the real of the complex conjugate multiplication of the signal y[n] with its derivative, a form which comes from simplifying the maximum likelihood timing estimator. The difference $y_{n+1}-y_{n-1}$ is an approximation of the derivative of y[n]:

$$real\bigl(y^*_n (y_{n+1}-y_{n-1})\bigr)$$

This error term would typically be the input into a loop filter in a timing recovery loop, and the output of the loop filter could select the appropriately timed output of a polyphase resampling structure such that the sampling rate is never actually changed (2 samples per symbol throughout).

Assuming the RC pulse is root-raised cosine (RRC) in the transmitter, and then again in the receiver, the Gardner TED has higher performance using the samples prior to the RRC filter in the receiver: this makes sense as the TED derives the timing error or S-curve based on the symbol transitions- if you observe the zero-crossing jitter in an eye diagram of RRC pulses compared to RC pulses with the same pulse-shaping parameter, you will see that the RRC has less jitter. The second RRC filter adds more jitter at the benefit of eliminating the inter-symbol interference at the symbol sampling locations.

That said, the waveform can be prefiltered prior to the timing error detection to minimize zero crossing jitter (this filtered waveform is used only for timing and is not the received waveform for data recovery due to ehanced ISI at the correct symbol decision locations) Here is an interesting paper on the use of such noise minimization pre-filters for QAM with the Gardner Ted:

Aldo Nunzio D’Andrea, Marco Luise, Optimization of Symbol Timing Recovery for QAM Data Demodulators, IEEE Trans on Communications, Vol 44 No 3 March 1996.

https://ieeexplore.ieee.org/document/486334

Further, here is an interesting graph I created showing how robust the TED is without carrier/phase synchronization. This shows how the TED performs as an SNR metric with carrier offsets normalized to the symbol rate, showing only 3 dB degradation with carrier offsets as much as a quarter of the symbol rate. (The "Gain" is from the slope of the discriminator, while the noise is the pattern noise of the discriminator. For the following two graphs these were just consistent metrics to compare various timing discriminators such as here the Gardner TED versus the M&M synchronizer.)

For comparison here is the same plot with the Mueller and Mueller (M&M) synchronizer. Notice the $10^{-4}$ scaling on the horizontal axis! This is overlaid on the plot above with the red line showing the same TED curve, just zoomed in on the horizontal axis; What this shows is that the M&M is much more sensitive to carrier offset (as expected since it operates on decisions) but has lower noise once the carrier is synchronized.

Here is a simplified plot showing how the Gardner TED estalishes a timing discriminator "S-curve" from the product of the signal with its derivative. This shows how the Gardner uses the zero crossing transitions to establish the discriminator curve (the actual curve is sinusoidal for raised cosine waveforms) and importantly how the lock point is at the zero crossings (blue samples) while the decision samples would be the other samples shown in red (with 2 samples per symbol). For higher order QAM the functional operation is the same with substantially higher pattern noise. (The timing loop BW would need to be reviewed to confirm the resulting discriminator noise is sufficiently averaged to not be an issue).


Reconstructing a signal from FFT by adding individual signal components

There are a couple of mistakes in your calculations:

The t array is incorrectly specified. It should be 0:1/8:7/8 as the time positions at which the signal is assumed to be sampled for the default FFT calculation performed by fft() are 0, 1, ..., (N-1) where N = 8 in your case, the signal length. Since in your formula you are using t as the ratio n/N (where n is the array of sample times) you get the aforementioned expression to use for t.
You wrote a(1) instead of a(5) at the last term of your recs expression.
You are missing 2* multiplying each term except the first (k=0) and last (k=N/2) ones, where k is the index of the frequency axis of the Fourier transform F (recall that this 2* comes from leveraging the symmetry of the Fourier transform associated to a real signal).

Therefore, your reconstruction expression should read:

recs =  (1*a(1)/8)*cos(0*pi*t + an(1)) + ...
        (2*a(2)/8)*cos(2*pi*t + an(2)) + ...
        (2*a(3)/8)*cos(4*pi*t + an(3)) + ...
        (2*a(4)/8)*cos(6*pi*t + an(4)) + ...
        (1*a(5)/8)*cos(8*pi*t + an(5));

which gives the original signal s as verified by computing:

recs - s

ans =

   1.0e-14 *

         0   -0.1332         0   -0.0444   -0.2665    0.3553    0.1332         0


Jamming resistance in spread spectrum systems

In order to work well, a jamming signal must be well-correlated with the signal it is trying to jam.  In the case of sine-wave modulated signals, it's easy to detect the signal one wishes to jam, to infer that it's narrowband (and, hence, riding on a real or suppressed sine-wave carrier), to infer its bandwidth, and from there to choose a suitable jamming signal.

To jam a CDMA spread spectrum signal you'd either need to overwhelm it with an uncorrelated signal (random, single tone, pseudo-random noise) or you'd have to find it, find its PN sequence, synchronize to its PN sequence, and then in most cases adjust your synchronization to match the time offset between your reception of the signal and the receiver's, then you'd have something to send -- but if I were trying to receive the signal, I could defeat you by having multiple receivers, located more than a wavelength of the chipping rate away from each other, each trying to receive the signal.  One receiver may get the jamming signal and the intended signal on top of each other, but unless your jamming signal is co-located with the intended transmission, other receivers would be able to pick out the good from the bad.

But keep in mind that civilian CDMA is intended to minimize unintentional jamming, and does pretty well.  It's not intended to defeat malicious jamming.

I suspect that there are smart military types out there that know a lot more than I do about jamming CDMA -- I'm just making this up as I go from first principles.  I also suspect that military systems will use something more complicated than civilian CDMA, for just that reason.


The notation is a little obscure but the answer is simple.

Since $\delta[n-8]^2$ can be represented as : $$\delta[n-8]^2 = \delta[n-8] \delta[n-8]$$ and since the multiplication of the impulses can be simplified as :
$$\delta[n-8] \delta[n-8] = \delta[n-8]$$ , then you can conlcude that the $\mathcal{z}$ transform of the signal $4 \delta[n-8]^2 $ will be:

$$ \mathcal{Z}\{ 4 \delta[n-8]^2 \} = 4 z^{-8} $$ whose region of convergence is all $z$ except $z = 0$.

There is a pole at $z = 0$ and a zero at $ z \to \infty $


Do we really need 10 samples of each word? or can we "combine" them into a single template? and if so, how?

You can use single template, 10 templates will make result more accurate.

  Should we run the SVM directly on the MFCCs vectors? On the MFCC vector of vectors (for each word)? On the DTW values? or should we combine all the MFCC vectors into a single one?

SVM has no meaning for your application because it is not invariant to time stretch so it's not suitable for classification of time series. There are SVM modifications to deal with time series, for example

Hidden Markov Support Vector Machines

but it's completely different story

  If we should look for a minimal match for the DTW between the input and the templates, should we give different weights to the templates depending on length, or even account for the differences between the templates (the templates of the word "1" are more distant from one another than the templates of the word "2").

No, it's not necessary. A simple comparison of DTW values should work.

  All three of these methods didn't work well (and gave almost everything wrong results).

Probably you made a mistake in DTW implementation first of all. You need to verify if you implemented everything correctly. You might make a mistake in DTW and go to incorrect conclusion

  or references to good sources.

Read a good text book. At least Rabiner "Fundamentals of Speech Recognition"