Is there a difference between filtering a signal before or after differentiating it?

So, first of all:

A differentiator is really just a high-pass filter. Digitally (and this is dsp.SE, so I presume this is the case), that typically means a differentiator is really just FIR with taps $[1;-1]$. A butterworth filter being a low-pass IIR, the combination of both does sound like you might want to build a bandpass, or are doing something that doesn't make very much sense; see the comparison of a very relaxed second order butterworth and a differentiator magnitude responses below:

now, mathematically, that's identical:

In time domain, you convolve ($*$) the signal $x[n]$ with the filter impulse response:

$$\begin{array}\\
y[n] &= x[n] * h_{diff} * h_{butter}\\
&\text{or, in freq. domain, thanks to the convolution theorem}\\
Y[z] &= X[z] H_{diff} H_{butter}\\
&\text{multiplication is commutative}\\
&= X[z] H_{butter} H_{diff} \\
&\text{but it's also associative}\\
&= X[z](H_{butter}H_{diff})\\
&= X[z] \mathcal Z\{h_{butter}*h_{diff}\}
\end{array}$$
which implies you can even without problem combine the differentiator and the low pass filter. To little surprise, you get the discrete derivative of the butterworth filter.


What is impulse response (or equivalently the frequency response) of the channel ?

This is a fundamental topic and it needs a lot of explaination.
I suggest you use the fading channel object of MATLAB, especially the Channel Visualization Tool.

There are lots of examples with visualization in this link. MATLAB channel OBJECT


Why does python's scipy.signal.dimpulse introduce delay in impulse response?

In lfilter the transfer function is described in decreasing powers of z, as shown below copied from the python doc for signal.lfilter:

                    -1              -M
         b[0] + b[1]z  + ... + b[M] z
Y(z) = -------------------------------- X(z)
                    -1              -N
         a[0] + a[1]z  + ... + a[N] z

While the parameters B and A in dimpulse are represented in increasing powers of z as described in the doc for signal.TransferFunction:

  If (numerator, denominator) is passed in for *system, coefficients for
  both the numerator and denominator should be specified in descending
  exponent order (e.g. s^2 + 3s + 5 or z^2 + 3z + 5 would be represented
  as [1, 3, 5])

Along with the following example given for:

$$\frac{z^2+3z+3}{z^2+2z+1}$$

Contruct the transfer function with a sampling time of 0.1 seconds:

H(z)=z2+3z+3z2+2z+1
H(z)=z2+3z+3z2+2z+1

>>>>signal.TransferFunction(num, den, dt=0.1)
TransferFunctionDiscrete(
array([ 1.,  3.,  3.]),
array([ 1.,  2.,  1.]),
dt: 0.1
)

So in your case, for dimpulse:

$$H_1(z)=\frac{1-0.9z^{-1}}{0.1}$$

And for lfilter:

$$H_2(z)=\frac{z-0.9}{0.1}$$

So

$$H_1(z)= H_2(z)z^{-1}$$

And mystery solved!


Method for determining probe angle by analyzing skewed sine wave

I will assume:

The probe is like a line, infinitely narrow.
We move the probe along the $x$ dimension, parallel to the large-scale surface, and from each position, probe the surface:
The probe is extended to the direction it is pointing to until it meets the surface. (An alternative assumption would be that probing is done to the direction $y$, that of the normal of the large-scale surface, in which case the waveform measured would be sinusoidal unless the probe stem meets the surface before the tip. I wonder if your edit describes this scenario.)

Equation of line starting at coordinates $(x_0, 0)$ and with a positive slope $k = \frac{\Delta y}{\Delta x} \ge 1$ (where $1$ is the maximum slope of the sinusoid of Eq. 2, to ensure continuity):

$$y = k (x - x_0)\tag1\\
\Rightarrow x = x_0 + \frac yk.$$

Equation of the first period of the constant-biased sinusoid, spanning $0 \le x \le 2\pi$:

$$y = 1 - \cos(x)\tag2\\
\Rightarrow x = \begin{cases}
\arccos(1 - y)&\text{if }0 \le x \le \pi\\
2\pi - \arccos(1 - y)&\text{if }\pi \le x \le 2\pi.
\end{cases}$$

By Eq. 2 the period starts at $(0, 0)$ and ends at $(2\pi, 0)$, which makes the corresponding period of the skewed sinusoid also start and end at these points. We obtain the waveform by solving $y$ as function of $x_0$ from the system of Eqs. 1 and 2. Unfortunately there is no closed form solution, but we can still make some plots:

Figure 1. The skewed waveform for different values of the constant $k$. At $k=\infty$, the waveform is equal to the constant-biased sinusoid. The waveform shown represents the probe tip $y$ position as function of the coordinate $x_0$ at which the probe shaft crosses the $y$ axis. If you desire the "probe length" waveform, you need to multiply this waveform by a factor $\tfrac{\sqrt{k^2 + 1}}{k}$.

I doubt that you will be able to find an analytical solution in frequency domain. What we can do is determine $k$ based on the waveform peak location $x_0 = x_p$ by setting $y = 2$, for which the system of Eqs. 1 and 2 gives:

$$k = \frac{2}{\pi - x_p}.\tag{3}$$

To get the angle $\phi$ by which the probe angle deviates (moderately) from the normal of the large-scale surface:

$$\phi = \operatorname{atan2}\left(1, k\right) = \operatorname{atan2}\left(\pi - x_p, 2\right).\tag{4}$$


An set of digitized intensity samples that varies in time produces a waveform (just graph the data).  If the waveform shows any periodicity, that will be an obvious frequency.  Even if it doesn't, Fourier theory says any (non-pathological) waveform can be decomposed into a set (perhaps infinite) of sinewaves.  And sinewaves are describe as having a frequency.  Thus frequency information is implicit.

And this frequency information can be changed by certain filtering processes or algorithms many different ways.

Even if the digitized samples don't vary at all (a set of constants), that is the same as DC, or a frequency of zero (0 Hz).


Can Convolution neutral network train(learn) separately (train different times)?

This is called curriculum learning , and not only this approach is possible but also it is shown to be quite effective.
A few years ago I also came up with similar idea. Basically the idea is, instead of flooding the network with different examples, we can categorize the training data from simple to hard and gradually feed data ordered based on difficulty.
After searching for a while, I encountered this paper by Yoshua Bengio, which i think is among the important ones on curriculum learning. Also check out
FlowNet 2.0 which is a network trained with this approach.


Guitar pitch detection with autocorrelation

Cons: Not as accurate

This is just compared to the other methods.  I was measuring frequency very accurately to look for clock drift, etc: 1000.000004 Hz for 1000 Hz, for instance.  For guitar pitch detection it will be fine.

  doesn't work for inharmonic things like musical instruments

I should have said "it can't find an accurate fundamental if there is inharmonicity".  It will be shifted up because the plucked string partials are slightly sharp.

In your case, you don't want to find the fundamental, though, you want to find the perceived pitch, which is what autocorrelation is good at.  (though I don't know how accurate, and if it needs weighting, etc.)

http://www.utdallas.edu/~assmann/hcs6367/yost09.pdf

https://ccrma.stanford.edu/~dattorro/Humans.pdf

  this implementation has trouble with finding the true peak

because autocorrelation will produce peaks at the correct frequency but also at all its sub-harmonics, and this function just naively chooses the highest peak.  You need to plot the correlation curve and find a way to choose the correct results for your waveforms.


Best beginner DSP book for developing plugins

There are a number of books in addition to Zölzer's DAFX, the one I would recommend is:

The Audio Programming Book by Boulanger...

Also have a look at these:

Designing Sound by Farnell

Designing Audio Effect Plug-Ins in C++ by Pirkle

Designing Software Synthesizer Plug-Ins in C++ by Pirkle


From a simple Google search...

http://sourceforge.net/adobe/adobedatasets/panoramas/home/Home/
http://www.ti.uni-bielefeld.de/html/research/avardy/index.html
http://sun360.mit.edu/
http://www.iiia.csic.es/~aramisa/datasets/iiiapanos.html
http://www.visualsize.com/mosaic3d/index.php


Magnification of a colored image through two color space: RGB & YCbCr

I have checked your code and it worked fine in my matab 6.0 implementation with the following lines of of modification:

A=double(imread('Sun.jpg'));  % convert uint8 data to double
...                           % further operations here
A1(A1(:)<0.0) = 0.0;          % check and correct out of range values
A1(A1(:)>255.0) = 255.0;      % check and correct out of range values
imwrite(A1/255.0,'Sun1.jpg'); % while writing the data of class double use [0,1] data range

According to this, as I said you must deal with numerical issues such as an out of range value.


How can I trace a path around an image consisting of only black and white lines?

Depending on the accuracy of the trace you wish to achieve, there are a couple of ways you can do this:

The most straightforward way is to:

Convert the image to grayscale
Threshold, to achieve a completely binary image
Apply a high pass mask to get rid of the constant areas

For example $\begin{matrix} 1 & 1 & 1 \\ 1 & -9 & 1 \\ 1 & 1 & 1 \end{matrix}$

The result looks like this:

This technique is not without limitations, because as you can see, Pink Panther's whiskers have been drawn boldly and therefore have been "traced" and the highlight in his left eye has been interpreted as "background", producing an inner trace. To reduce these artifacts you can use a "wider" filter (for example, a $5 \times 5$ square of $1$s with a $-25$ in the middle. This produces a stronger wider trace line), morphological operations or adaptive application of the mask depending on features of the patch you are looking at.

The "path based" tracer you seem to suggest is how active contours work. In this case, you would initialise the active contour on the image rectangle and then iterate it towards the outer edges of the cartoon. Scikit-image has an implementation of active contours.

But based on the work you have done already, you could switch from a "Traveling salesman" route discovery to a convex hull of the points you collect.

The convex hull will trace an outline that encloses the set of points and give it a "cartoonistic" sketch view as well. Again, scipy has a convex hull implementation and so does Octave, from which I am posting this indicative output.

You will need to have the image package loaded for the following code to work and I am re-using the last image of the traced Pink Panther from above to extract the points, but that is just because I had it already, you can simply feed the convex hull the points you have extracted right after thresholding. It might be a lot of data indeed but convex hull is not solving a TSP and returns much quicker.

I = imread("thetracedimage.png");  % Load the image
Z = find(I);  % Find the outline (anything that is nonzero)
[Qy, Qx] = ind2sub(size(I),Z);  % Convert the index values returned by find, to image coordinates (in Qx, Qy) NOTICE ORDER OF X,Y!
H = convhull(Qx, -Qy);  % Obtain the convex hull
% H is the index of the existing points that make up the convex hull and we are assuming that these are connected by lines. To retrieve the actual points:
Fx=Qx(H);
Fy=Qy(H);
% Now plot everything. Notice the reversal of y purely because of the plotting coordinate system that puts 0,0 at the bottom left corner.
plot(Qx, -Qy,'.',Fx,-Fy);axis image;

And this results to:

Which of course is less accurate, but it depends on what you are trying to do and you seem to indicate that the approximate output of the TSP is sort of satisfactory (?).

Hope this helps.


Since it is very hard to guess what result you need, proposing any automatic segmentation method is very hard.

Let me instead propose two different semi-automatic segmentation with markers.

Segmentation using Minimum Spanning Forests with markers

This method is based on a version of watershed [Cousty et. al.: Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle, 2009].

Basically, you calculate a version of watershed (there are many versions), called watershed cuts, which gives you a (pretty severe) over-segmentation.

After that, you mark your one object (or all of your objects, 1 marker per object + 1 for the background if there is any) you want to find, and it calculates the most likely segmentation based on pre-calculated watershed.

I don't want to go into details, since there is pseudo-code and example images provided (one is even a cell image) in the article, but if you want me please leave a comment and I can add some more information.
Segmentation using $\alpha$-trees (or $\alpha$-connected components) with markers

Here, the "starting over-segmentation" is actually simply the segmentation to flat-zones (equal-intensity areas in the image). I am going to explain a little bit of theory here.

An $\alpha$-connected component of a pixel $p$, or $\alpha$-$CC(p)$ is obtained as follows:

first you choose a connectivity (4- or 8-conn.; 4-conn simply means the neighbors of a pixel are up, down, left, and right)
if the intensity difference between two neighboring pixels is $\leq \alpha$, those pixels are connected
$\alpha$-$CC(p)$ contains all the pixels connected to pixel $p$ as well as all the pixels that can be reached from pixel $p$ through a connected path

Now, for various reasons (I can expand the answer again if somebody asks), if you follow the $\alpha$-$CC$ around one pixel, its area usually grows slowly as you increase the $\alpha$, then at some point you have a sudden and rapid increase in the area, and then it continues to grow slowly. Usually, if you look at the $\alpha$-$CC$ just in the moment before the rapid size increase, it will give you your object.

Again, you can use markers to mark which $\alpha$-$CC$'s you want to "follow". The advantage here is that you can calculate the $\alpha$-hierarchy once and then do all the calculations on the same hierarchy.

I know it sounds a bit confusing with all the new terms unless you know something about $\alpha$-trees, but here's some references to help:

Najman, Soille: On Morphological Hierarchical Representations for Image Processing and Spatial Data Clustering, 2012 - a good reference for $\alpha$-hierarchies in general
Serna, Marcotegui: Attribute Controlled Reconstruction and Adaptive Mathematical Morphology, 2013 - these guys have done something very similar to what I am proposing. One of the proposed and explored usages in their article is segmentation (they're just using more attributes, not just area).
Merciol, Lefevre: Fast Image and Video Segmentation Based on Alpha-tree Multiscale Representation, 2012 - this is more or less exactly the usage I described.


Let's do the following, in this order:

figure out how much power your signal has at the beginning of the stop band (i.e. 5 MHz away from the lowest or highest of the 1200 used carriers)
compare that to how much attenuation you need relative to the peak power region
convert the result to a filter specification

How much power does this OFDM system have 5 MHz from the used carriers?

I hope that number's right, but as far as I could find, the LTE subcarrier spacing is 15 kHz.

That means that within your 5 MHz transition width, 333 subcarriers fit.

Let's look at each sinc (the spectral shape of an OFDM subcarrier!) individually first:

The energy in its sidelobe monotonously falls with the "number" of the side lobe.
The exact value is rather hard to come up with, so I just wrote a python script that simply calculates the energy for the sidelobe in the middle between the sinc zeros; it's rather compact, so here it goes:

def f(k, N=1200):
    positions = numpy.arange(k, k+N, 1)+0.5
    return numpy.sum(numpy.sinc(positions)**2)

So, calling f(13) will give you the sum energy of all the 1200 sidelobes that fall into the 13. sidelobe centre from the edge of the band of used subcarriers; that's a slight miscalculation for the maximum energy in the sidelobe (as the maxima don't exactly hit the center between zeros), but the error should be negligible. (If you think it's not: calculate the probability of the event "one of the two symbols that lead to the actual maximum energy is transmitted" and see that we really don't have to be afraid of regulators knocking on our doors.)

So, plotting that in a dB scale yielded this fine figure:

The orange vertical marks the 333. sidelobe – which, if the 15 kHz subcarrier spacing is right, is the sidelobe at 5 MHz from the highest used carrier.

Compare to how much attenuation you need

From the simulation, the 333. sidelobe has -36 dB power relative to main lobe.

You need -34 dB attenuation.

Hence, to little surprise, the number of occupied carriers (1200) relative to the full bandwidth of the system was chosen exactly to guarantee sufficient out-of-band suppression.

Convert the result to a filter specification

As predicted, an LTE signal fulfills the LTE spectral mask by itself.

You hence don't need any filters to fulfill the spectral mask.

Hence, you also shouldn't use any filters, to allow the flexibility of the OFDM system to correct the channel, not your superfluous filtering.

Hence, the filter specification is: do a FIR filter with taps [1.0]. Or don't do one at all ;)


Matched Filter BER vs SNR curve does not match the theoretical curve

I dont know, if this is exactly what you wanted to do, but here's a script that creates your theoretic curve. I used your mirrored coefficients as the transmit pulse, and the coefficients as the filter coefficients. This way, the pair becomes a matched filter pair.

N=10^4;        %No: of bits
 A=10;          %amplitude of S(t)
 T=length(v);          %duration of S(t)
 S=fliplr(v)*A;       %the transmit pulse
 E=norm(S)^2;         %S(t) Energy
 input = rand(1,N)>0.5;   % generating 0,1 with equal probability
 bits=(2*input-1);      %mapping the bits with S(t)
 x=kron(bits,S);

 SNR_db=0:0.25:10;
 n_var=(E/2)*10.^(-SNR_db/10);

 %Noise that gives SNR=0dB to SNR=7dB

 for i = 1:length(n_var)

  n=randn(1,length(x))*sqrt(n_var(i));   %noise
   r=x+n; %received signal

  %v=ternary coefficients, filtering the received signal and coefficients
   filtered=filter(v,1,r);

   %sampling at T & using thresholding operation
   z=sign(filtered);
   z=(z(T:T:end));

  %recovering the sent bits
  received_bits=(z+1)/2;

  %calculating the bit error
  Pb(i)=mean(abs((input)-(received_bits)));

 end

 %Plotting Pb Vs SNR
  figure(1);
 SNR=10.^(SNR_db/10);
 P_theory=(0.5)*erfc(sqrt(2.*SNR)./sqrt(2));
 semilogy(SNR_db,Pb,'o',SNR_db,P_theory,'r-')
 title ('SNR Vs Probabolity of Error');
 xlabel ('SNR (dB)');
 ylabel ('Probability of Error');
 legend('practical curve','theoratical curve');
 grid on;


What is the meaning of the sum of kernel components in digital image filtering?

Considering an LTI filter impulse response (aka filter kernel) $h[n_1,n_2]$, then its frequency response will be

$$ H(\omega_1,\omega_2) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j \omega_1 n_1}e^{-j \omega_2 n_2} $$

When an image is filtered, its average brightness is mostly affected by the DC response of $h[n_1,n_2]$, which is given by $H(0,0)$. When you want the brightness not to change after filtering, then you set $H(0,0) = 1$, which means that

$$ H(0,0) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j 0 n_1}e^{-j 0n_2} = \sum_{n_1} \sum_{n_2} h[n_1,n_2] = 1 $$

Which shows that the sum of impulse response (aka filter kernel) samples is one.


What image-processing techniques are ideal for this particular shift/scale invariant template matching?

Broadly speaking there are two primary approaches to solving this problem:

a. template matching or
b. matching with a parameterized model.

Personally, given the simplicity of the model, i would prefer the second approach for this problem.

Step 1: Identify the ping blobs

As a first step, extract the "Boxes" - essentially the yellow/blue squares. A simplest method here could be to just threshold the image. There doesn't seem to be much variations in the background except that as the noise increases it tends to be quite similar to the pings. Hence, the threshold can be global across the image - and i think you will be able to get a decent threshold that would work across the images. The threshold, however, should be smart so you can use something like Otsu's method. See this and this for more reference.

Improvement of the blocks

One of the good thing you can do to remove spurious points which looks similar to ping foreground where as some pixels inside the pings are also blue. You need here a morphological operation called "Opening". Here is one of the reference from HIPR. These type of operations requires smart shape that helps preserving similar shapes and removing others. In your case you can choose horizontal lines.  By the end of this, you have background and foreground pixels neatly separating the pings without touching each other.

Step 2: Mark/Register blocks

Now that you have got the cleanest possible blobs, time to mark them as 1,2,3,4 or A,B,C,D etc. This is usually accomplished by what is called some simple algorithm. Run through each pixel and if it is touching with any marked region, and is also foreground then assign it to that segment else create a new one. If you happen to find more than 4 important segments you need to choose the most relevant ones. Apply some heuristic here rather than any theoretical algorithm.

Step 3: Modeling of parameters

Now, once we get most likely top 4 candidates of the 4 corresponding pings, you can identify the bounding boxes - essentially the top left and bottom right corners (or left and right most pixels, as well as top and bottom rows of the blob). Find the bounding box, you can fined centroid the of each ping box.

Hence you can have four centroids $C_i \text{ where } i \in \{1,2,3,4\}$
Based on the description, the absolute location of $C_1$ can be arbitrary. However, following to $C_1$, the other centroids should be within certain range,

hence inter centroid distances can now be calculated:

$$

C_{1i}[x,y] = \{ C_1[x] - C_i[x],  C_1[y] - C_i[y] \}

$$

The other critical parameter is the length of the pings - which is $L_i \text{ where } i \in \{1,2,3,4\}$

Essentially you have 8 parameter vecoter:

$$
1.    \{ C_1[x],C_1[y] \}  \\
2..4  C_{1i}[x,y]  \\
5..8  L_i  \\
$$

Step 4: Classification
In the last step, now we need to have classification whether there is a ping-pattern yes or no.

For something like this, you can have a nice Bayesian classifier framework.

There are two Hypothesis you have

  H0 : The ping pattern exist.
  H1 : The ping pattern doesn't exist

Now, i am putting only 1 parameter $C_{1,i}[x]$ case here but you can extend this.

First of you study a lot of images where you know they belong to the image. Now, you can plot various histograms or apply some method of density estimation method. Read Pattern Classification by Duda Hart.

For simplicity you found that in cases where pings are present, the density function of parameter $C_{1,i}[x]$ is Gaussian with $ \mu = 20, \sigma = .1$ The sigma being small here could imply that there cannot be so much variations from the original.

This is actually conditional probability of $C_{1,i}$ given $H0$ is true. So you have values $P[ C_{1,i}[x] | H0 ]$

Now, from the test data

$P[ H0 | C_{1,i}[x] =x_1]$ is = given the reading of $C_1[x]$ being $x_1$, what is the probability of $H0$ ?

Hence, you can compute

$$

P[ H0 \space | C_{1,i}[x] ] = { P[ C_{1,i}[x] | H0 ] * P[H0] \over {P[ C_{1,i}[x]]  } }

$$

You need to combine this expression for the vector expression and need to fill in a lot of effort for exactly putting the parameter estimation to be done. But i am leaving you with basic approach here.

Please do your own math, this will be lengthy, but should be still quite intuitive.

Once, you combine expression, if $P[H0 | \text{all parameters}] > 1/2 $ then the ping pattern exists.

EDIT
Since you are already defining the specific criteria of 110 +/- 5 pixels, the last step can be simpler. You probably don't quite need to calculate all these probabilities if your classification criteria is fixed.


Analog to digital conversion using Python

You see this because you're sampling too close to Nyquist rate, but do not worry, although it seems to be aliased however, your sampling system works well and the signal is not aliased. To be more sure, just take a look at sampled signal in frequency domain:
As you see, only 9th bin is active, all others are zero, so no aliasing is occurred in sampling.The code (MATLAB) for FFT calculation is:

Fs = 20;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = length(sampling_signal);  % Length of signal
t = (0:L-1)*T;                % Time vector
y = sampling_signal;     % Sinusoids plus quantization noise
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

Note! I've changed time_of_view to 10 secs, to avoid FFT leakage (the side lobes you see around 9th bin).

It is O.K to increase signal rate this close to Nyquist rate, as far as you're not intended to return the signal back to the analog domain,
 because then your DAC requires a very sharp Low-pass filter in output and its realization is too hard (if not impossible). This paper intuitively expalins sampling problems, Sampling: What Nyquist Didn't Say, and What to Do About It


filter bank of gabor filters with different frequencies and orientations

Look on wikipedia. https://en.wikipedia.org/wiki/Gabor_filter The $\theta$ controls the orientation.


Apply decay envelope to waveform without multiplication

There are other methods if you can't support direct hardware multiplication, but you're going to have to sacrifice something. For instance, if you have multiple clock cycles to spare, you could use an iterative approach. C-like pseudocode:

accum = 0
for (i = 0; i < LENGTH(envelope); ++i)
{
    if (((envelope >> i) & 0x1) == 0) accum += signal >> (i + 1)
}
// "accum" now holds the multiplication result

Translating the above functionality to Verilog is an exercise left to the interested reader (it's been a while since I've written in that language). This scheme essentially performs long base-2 multiplication using additions and shifts. For each bit that is set in the envelope scaling factor, an appropriately-right-shifted version of the signal sample is added to the accumulator. The result of the multiplication is held in the accumulator after the loop is complete.

As I said, this isn't the most efficient approach in the world; it's probable that a good synthesis tool would come up with something better even on a device without hardware multipliers. This method trades resource usage for latency in the form of extra clock cycles required to perform the iterations of the unrolled multiply. I'm left to wonder what drives the no-multiplication constraint in your application. Could you add more details on that?

Note: You may need to tweak the amount of extra shift that you apply to the signal when adding to the accumulator based on the numeric format of accum and signal. If they are both of the same fixed-point format (like Q15), then the above approach should work properly.


Implement a software PID for the first time in real time software, I can find P what about I and D?

Once you get the basics covered (Jason has covered this pretty well), you will also want to research Integral Windup and filtering of the error signal and its estimated derivative.  The former is very important if your control point changes discontinuously (i.e. you don't have a separate controller 'ramping' the control signal for your PID.)  The latter is important in the presence of real-world sensor noise, in order to avoid instability.  Many sensors can introduce large noise 'spikes' that can greatly perturb a control system without filtering.

Many physical systems also exhibit resonance effects that can lead to run-away oscillations (the Tacoma Narrows Bridge is the classic example of how bad this can get.)  A feedback system can exacerbate these effects, so an important part of 'tuning' your PID is characterizing its stability at the various feedback and control signal frequencies you will encounter.  This is another area where filtering of the feedback signal is helpful; in particular, notch filters at 60Hz are common for systems using AC power in the US.  In practice, the resonances are application-specific, so most controller manufacturers supply a general toolkit of filters that can used and then leave it up to the process engineers to configure them appropriately for their application.

Next steps after you've covered stability:

'Ramping' the control signal by applying constraints on its rate of change.  For motion systems, this typically involves limiting velocity and acceleration.  This makes the job of the PID easier, as the control signal is moving more 'smoothly'.  It also generally results in a smoother output signal for your actuator, which can be important, since most actuators have practical limits on their ability to respond to rapid changes in their control signal.
Adjusting the feedback/control signals outside the PID to get the best response from your particular system.  For instance, there may be known hysteresis effects in the actuator or sensor that you want to correct for.

All that said, there's a reason why most integrators buy off-the-shelf controllers for the devices they're integrating.  Control systems are complicated and there are a lot of proprietary tricks that vendors have developed over the years to get the best performance out of their particular components.  It's still valuable to understand these systems to get the best results from them, but there is more going on here than you may realize.


If your environment's noise level allows it, you can use a quiet no-brainer solution. Just make sure that the impacts have greater peak amplitude than the noise floor in every case.

Specify a threshold level based on previous experiments. This level has to be great enough for all the impacts you generate and you want to detect. If your signal crosses the level, you can assume, that an impact has happened.
If a detection has occurred, you can start a timer, that prevents unwanted re-detection for the same impact because you can't expect that your impact will be below the threshold for the next sound sample. For this reason, the timer has to run longer than the maximum detectable length of any impacts. Until the timer runs, you disable the threshold detection.
After the timer has fired, you can enable the threshold detection again.

This method requires previous experiment and knowledge about the detectable impacts. You have to determine the minimum sampling rate that allows you to catch at least a few sample of the impact. Standard audio quality sampling rate would be good enough.

You can of course apply signal conditioning by filtering the signal with a digital filter, but this necessity is highly depends on your environment and noise floor.

Note that this is possibly the most no-brainer solution. There are many more sophisticated methods for this particular problem.


Calculate amplitude in dB and dBA from audio signal

Q1: Is this the correct way (or a valid way) to get a dB amplitude from a sample?

Assuming you want to build a sound pressure level meter the answer would be a resounding "no". The correct way is to build a running energy detector with a proper time management as defined in the IEC standard IEC61672. https://webstore.iec.ch/publication/5708. This is actually pretty simple. By far the most difficult part is the calibration of the microphone. As ZR Han in their answer already pointed out: you get to a sound pressure level in dB (dBSPL) by using a reference pressure of $20 \mu Pa$, but that means you have to convert your input signal from Volts or "number" to Pascals as well and this is very device and settings dependent.

Q2: How do I interpret this bin?

That's a complicated question. Each digital bin contains energy from ALL analog frequencies EXCEPT the center frequencies of the other digital bins and the exact relationship isn't trivial. In your case you should just ignore this since there is no reason to do go into the frequency domain in the first place.

Q3: I have 256 bins, ...

Again, that's a very complicated question and I won't get into the answer since its completely unnecessary for what you want to do.
Building a sound level meter is relatively straight forward

Calibrate your input signal to represent actual sound pressure in Pa (Pascal)
Apply time domain weighting filter (A,C, Z, none)
Square the signal
Apply exponential (first order lowpass) averaging with proper time constant
Convert to dB using $p_{ref} = 20 \mu Pa$

Step 1 is by far the most complicated.


By definition, the magnitude of the analytic signal $s_a(t)$ can always be considered the envelope of the signal in a purely mathematical sense, simply because the real-valued signal is given by
$$s(t)=\textrm{Re}\big\{s_a(t)\big\}=|s_a(t)|\cos\big(\arg\{s_a(t)\}\big)\tag{1}$$
However, this does not mean that $|s_a(t)|$ has necessarily a physical significance as an actual envelope of $s(t)$. It does have physical significance if the term $\cos(\arg\{s_a(t)\})$ in $(1)$ has the function of a carrier. This carrier doesn't need to be sinusoidal but it must be a narrowband signal, i.e., its center frequency must be large compared to its bandwidth.
If you want a strict mathematical definition of the conditions under which $|s_a(t)|$ is an envelope in a physically significant way, then you must define what you mean by 'envelope' and by 'physical significance'. This is not trivial and probably also not very useful.
If your question is just about the feasibility of envelope detection, then that's a very different question, and one that's much easier to answer. If the modulated signal is given by
$$s(t)=m(t)\cos(\omega_ct)\tag{2}$$
then the magnitude of the corresponding analytic signal equals $m(t)$ if

$m(t)$ has a maximum frequency smaller than the carrier frequency $\omega_c$
$m(t)\ge 0$

These conditions are met for standard AM, but not for DSB-SC.

EDIT: For the signal $(2)$ it's straightforward to show that $|s_a(t)|=m(t)$ if above conditions are met. If $m(t)$ is a low pass signal with a maximum frequency smaller than the carrier frequency $\omega_c$, then the analytic signal is given by
$$s_a(t)=m(t)e^{j\omega_ct}\tag{3}$$
with magnitude
$$|s_a(t)|=|m(t)|\tag{4}$$
Clearly, if $m(t)\ge 0$ then $|s_a(t)|=m(t)$.


Mapping Saw Wave That Is Between -1 And 1 To A Sine Wave

You can scale your sawtooth before you put it into the sine function, like so:

f=(-1:1/128:1)';

ramp = f * (length(f)-1)/length(f);  % Scale down to not include the last 1 into the sine
c1=repmat(ramp,1000,1);

c2=sin(pi*c1);   % -1 of saw wave maps to 0, 1 maps also to 0

plot([c1 c2]);
%sound(c2, 44100);


Estimating Velocity from Fuzzy Position Data with Known Uncertainty

Yes it is. That what's Kalman filter is all about. Your measurement includes position data, and there is a hidden velocity parameter. You will need some kind of model, like constant-velocity. You can input your idea of uncertainty into Kalman noise matrixes to improve results.
You can use Kalman filter, which is good to predict the future. If you just want to estimate the velocity in post-processing step (after you have all data), you can use something simpler, like estimating the derivative. This will have the advantage of being smoother and more accurate.


Method for calculating musical note-based magnitude spectrum

Have you looked into a Chroma/pitch-class analysis? That sounds more like what you want overall (expect perhaps for the octave-folding/collapsing).

Check out librosa's CQT-based chromogram. Code to generate plot below on github:

You can use the tuning parameter to tune the musical scale, but study the full paramter set.

tuning : float
    Deviation (in cents) from A440 tuning

Finally, checkout these examples showing "enhanced chroma", to say isolate harmonic content:

https://librosa.github.io/librosa_gallery/auto_examples/plot_chroma.html#sphx-glr-auto-examples-plot-chroma-py

Regarding performance of the CQT computation, the FFT + a kernel can be used to make the computation much faster. See Judith Brown & Miller "Super" Puckette's An efficient algorithm for the calculation of a constant Q transform.


Spectral peak location estimation using complex DFT

I admit that I failed my first attempt.
The text suggests that we would take $g(m) = am^2 + b m + c$, interpolating $g(k) = X(k+l)$, for $k \in \{-1, 0, 1\}$, interpolating this would give us polynomial parameterized in $X[k], X[k-1], X[k+1]$, derive, solve and it gives the optimum value at $m = k+\delta$. Not that easy, and it gave me three roots, none of them seemed to be the desired value :(
The comment on the paper is not helpful but they give us a reference [4], from the one of the authors, and there it gives more details. There you will find the equation (3.7) and they give a concise derivation on the Appendix A.
Derivation
Here I will fill some details of the proof.
You express the interpolator polynomial centered at $k + \delta$, so delta will maximize it by construction (more clever the attempt described above).
$$g(m) = a (m - k - \delta)^2 + b$$
Then $g(m)$, determine delta so that interpolates the three samples. Since the quadratic polynomial interpolating three points is unique, we know that $\delta$ here is unique as well.
$$\begin{eqnarray}
X[k-1] &=& g(k-1) &=& b+((k-1) - k - \delta)^2 &=& b + a\cdot(\delta^2+2\delta + 1) \\
X[k  ] &=& g(k  ) &=& b+((k  ) - k - \delta)^2 &=& b + a\cdot(\delta^2) \\
X[k+1] &=& g(k+1) &=& b+((k+1) - k - \delta)^2 &=& b + a\cdot(\delta^2-2\delta + 1)
\end{eqnarray}$$
Start looking at the coefficient of the rightmost terms, and note that
$$2 X[k] - X[k-1] - X[k+1] = 2 a$$
Also note that
$$ X[k-1] - X[k+1] = 4 a \delta$$
$$\delta = \frac{X[k-1] - X[k+1]}{4a} = \frac{X[k-1] - X[k+1]}{2(2X[k] - X[k-1] - X[k+1])}$$
With the sentence

A simple adaptation provides an estimator suitable for estimating tone frequencies from complex DFT outputs.
They introduce the formula we see in the first paper you mentioned.

The adaptation they mention is because $\delta$ is complex, so he takes the real value. Maybe a justification for that would be that $\min_{f} |f - \delta|$ is $Re[\delta]$.


How do I calculate the noise variance in this problem?

Since the signal $r(t)$ is expressed as the low-pass equivalent of a bandpass signal, in general it will be complex. That means that the noise is also complex.

White Gaussian noise has a power spectral density commonly denoted by $N_0/2$. This means that if you filter the noise with an ideal filter of bandwidth $B$, the noise at the filter's output is $\sigma^2=BN_0$.

Note that this is the power "per dimension". That is, the in-phase and the quadrature components of the low-pass signal $r(t)$ each suffer from noise of power $\sigma^2$.

You define the SNR as $E|\alpha|/\sigma^2$ and you need it to be equal to 10. One way to do this is to configure your simulation so that $E|\alpha|=1$, and then choose $\sigma=\sqrt{0.1}$. If you later need a different SNR, you can just change one value ($\sigma$) in your simulation.

You can generate the complex noise as follows:

N = length(s);
noise = sqrt(0.1)*(randn(1,N)+1i*randn(1,N));
r = alpha*s + noise;

Note that in this setup, the total noise power is $2\sigma^2$. The in-phase noise and the quadrature noise each have power $\sigma^2$.


How to build a simple model for voice-controlled system?

DTW does not really work in speaker-independent systems, you'd better use a proper speech recognition toolkit.

For 10 recordings of a single speaker, you compare to all one by one and select the best match.


Given some parameters, how do you determine how much of an image sensor an object will consume?

The basic idea is that the ratio of the distance/focal length will be the same as the physical/image ratio.

So, assuming a focal length of 100mm, this gives:

$mag=10m/100mm=100$
$size_{image}=1m/mag=10mm$

Note, the distance is the distance measured from the focal point, which will be 1 focal length from the sensor, so that might factor into the equation as well.

See this calculator to play around with things a bit more.


How to implement path loss in wireless communications?

A good text is by WIlliam Y. C. Lee. He did alot of the original work on
propagation models. There are many others, just google "path loss model".
This  equation you selected is a very basic one for a cellular network. It doesn't take into account many of the fActors like environment(urban, suburban, rural), clutter or morphology, line of sight propagation or non line of sight, statistical factors and others. There are empirical models for various environments that have been used and modified (Hata, Okumura cost 231 etc) and can be iteritively adjusted with real field data.

With regard to your question about the reference distance. I'm sure you can find the answer in the references I mentioned above and the Goldsmith reference mentioned by ...

There might be a relationship to that distance and the far field or it may be justa standard reference, I'm not sure. In any case you could use a different model if the user is too close if necessary.

The path loss exponent is also a very important parameter(alpha in your equation)

Hope this helps


Alignment of 2 Set of Samples from Different Sensors

There are 2 main factors here:

Aligned timeline
Namely each sample has time stamp taken with the same reference time.
The lower sampling rate is high enough according to the sampling theorem
Namely even the lowest sampling rate, in your case 32 [Hz] is high enough to reconstruct the original signal.

If both (1) and (2) hold you'll be able easily to compare them by interpolating both into a shared time grid.
If (1) doesn't hold you may be able, by interpolation and some other DSP techniques (Such as Matched Filter and Dynamic Time Warping), to align both to the same grid.
Yet it would be harder (And in some cases, not achievable).


The easiest way to do this is to keep track of the parameters a and h from detection. How did you pick out a car from the scene? a simplified method would be: background subtraction, followed by clustering, data association, find centroid of clusters => x,y. Also store the bounding box of the cluster, easiest way to keep an accurate a,h.
But I'm guessing you want to start with initial conditions a0,h0 and based on x,y,u,v update a,h. The problem is that a,h cannot be completely known using only x,y,u,v,T. You also need information about the imaging system or the environment.
Assuming we use a Cartesian coordinate system x=horizontal, and y=vertical
there are a few different ways to view this problem. If you are tracking cars from a bird's eye view (think google maps for instance)

then for all cases

a = a'
h = h'

realistically there can be slight variations. for instance if the car is facing this way __ vs this way | the bounding box would change, but given a set of x,y,u,v we have no way of knowing the cars orientation. Concretely, if the car is headed in the same direction (red and green) the bounding box shouldn't change. If the car changes heading (blue car) the bounding box changes drastically,

h = width' and width = h'

recall

a = w /h

so

h = a'*h'
a = width / h = h' / (a'*h')

unfortunately, we have no way of knowing with absolute certainty that the car changed directions. We can make up some kind though. Since cars are non holonomic vehicles ( a car can move forwards and backwards instantaneously but not left and right instantaneously) we can come up with some criteria to figure out a cars direction. For the red car we see |v| > |u| and |v'| > |u'|, or for the green car |u| > |v|  and |u'| > |v'| but for the blue car |v| > |u| and  |v| < |u| in this case we could guess that the car changed directions and the bounding box needs to be rotated using our second set of equations. However, if u~=v all bet's are off and we can't safely say which direction the car is pointing.
More complicated still is an image with perspective. For the following examples please note the origins in my picture are wrong, we should use the standard image convention where (0,0) is the upper left corner and (m,n) is the lower right corner where m=number of rows and n = number of columns. For this example I will use a image with a single vanishing point outside of the image.  We can split this into a few cases
CASE 1: car is moves horizontally along x

assuming the field of view is small we get

a = a'
h = h'

why? the height will only change as an object gets closer or further from us, since it is moving only along x, v=0 so our height shouldn't change. our width also shouldn't change (again the viewing angle must be small) because if we have a car traveling at velocity u, all points on the car move the x direction with the same u. More concretely, if delta_x = (x(T) - x(T-1)), delta_x should be the same for all points on the car. If the front tire moves delta_x (orange in pic), so should the back tire(light blue in pic), otherwise the car isn't moving as one solid body. IF both points moved the same delta_x, they should still be in the same position relative to one another (purple line in pic) so that means width doesn't change. we already said h doesn't change if the car moves only along x, and not y, so a shouldn't change either
CASE 2: car moves vertically, along y

this is a bit more difficult. As far as I know there is no way to determine this without performing some kind of offline calibration. There are a few issues, first is that size alone can't determine distance in an image. Think of an image with a person, and mountains in the background. The person looks much larger than the mountains, but we know this isn't true, however, a computer would have more difficulty figuring this out.
still using the single vanishing point outside of the image. We need to know something about the image, lets say we know the street is 6 meters  wide (call this Xstreet). Now we measure the street at various y locations. at the top of the image Xstreet0 = 100 pixels. At the bottom Xstreet1 = 180 pixels. Lets also say that Xstreet0 and Xstreet1 are 100 pixels apart vertically. Even though the width of the street appears to be changing we know in reality the street is really 6m at both places. So we come up with some kind of equation to mirror our perspective. Remember this equation

y = Mx + b

lets use that, so to find our slope (we aren't really interested in b), which models our perspective

M = rise / run = (y1,y0) / (x1 - x0) = 100 / 80 = 1.25

that means for every 1 pixel we move downwards in the y direction our x pixel gets 1.25 times bigger! With this knowledge we can now do some modeling

h ~= h' (this is really difficult to determine, it gets into perspective and turns very messy, since we said there is only 1 vanishing point outside of the image we can  say this doens't change much)
width = ((v*T)/M) + width'

looking closer we see that

v*T = delta_y
delta_y = M * delta_x ==> delta_x = delta_y / M ==> delta_x = (v*T)/M

so the term

(v*T)/M

gives us the change in X (delta_x), we simply add this to the old X to get our new X. but we aren't done yet a = width/height not simply width! so lets keep going

h = h' + delta_h
width = width' + delta_width
a = width / height = (width' + delta_width) / (h' + delta_h')

recall h ~= h' so delta_h = 0

a = (width' + delta_width) / h' = width'/ h' + delta_width / h' = a' + delta_width / h'

well we just calculated delta_width == delta_x, we know h' so we have our final equation

a = a' + (v*T)/(M *h')

where M was our slope (perspective compensation) we calculated offline.
CASE 3: a movement in x,y
is really just a movement in y, we have already said any change in x doesn't change the parameters at all. only a change in y

h = h'
a = a' + (v*T)/(M *h')

This problem gets a little more complicated when the vanishing line is actually in the scene. In this case, we have to keep track of if we are below or above the vanishing point and signs become very important.
This method also assumes that the objects we are tracking are in a constant z-plane. If our scene has a 3d axis, in our example, the street would be at one z value and cars travel in the x,y plane. We assume the cars don't levitate up and down, or our assumptions about y movements become invalid.
I mentioned it briefly, but this also won't work with wide FOV (field of view)cameras because our assumption that moving along x doesn't change width or height becomes invalid. If this doesn't make sense, just watch some fisheye video clips, like this video around 0 minutes 50 seconds
conclusion
Overall, my methods are highly approximate, and only work under certain conditions (assuming a camera doesn't move often, they are fairly reasonable restrictions). The absolute best way to keep an accurate a,h would be to save the bounding box from the original car classification steps. Hope this helps someone in the future


Time series forecasting on dsp embedded systems

I'm not sure why you oppose statistical signal processing to digital signal processing, they are not mutually exclusive. If you perform statistical signal processing with non-continuous quantized samples, you perform digital signal processing.

Second of all, if you are using a microcontroller, you are likely to have limited RAM resources. Can you clarify what kind of algorithms you are using? You need to estimate the RAM resources required for your algorithms.


Sine wave in frequency domain on signal analyzer

in frequency domain each spike/impulse represents to a sinusoid in time domain. so if in frequency domain the signal shows a single spike/impulse it means the input signal has a sinusoid in time domain.


Calculate the Derivative of Gradient Field of an Image

Consider the expansion of the term below
$$\|\nabla x -p\|_2^2 = (\nabla x -p)^T(\nabla x -p)$$
$$\|\nabla x -p\|_2^2 = (x^T\nabla^T -p^T)(\nabla x -p)$$
 $$\|\nabla x -p\|_2^2 = (x^T\nabla^T\nabla x - x^T\nabla^T p  - p^T\nabla x +p^Tp)$$

Now consider the following basic definition:

$$\frac{\partial(A x)}{\partial x} = A^T$$
$$\frac{\partial(x^TA)}{\partial x} = A$$

Now applying the above two definitions together with the expansion of the objective above to differentiate the objective we have

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla^T\nabla x - 2\nabla^Tp$$

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla^T(\nabla x - p)$$

since $\nabla^T = \nabla$, therefore we have

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla(\nabla x - p)$$ the constant 2 is just a constant, so the result we have is consistent with the ones that authors are using, its simply a consequence of vector differentiation


Allpass filter feedforward/feedback design and code

First of all make sure that the magnitude of $a_M$ is strictly less than $1$; this is necessary for the filter's stability. Second, you use several unnecessary variables. This is no serious problem, but it makes it harder to see what's going on. For $M=1$ a simple pseudo-code would be

v = x - aM * v_old;
y = v * b + v_old;
v_old = v;

The important thing is that the state variable v_old must be remembered between function calls. Are you sure that your variables allPassFeedback and allPassDelayedSample are not reset to $0$ every time you call the function?


The magnitude of the frequency response will remain unchanged if you reflect any poles outside the unit circle - these are the ones causing instability - back inside the circle. I.e., a pole $p$ (with $|p|>1$) is replaced by the new pole $\tilde{p}=1/p^*$, where $*$ denotes complex conjugation. This will not change the magnitude of the frequency response (up to a scaling factor), but it will change the phase response.

The correct scaling is obtained if each factor $(1-p_iz^{-1})$ with $|p_i|>1$ of the denominator polynomial $A(z)$ is replaced by $|p_i|(1-z^{-1}/p_i^*)$.


It would help if you could state what exactly you are trying to do.
The "relative loudness of the harmonics" is normally a function of time because of the different decay rates and/or temporal envelope of the harmonics. There are different ways to condense this into a single number: initial amplitude, peak amplitude, total energy, 250ms into the decay, etc. The most suitable method would depend on your specific goal.

The look up table looks questionable to me. None of the loudness curves that I'm aware of (A-weighting, C-weighting, Fletcher Munson, Robinson Dadson, ISO 226) look like that. -9dB at 1500 Hz seems just plain wrong.

The shape of the loudness curves is a function of absolute level. So whatever curve you pick it can only be correct at one level (or a small range of levels). In order to get this really right, you would need measurements that are calibrated in Pascal or dBSPL.


how to extract 45 cycles from bunch of waveforms

For best performance in the presence of noise, I recommend "locking" onto the cyclical signal using a Phase-Lock-Loop (PLL) approach. Then you can feed the PLL output (which will be a cleaner representation of your signal and can be a digital square wave if desired) into a digital counter, which you can then count to 45 after resetting the counter at any position desired.


Creating uneven sidebands with AM + PM modulation?

If I understand correctly, you want to use each side band as an independent information-bearing signal. In other words, you want to be able to generate $$a_k e^{2\pi f_i t}$$ with the information in carried by the amplitude $a_k$. There are several ways to do this, but I don't know which are feasible using optical processing. All of the following methods are feasible if you start with discrete signals and then convert them to optical.

Single-sideband modulation. The idea is that $$s(t)\cos(2\pi f_i t) \pm j\hat{s}(t)\sin(2\pi f_i t),$$ where $\hat{s}(t)$ is the Hilbert transform of $s(t)$, is a single side band signal or, in other words, one of the two sidebands you want to create. By choosing $+$ or $-$ you select the lower or the upper sideband.
Use the IFFT to go from any given spectrum to the time domain. This is the way OFDM works.
Use quadrature modulation (QAM). With this technique, you don't get to control each sideband independently; rather, you transmit two independent bitstreams by mixing two different versions of the sidebands. I have heard that QAM is possible to implement optically.

I hope this answer points you in some fruitful direction.


Tracking movement of audio speaker driver using images

A single frequency coming out of speaker is very much a stationary process, then why go for a Kalman filter? Use a Weiner filter instead, much more easier to setup and makes sense for your problem.

Just go for a standard weiner filter, you could go for recursive least squares as well to model this.


For what you intend to do, a low-pass filter is the way to go. Your statement about filtering frequencies vs filtering amplitudes is incorrect. Your signal contains many components at many frequencies, the amplitude of which varies in time, and the high frequency components are those causing the "jaggedness" and you want to get rid of them. Not sure why you say your signal is "constant frequency" - maybe you are getting confused by the sample rate?

What you have tried (averaging) is indeed a special case of low-pass filtering, but one with a frequency response far from being ideal. You should try a properly designed IIR or FIR filter. In particular, the FIR filter is not very different from what you tried - this is just a weighted combination of the samples neighboring each sample. But the choice of the coefficients is important and ensures that only unwanted components are eliminated. Note that an FFT is not the way to go. This question crops up here quite frequently under different forms, but in short - FFT, messing with coefficients, IFFT - is a bad idea.

By design, the output of a moving average filter (what you implemented) has less energy than the input. It is thus impossible for a moving average filter to cause distortion. If the input signal is in the [-1, 1] range, there is no way for an averaging filter to yield values outside this range. The "distorted sound" you observed was probably due to an implementation error on your side (overflow / clipping of integer values, signed values treated as unsigned value, or maybe incorrect in-place processing)...

EDIT: one thing worth mentioning is that there are situations in which a speech signal actually has high frequency components (appears "jagged", is noisy) - for example during a sss or shhh ; and removing those with a low-pass filter will affect its brilliance. Ideally, you'll want your low-pass filter to be active only when you detect that your speech signal is voiced - and inhibit it when you detect an unvoiced, noisy consonant.


For initializing window you need the start point would depend on your overlap. If the overlap is 50% then you would start at window_length/2
If you are using a rectangular window and with 50% overlap in that case you want to use a square root rectangular window on the input and output so that there are no scaling issues. so basically you window will be scaled also dependent on your overlap.

Also it would be worthwhile to consider using other windows dependent on your application. As in a lot of other windows the scaling is taking into account.


What you are seeing are indeed numerical artefacts. The conversion from polynomial roots to polynomial coefficients (and vice-versa) is known to be an ill-conditioned problem. In your case the problem occurs at a relatively low polynomial order because zeros and poles are all clustered around the same point (DC), and the poles and zeros almost cancel each other. By removing one zero, as you tried, the problem's numerical condition is improved, simply because the polynomial order becomes smaller. Another way to avoid the problem is to distribute the poles and zeros more evenly around the circle in the complex plane. Note that the available frequencies extend up to $f_2/2=22050\text{ Hz}.$ That means that just decreasing the sampling frequency would also improve the numerical condition.

Apart from the numerical effects, your approach for converting zeros and poles to numerator and denominator coefficients is correct. However, it can be done in a much simpler way by using the function zp2tf, which transforms given poles and zeros to transfer function coefficients.

Also note that with your approach you generate a double zero at zero frequency because you produce complex conjugate pairs, even for real-valued zeros.

Finally, the filter you computed is not causal and stable. This has nothing to do with the effects you see, but I believe it's still worth mentioning. So you wouldn't be able to implement that filter with the coefficients you computed. The reason for this is that the transfer function is not proper, i.e., its numerator order is greater than the order of its denominator. This means that apart from the given poles, there occur a number poles at infinity, and that number equals the difference between the numerator and denominator orders. This problem can be solved by adding poles at $z=0$, i.e., by multiplying the denominator by $z^K$ where $K$ is the difference between the polynomial orders of numerator and denominator. This will not change the frequency response, it will just make the filter causal and stable.

EDIT:

For maximum precision you shouldn't compute the coefficients of the numerator and denominator polynomials. Instead you should work with second-order sections. Below is a simple example in Octave/Matlab. Note that I didn't optimize the sequence of sections, neither did I optimize the pole-zero pairing per section. But, as you can see in the figure, the numerical artefacts are greatly reduced.

fs = 44100;
rz = 1; rp = 0.9999;
phi_z = 2*pi*[30; 45.8; 66] / fs;
phi_p = 2*pi*[41.7; 60; 95.2 ] / fs;

m = length(phi_p);
L = 2^15;

a = zeros(m,3); b = zeros(m,3);
H = zeros(m,L);

for k = 1:m,
    b(k,:) = [1, -2*rz*cos(phi_z(k)), rz^2];
    a(k,:) = [1, -2*rz*cos(phi_p(k)), rp^2];
    Htmp = freqz( b(k,:), a(k,:), L );
    H(k,:) = Htmp(:).';
end

% add real-valued zero at DC
[H0,w] = freqz( [1, -rz], 1, L );
H = [ H0(:).'; H ];

Htot = prod(H);
semilogx( w/pi+1e-5, abs(Htot) );
axis([1e-5,.01,0,.0005])


What is the threshold to prefer polyphase channelizer over FIR filtering?

A polyphase channelizer is not a special kind of filter. It is a structure that works well when using filters in multi rate settings. Polyphase is a sampling rate conversion method that leads to efficient implementations that are useful for building filter banks. The efficiency comes from only having to design one filter. The downside is that the extracted channels have to be uniformly allocated and the signal sample rate has to be an integer multiple of the channel's bandwidth.
Use a channelizer for 2 more channels that are uniformly allocated.


Sine wave amplitude modulation - Analysis using Audacity

Your samples have a finite bitwidth. For example, assume your samples have 16 bit width; if you increase their amplitude beyond what these 16 bit can represent, you end up with an overflow.

In case your audio samples being signed integers, that overflow typically first manifests as a sign inversion – exactly what you're seeing here. (That's because only the lower 15 bits in your 16 bit signed integer are actually the number – the uppest bit is the sign; you shift a 1 there, the result becomes a negative number in 2-complement represenation.)

This isn't a limitation of audio processing only – in fact, integer limits is something basically (nearly) every programmer needs to be (slightly) aware of. Your number variables can't hold arbitrarily large values!


Training Fractionally Spaced Channel Estimator and Equalizer

The training sequence should be at the same spacing as the equalizer when considering its sampling at the input to the equalizer. Adaptive algorithms converge to the least square solution based on the error between the received sequence and the transmitted sequence (known when a training sequence is used). Further, the equalizer can only determine a solution where channel energy is present (for example, if you transmitted a single tone, you would only find the single phase/magnitude weight for that frequency)

If you use only one sample per symbol to derive the channel, the solution will suffer from "band edge aliasing" as I depict in the figures below, and the channel solution will suffer from the folded energy at the band edges.

For further details see this post where I have derived the operation of the LMS Equalizer and the solution using the Wiener-Hopf equations (for when you have the luxury of block post-processing, but also helps illuminate the operation since this is the same solution that the adaptive algorithms would converge to).

Compensating Loudspeaker frequency response in an audio signal


For the CWT the values in your definition can be any $a\in \mathbb R^+$ and $b\in \mathbb R$ if applied to functions. The CWT can also be applied to discrete signals (discretized CWT) with $b$ restricted to appropriate discrete values (and any $a$, but typically $2^{j/v}$ where $v$ is the "voices/octave" of the transformation). Now when we restrict $a$ to some discrete values $2^j$ we get the DWT, which is an orthogonal transformation.

(Small edit to correct the scale : $2^{j/v}$ instead of $a^{j/v}$)


How to compute transfer function $G(s) = \exp \left( - \sqrt{s} \right)$ in Matlab / Simulink?

The inverse Laplace transform of $G(s)=e^{-\sqrt{s}}$ can be computed analytically:
$$g(t)=\mathcal{L}^{-1}\left\{e^{-\sqrt{s}}\right\}=\frac{1}{2t\sqrt{\pi t}}e^{-\frac{1}{4t}},\qquad t>0\tag{1}$$
[I've cross-checked this result by looking it up in one of my old math books.]
So now you have an analytical expression for the impulse response $g(t)$.
The step response can also be computed analytically. It is given by
$$s(t)=u(t)\int_0^t g(\tau)d\tau=\textrm{erfc}\left(\frac{1}{2\sqrt{t}}\right)u(t)\tag{2}$$
where $u(t)$ is the unit step function, and $\textrm{erfc}(\cdot)$ is the complementary error function.
The result $(2)$ was obtained with a little help from WolframAlpha.
And that's how the impulse response and step response look like:


Frequency response of numerical differentiation by polynomial interpolation / finite difference

Polynomials are not integrable, hence their Fourier interpretation is complicated. However, locally, there are a few works. Under least-squares, polynomial interpolation's linearity works well with frequency analysis.

The class of Savitsky-Golay filters performs least-squares polynomial interpolation, and yields different derivative orders. As they seemed to be little known in signal processing, R. W. Schafer dedicated them What Is a Savitzky-Golay Filter? (public version) in the Lecture Notes section of IEEE Signal Processing Magazine in 2011. There, he especially covers their properties, including frequency responses. Analysis was also performed in Smoothing by Savitzky-Golay and Legendre filters by P. O. Persson and G. Strang in 2003. One can combine derivative accuracy and noise shaping or frequency response fitting.

More recently, there were derivations with IIR filters and fractional differentiation, for instance Maximally Flat IIR Smoothers With Repeated Poles and a Prescribed Delay, and even non-linear versions.

In image processing, there has been a lot of works as well, I can dig into my reference database if needed.


Unit impulse vs Kronecker delta vs Dirac delta?

Kronecker Delta used mainly in context of discrete time signals is defined as $$\delta_{ij} $$ is 1 only when $i=j$ or $i-j=0$, $0$ otherwise. So alternatively it is also written as $\delta[i-j]$. This value is exactly 1 and finite (refer below for the distinction with Dirac Delta)

The unit impulse is just the kronecker Delta with $j=0$ hence we only refer to unit impulse with one parameter $\delta_i$. Since $j=0$, this is alternatively written as $\delta[i ]$.Hence is 1 at $i=0$, $0$ otherwise.

The continuous time unit impulse can be thought of derivative of continuous time unit step (modified to have a slightly less steep than the sudden step at $t=0$). Refer to Oppenheim's book on signals and systems Chapter 1.

The Dirac Delta is used in context of continuous time signals to define impulses. It is defined as having an infinitely small width and hence a large magnitude, however the area under the curve integrates to 1. Refer to Oppenheim's book on signal and systems for various shapes that a Dirac Delta function could take, in all cases the area being 1.


FFT giving a huge magnitude of first frequency and pretty much zero after that

That's correct. It's simply a property of your signal: your signal has very small "wiggles" with a very large mean value. Your mean value is about 1 and your wiggles are around 0.01 we would expect the first frequency to be about 100 times larger then all the frequencies.

The first frequency represents the DC component (at 0 Hz), it's simply the sum of the signal.


How to determine the frequency of each row in a MATLAB spectrogram

The physical frequency $f$ that a given Discrete Fourier Transform (DFT) frequency bin $k$ corresponds to is $f = \frac{k}{N_{DFT}} \cdot Fs$.

For an $Fs=44100$ Hz, the middle point of a given $N_{DFT}$ point DFT would correspond to $\frac{1}{2} \cdot Fs$, also known as the Nyquist Frequency.

EXCEPT MATLAB, which maps $1$ to $\frac{Fs}{2}$.

For more information, please see here.

Hope this helps.


The processor is stalled. For example you are trying to read a peripheral register and there is no data available. This can have many different reasons so without further detail that's impossible to assess. First step would be to look at program memory at address 0x20a5a and see  what the processor was trying to do. Take a look at the disassembly window.


What are the mathematical consequences of a discrete signal living in Hilbert space?

whether or not a signal lives in a Hilbert space depends only on if you say it does.  if you say so, then you need to specify an inner product $<x,y>$ (which is pretty easy for finite-length discrete-time signals)

$$ \langle x,y \rangle = \sum_n x[n] y^*[n] $$

and you need to use that inner product definition to define a norm

$$ |x| = \sqrt{\langle x,x \rangle} $$

which should also be the same as the distance measure from $x$ to the $0$ element in that Hilbert space.

and you need to be able to add and subtract elements in the Hilbert space (easy to do for discrete-time signals).  why this might be useful would be for defining a signal $x$, an estimation of that signal $\hat{x}$, and an error signal, $e_x$ such that

$$ \hat{x} = x + e_x $$
or
$$ e_x = \hat{x} - x $$

we can think of the norm of $e_x$ which is

$$ |e_x| = \sqrt{\langle e_x,e_x \rangle} $$

as a measure of the distance that $\hat{x}$ is from $x$.

if gives you a formal framework for defining error and then maybe help you with algorithms (like a matched filter) for minimizing error, so that a received signal most closely approximates what was transmitted.  or so that some other  "thing" most closely hits some specified target "thing".


Why MATLAB only allows cpsd and mscohere between vectors with the same length?

The calculations are done in the frequency domain using pairs of FFT's which is why the vectors need to be the same length.  You could truncate the longer vector to the length of the shorter. If the longer is much longer, you could select multiple (different) truncated segments.

You can do things in the time domain like cross correlation, but mscohere and cpsd average (the default of 8 segments).  A single crosscorrelation isn't equivalent.

When reading "does it make sense" (you're not the only one who asks), Internally, I think this is a lack of confidence and I don't want to discourage this person from trying, but this is whacko .  Try it. MATLAB is easy.  It's designed to explore the whacko.


Why does alpha of 0.5 in detrended fluctuation analysis indicate randomness?

The application of the DFA method starts with a conversion of a bounded time series under test into an unbounded process, the elements of the unbounded process being partial sums of the bounded time series sequence.
Let we have a (discrete-time) realization of WGN of length N. It is a bounded time series. The starting conversion procedure of the DFA method transforms this WGN time series sequence into a realization of the random-walk process by the definition of the random walk process. For the RMS deviation's log-log plot of this random-walk process realization (this RMS deviation named the fluctuation in the DFA parlance), we draw a straight trend line using the least squares method. The slope of this trend line approaches 1/2 when N tends to infinity.
Hence the exact equality of DFA's alpha to 1/2 for the random process. You may have failed to find an explanation in scientific literature, because it should be read from textbooks, for example, read the Wikipedia article on DFA.
UPDATE
For an uncorrelated random walk of length N, fluctuation F(N), introduced in the DFA Wiki article, grows like a square root of N, i.e., its Hurst exponent is 1/2. Recall that the "fluctuation" of the uncorrelated time sequence such as a white noise realization does not depend on the sequence length by definition, its Hurst exponent is 0, and a white noise realization has its self-affinity at negative maximum -- the stretches of steady increase/decrease must be compensated by the stretches of the opposite direction, because the mean is zero. The alphas between zero and 1/2 correspond to an "anti-correlated" behavior from the POV of self-similarity. Likewise, the alphas between 1/2 and 1 correspond to correlation, and the alphas greater than one bear witness of the presence of non-stationary, unbounded process. The alpha=1 process is the stationary process, correlated at the maximum amount of self-affinity, it is the "1/f noise", because its PSD(f) is proportional to the inverse of frequency. The variance of an uncorrelated random walk is proportional to the square root of the series length; for the time series characterized by the Hurst exponent H, the variance and higher order moments are proportional to the series length raised to power H. This statistics results in the scaling behavior
$$
x(at) \sim a^Hx(t)
$$
where x is a random variable of the time series generating distribution, a is the scaling coefficient, H is the Hurst coefficient.
With H=1, the scaling behavior is closest to the concept of pure self-similarity, x(at) ~ ax(t), reproducing the system behavior at every time range.
These explanations do not prove the self-similarity claims, as mathematics requires it. But it would greatly help you learn the signal processing theory and master developer skills, if you compute and simulate the problems and processes. You can create an array of data and initialize it with a WGN samples in Matlab or Python; transform this array into a random walk process realization; compute trend lines with least squares and plot log-log graphs. Then, you can generate 1/f noise samples (or whatever noise spectrum you want) and plot trend lines for these distributions. Changing time scales, you will see for youself the claimed self-similarity properties. And in the process, you will master the necessary math.


With your edit it becomes clear that you've modeled your problem incorrectly:

While the offending signal appears shortly with a frequency of 50 Hz, that is by no means the frequency content of the interference!

(also, your filter isn't well-designed, probably too short, judging from the impulse response it displays, to filter out 50 Hz)

You'll find that interference to look more or less like dirac impulses – and these have a very high bandwidth! So, in spectral domain (plot a PSD!) you'll notice that there's not only a peak at 50 Hz, but wideband interference.

Technically, you seem to have a comb of such, and that would actually lead to a line spectrum again, which you could "kill" with a multi-notch filter.
But then again, these aren't really diracs, and you'd need a very long filter to even come close to being able to notch these out.

I think a time-domain solution is much more promising than filtering:

Simply build a threshold that detects when a burst interferer starts and ends, and zero out all samples in between.

If you want to be really elegant, you could try to use for example linear prediction to fill these gaps instead of zeroing, but my guess is that you'd only win a tiny degree of signal quality.


Recovering a signal after nonuniform sampling

I agree with your result for $Z(j\omega)$. Apart from scaling, what you have is the original spectrum with positive and negative frequencies swapped and multiplied with factors $1-j$ and $1+j$, respectively. In order to restore the original signal, we need to add right-shifted and left-shifted versions of the spectrum, while getting rid of the complex factors. Final lowpass filtering removes the redundant components in the frequency range $[M,2M]$.
The necessary shifting and complex scaling is achieved by multiplication with
$$e^{jMt}(1+j)+e^{-jMt}(1-j)=2\cos(Mt)-2\sin(Mt)\tag{1}$$
Modulation with $(1)$, lowpass filtering with cut-off frequency $M$, and appropriate scaling restores the original signal.


Stationarity is a multi-fold concept in signal processing. It can denote a wide range of behavior, encompassing deterministic or stochastic aspects. Beyond that, the main question is: do you know if your signal is stationary, and how?

If you actually know how, it is probably wiser to use the generation process to build a custom, adapted model or transformation, and use it for the analysis.

Even in that case, I strongly advocate using different analysis methods in parallel, to help you detect artifacts, issues than you would not detect with a single model. For instance, let us remind that one usually observe only a few realizations of a "signal", and that acquisition issues, outliers, etc. may occur.

Finally, analyzing in first intention a signal  with time-frequency or time-scale transforms is a good idea, as it can help you detect the useful scales of interest, estimate parameters of stochastic events, etc.

The drawbacks are:

The difficulties in choosing the appropriate wavelet (real or complex), and the associated sampling (and the resulting speed)
The difficulties in interpreting the scalogram, as a knowledge of the underlying processes could be useful


Channel estimation using bit reversed order of FFT matrix

So, according to your text,
$$ H := F^bh,$$
i.e. what you call the bit-reversed order FFT of the channel impulse response.
Since "bit-reversed order" just implies you're permuting outputs of the FFT, $F^b$ is just a row permutation $\Pi_{BR}$ of the standard DFT matrix $D$.
Since the "proper" DFT actually leads to a diagonal $H$ (as a linear channel never "mixes" any other frequencies in), you already know exactly what to do to convert your $H$ to a diagonal matrix:
Reverse the reverse bit ordering:
$$\Pi^{-1}_{BR}H=\Pi^{-1}_{BR}\left(F^bh\right) = \Pi^{-1}_{BR}\left(\Pi_{BR}Dh\right) = Dh,$$
which, in best classical OFDM nature, diagonalizes the channel.
(Using the full rank of any permutation matrix, you can also show this is the only way to diagonalize the result of the bit-reverse ordered FFT; putting this into a ray of optimism: Those who set out to invent OFDM will in expectation invent OFDM!.)
Channel estimation relates to the way you want to correct the channel influence. Hence, once you diagonalized your system, it becomes "normal" OFDM, and all the normal OFDM methods applies.


Difference between sticky channel Poisson-repeat channel

The Poisson-repeat channel is a binary channel where each transmitted bit is repeated at the channel output anywhere from $0$ to $\infty$ times. The repetition distribution is Poisson.

A "sticky channel" is a repeat channel where deletions are forbidden; that is, every transmitted bit is repeated from $1$ to $\infty$ times at the channel output, according to some specified distribution.


What happens to the phase spectrum when I resample?

1- Is this resampling correct?

No. Resampling is almost never 100% correct but it usually requires a complicated trade off between different artifacts and properties: passband extension and flatness, stop-band attenuation and size, aliasing rejection, phase distortion, onset & transient preservation, causality, CPU and memory requirements etc.
Which trade off is the best, depends a lot on the requirements of your specific application.

2- Is there a solution to create any kind of corresponding phase information for higher frequencies?

Why would you need that? Your original signal cannot contain any energy at frequencies equal or above the original Nyquist frequency. When you re-sample you extend the frequency range, but the energy at these new frequencies (between 11 MHz and 20 MHz) is close to or exactly zero. That's the whole point of a "good" interpolation: make sure you don't add stuff that wasn't there before.
That means the phase at the new frequencies is either undefined or doesn't really matter since there is no relevant energy. If the magnitude is zero, it doesn't matter what the phase is.


How to Select Point Spread Function Empirically for Image Deconvolution?

One way, though not directly visual, is to observe the image statistics.
We have a pretty good idea about the statistics of Natural Images, more specifically, their Gradient Distribution (See Statistics of Natural Images and Models by Jinggang Huang, D. Mumford, What Makes a Good Model of Natural Images by Yair Weiss, William T. Freeman, The Statistic Distribution of Image Gradient - Mathematics StackExchnage and UCLA Stat232A-CS266A Statistical Modeling and Learning in Vision and Cognition - Chapter 2 - Empirical Observations: Image Space and Natural Image Statistics).

So what you can do is see how well the statistics of the result matches the models.
Actually many modern (Prior to the Deep Learning era) Deblurring Methods work according to this approach with very nice results.

Look at Blind Motion Deblurring Using Image Statistics by Anat Levin.


How does both shifting and scaling a function affect the Fourier Transform?

Your logic is correct. There are two ways to determine the answer. Yours is one of them:

$$\begin{align}v(t)=x(at)&\Longleftrightarrow \frac{1}{|a|}X(\omega/a)\\v(t-t_0/a)=x(at-t_0)&\Longleftrightarrow\frac{1}{|a|}X(\omega/a)e^{-j\omega t_0/a}\end{align}$$

The other option is

$$\begin{align}w(t)=x(t-t_0)&\Longleftrightarrow X(\omega)e^{-j\omega t_0}\\w(at)=x(at-t_0)&\Longleftrightarrow \frac{1}{|a|}X(\omega/a)e^{-j(\omega/a) t_0}\end{align}$$

Of course, both options lead to the same result.


Canny's original paper addressed the issue of ridge edges and roof edges in addition to step edges, but I'm unaware of any implementation of "Canny" edge detection that detects anything other than step edges.  Canny's paper is a pretty easy read, but, unfortunately, paywalled.  J. Canny, "A computational approach to edge detection", IEEE Trans. Pattern Analysis and Machine Intelligence, 8:679-714, 1986.

Canny's paper is really about how to use calculus of variations to derive a 1-dimensional filter that has optimal response to each type of edge.  It turns out that the ridge-edge filter is not scale invariant (you have to know how wide your ridges are).  Here is the picture from Canny's paper of the optimal ridge and roof detection operators (page 683, © IEEE, 1986, this is "fair use"):

What is typically called "Canny edge detection" is as follows:

Perform an approximation of a derivative of Gaussian function in X and Y.  (This is usually done with a Gaussian blur followed by the Sobel operator.)  (The derivative of Guassian is actually not the optimal step-edge detection operator, but the optimal step-edge detection operator "kinda-looks-like" a derivative of Gaussian.)
Convert from $\langle\Delta x, \Delta y\rangle$ form at each pixel to $\langle$magnitude,"angle"$\rangle$.  (Where "angle" is a 2-bit approximation that tells you whether the edge is closer to horizontal, vertical or one of the two 45-degree diagonals.)
"Thinning": For each pixel look at its two neighbors in the "angle" direction and keep the pixel only if it is the maximum of itself and its two neighbors.  (Otherwise set the magnitude to 0.)
"Thresholding": For each pixel with magnitude less than some "weak threshold" suppress the pixel to 0.  For each pixel with magnitude greater than some "strong threshold" set the pixel to whatever value means "definitely an edge."  For pixels between the two thresholds set the pixel to the value "maybe an edge".
Hysteresis: For each "maybe an edge" pixel change it to "definitely an edge" if any of its 8 neighbors in the horizontal, vertical or diagonal directions is "definitely an edge."

Theoretically, you should be able to turn a Canny edge detector into a ridge-edge detector by changing just step 1, but I've never tried it.


Practical implementation of Expected Value?

The answer is yes. Even if we could generate several realizations of the random process by conducting the random experiment a number of times. We could only have a limited number of realizations and for a limited period of time practically. Now to take expectation of the random process we will have to take freeze time and look at the distribution of the values across the several realizations we generated. But that would be an estimated expectation because we can only estimate the distribution, because although we have conducted the experiment several times, we cannot practically be sure to exhaust the sample space.


A basic question about the use of moving average vs low-pass filters in DSP

A moving average filter can be thought of as a type of low-pass filter that doesn't have any control over its bandwidth for a fixed number of taps.

For a finite impulse response (FIR) filter, the output signal $y[n]$ is given in terms of the input signal $x[n]$ and the filter taps $h[n]$:
$$
y[k] = \sum_{n=0}^{N-1}h[n]x[k-n].
$$
The filter length in this case is $N$ "taps" (or $N$ coefficients). A moving average filter has coefficients that are all equal:
$$
h[n] = \frac{1}{N}, \qquad n = 0, 1, \ldots, N-1,
$$
whereas in general, a low-pass filter (LPF), can have different values for each tap. This allows you to control the frequency selectivity of the filter.

Here is an example showing the difference, for a system with a sample rate of 200 Hz and a filter with 32 taps:

The LPF was designed to have a cutoff frequency of 50 Hz, but it can be set to any value you want. With a moving average filter the filter is narrowly focused around the 0 Hz component ("DC"), and the peak gets narrower the more taps you have in the filter. Another problem with using a moving average filter as an LPF is that it has high sidelobes (the ripples to either side of the main peak) compared to a "properly designed" filter.

The taps for the LPF are:

[-0.00116417 -0.00139094  0.00195951  0.00293134 -0.00437535 -0.00637313
  0.00902803  0.01248116 -0.0169409  -0.02273977  0.03045372  0.04118039
 -0.05729369 -0.08500841  0.14720004  0.45005217  0.45005217  0.14720004
 -0.08500841 -0.05729369  0.04118039  0.03045372 -0.02273977 -0.0169409
  0.01248116  0.00902803 -0.00637313 -0.00437535  0.00293134  0.00195951
 -0.00139094 -0.00116417]

Python code for generating this figure:

from pylab import *
from scipy.signal import windows, firwin

Fs = 200.

N = 32
h_lpf = firwin(N, 50, nyq=Fs/2., window='hamming')
h_ma = ones(N)*1./N

print h_lpf

M = 512

X_lpf = fftshift(abs(fft(h_lpf, M)))
X_lpf /= X_lpf.max()
X_ma = fftshift(abs(fft(h_ma, M)))
X_ma /= X_ma.max()

f = arange(M)*Fs/M - Fs/2.

figure()
plot(f, X_lpf, f, X_ma)
xlabel('$f$ (Hz)')
ylabel('$|H(f)|$')
legend(('LPF', 'Moving Average'))
grid()
show()


Can some one give me an example of cross correlation's high value where template matches with image?

Please read this document here. If I understood your question clearly, then you need to take out the mean from the image. After that probably you would find cross-correlation for the matching patch.


Accuracy of time domain signal when one frequency is absent

When a signal is expressed as a linear superposition of orthogonal basis functions, like a DFT, when you remove one orthogonal component, that error is orthogonal to the remaining terms, so any perturbations to those remaining terms will not have an effect on the error due to the removal.

Nevertheless, missing "data" is an active area of research.  If you can make some assumptions about what is missing, you can sometimes, with varying success recover that data.

There is a very popular family of algorithms known as the "EM algorithm"

Another more recent area is known as "Matrix Completion"

There is a mathematical property called "restricted isometry" that is used in sparse sensing that may of use in your problem.

When you google scholar the term "incomplete spectrum", you get hits like:

  E. Serpedin, "Subsequence based recovery of missing samples in
  oversampled bandlimited signals," in IEEE Transactions on Signal
  Processing, vol. 48, no. 2, pp. 580-583, Feb 2000. doi:
  10.1109/78.823989 Abstract: A new approach for recovering an arbitrary finite number of missing samples in an oversampled bandlimited signal
  is presented. This correspondence also proposes an approach for
  recovering the original signal's spectrum from the spectra of a
  certain number of subsequences, obtained by downsampling the original
  sequence. Closed-form expressions for the missing samples in terms of
  the known samples are obtained by exploiting the linear dependence
  relationship among the spectra of downsampled subsequences keywords:
  {bandlimited signals;signal reconstruction;signal sampling;spectral
  analysis;closed-form expressions;downsampled
  subsequences;downsampling;known samples;linear dependence;missing
  samples;oversampled bandlimited signals;spectrum;subsequence based
  recovery;Bandwidth;Closed-form solution;Equations;Fourier
  transforms;Frequency domain analysis;Iterative algorithms;Iterative
  methods;Numerical stability;Sampling methods;Signal restoration}, URL:
  http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=823989&isnumber=17842

If you can express your problem in terms of any of these paradigms, you can arrive at solutions, and some have bounds on errors.

As was suggested in the question, an often useful simple heuristic is to interpolate the missing data. The most straightforward way to get an handle on error, is trial and error on actual data.  A lot of times, acts of technical desperation work out.


Following, a way to remove the 3 exponential pulses.

1.- Rebuilding signal + interferer from graph provided

clear all;clc;close all

A=imread('001.jpg');
A1=A(:,:,1);
figure(1);imshow(A1)

A2=~imbinarize(A1);
figure(2);imshow(A2)
% 001

problem: bwmorph breaks causality and splits signal

A3=bwmorph(A2,'skel',Inf);
figure(3);imshow(A3);
% 002

[sz1 sz2]=size(A2)

for k=1:1:sz2
    L1=A2(:,k);
    n1=find(L1);

    A2(:,k)=zeros(1,sz1);
    if ~isempty(n1)
        A2(floor(mean(n1)),k)=1;
    end
end
figure(4);imshow(A2)
% 003

Now signal is correctly mapped.

Since the contaminated signal is available to the question originator, the original data should be inserted here:

s=0
for k=1:1:sz2
    L1=A2(:,k)
    s=[s find(L1)];
end
s(1)=[]

2.- where is mean(s)?

don't do just s=s-mean(s)

Left behind the y=0 axis tick on purpose

% 004

mean_s=287

s=-(s-mean_s)

3.- Rebuilding time reference

t=linspace(0,10,sz2)
dt=mean(diff(t))

4.- Plotting signal + interference

figure(5);plot(t,s);grid on
xlabel('t[seconds]');ylabel('s+i')  % signal+interferer

Comments:

The VLF bending within the initial 5 seconds may come from conducted interference; contact, applying an ECG sensor on a mobile phone?

The base station or mobile seeking answer, starts with certain power and
then as the packets are exchanged Base Station Mobile Terminal tell one another to increase
or reduce power, usually the Base Station starts higher than needed
to make sure reaching cell edge, and then reduces power until MS says
level ok.

The negative sharp spikes show there may be something metallic reflecting the interfering signal.

You have provided a reference signal, what it should be like, but there's
no time reference, just the amount of samples, without further
information, it's not sure the reference signal duration is also 10 seconds, is it?

[pks,locs]=findpeaks(s,'MinPeakHeight',50,'MinPeakDistance',50)
s_peaks=s(locs)

5.- When do the 3 exponential decay interfering pulses start?

t_peaks=t(locs)

6.- What is the mean delay between those 3 pulses measured in time?

diff(t_peaks)

7.- What is the mean delay between those 3 pulses measured in amount samples?

diff(locs)

8.- What is the real amplitude of the pulses? approximately the 1st sample

ds_peaks=s(locs)-s(locs-10)

9.- what is the amount of samples per interfering pulse?

nT1=floor(mean(diff(locs)))

10.- What are the respective attenuation constant of the 3 exponential pulses?

a=-1./(dt*nT1)*log(1./ds_peaks)

11.- 1st interfering pulse start time measured directly on the graph

t0=.688

12.- What is the sample numeral of the start moment of the 1st interfering pulse?

nt0=[1:1:floor(mean(intersect(find(t>t0-.01),find(t<t0+.01))))]

13.- Values of each interfering pulse:

e1=ds_peaks(1)*exp(-a(1)*dt*([1:1:nT1]+t0))
e2=ds_peaks(2)*exp(-a(2)*dt*([1:1:nT1]+t0))
e3=ds_peaks(3)*exp(-a(3)*dt*([1:1:nT1]+t0))
e_interf=[zeros(1,numel(nt0)) e1 e2 e3 zeros(1,sz2-3*nT1-numel(nt0))]
hold all
plot(t,e_interf)

s2=s-e_interf
figure(6);plot(t,s2);grid on
xlabel('t[seconds]');ylabel('s+i-3pulses')  % signal+interferrer

If interested I can show how to use findpeaks to remove the sharp negative pulses.


What Is the Difference between Parseval's Theorem and Plancherel Theorem?

According to the links you shared, and using Wikipedia's notation, Parseval assumes $A(x),B(x)$ are two square integrable (with respect to the Lebesgue measure), complex-valued functions on $\mathbb {R}$  of period $ 2\pi $  with Fourier series

$$A(x)=\sum_{n=-\infty}^\infty a_ne^{inx}$$
and
$$B(x)=\sum_{n=-\infty}^\infty b_ne^{inx}$$
Then
$$\sum_{n=-\infty}^\infty a_n\overline{b_n} = \frac{1}{2\pi} \int_{-\pi}^\pi A(x)\overline{B(x)} \, \mathrm{d}x$$
Plancherel provides a continuous version of Parsevals theorem (in a sense that $a_n$ and $b_n$ are now the $f(x)$ and $g(x)$, respectively and $A(x),B(x)$ are replaced by $ \widehat{f}(\xi), \widehat{g}(\xi)$)


Weird note structure on guitar audio waveform

Here are some possibilities:

Transverse waves on a guitar string have two polarized components, one parallel and one perpendicular to the guitar's soundboard. User hotpaw2 mentioned energy exchange. It is imaginable that energy transfer could happen between the components. The component parallel to the soundboard is largely inaudible, so it could stealthily store energy until released later.

Or it could be that you have two strings with modes that are slightly out of tune, and they form a beat frequency.

EDIT: Having analyzed the audio sample, the things you have marked are just some low-frequency artifacts, I think from things like someone touching or shaking the mic or moving around and bumping the guitar body. Adding a 100 Hz high-pass filter mostly gets rid of those:

Figure 1. Sample after 100 Hz highpass filtering.


Understanding the constraint on distance between microphones in microphone array for Direction of Arrival estimation

Say the microphones $A$ and $B$ are omnidirectional and exactly half a wavelength apart. If the sound field propagates in the direction of the vector $A\to B$, then a chosen part of the waveform will arrive at $...,$ $A,$ $B,$ $A,$ $B,$ $...$ at intervals that equal half the period of the waveform. If the sound field propagates in the direction of the vector $B\to A$, the chosen part of the waveform will arrive at $...,$ $B,$ $A,$ $B,$ $A,$ $...$ at identical intervals. You won't be able to tell the difference between the two cases from the arrival times because for a periodic waveform there is no "first arrival". This ambiguity is resolved by the constraint you mention, as you then know that you can safely choose for comparison arrival times at $A$ and $B$ that are separated by at most some known time less than half the period of the waveform.


Splitting a magnitude across multiple bins in an IFFT

Your frequency in term of bin number is approximately:

SignalF / binDelta = 110 / 86.1 = 1.28

So, for the positive frequency image :

w0 = pi * (0 - 1.28)
w1 = pi * (1 - 1.28)
w2 = pi * (2 - 1.28) = pi * 0.72
wx = pi * (x - 1.28)
binmag(x) = sin(wx) / wx

Also, one needs more than 3 points to avoid windowing artifacts.

But the big problem with Sinc interpolation for any f near bin 0 (DC) is that the signal plus its complex conjugate mirror image at the nearby negative frequency -f both wrap across bin 0, and thus constructively or destructively interfere within all those nearby bins, depending on the phase of f measured relative to the center of the window.  See http://www.nicholson.com/rhn/dsp.html#4 for a more complete Sinc (or Dirichlet) interpolation formula.  A similar problem occurs for frequencies near bin N/2 for an IFFT of length N. This creates a problem with IFFT (re)synthesis of any fractional frequencies near DC or N/2.


We want to find a representation in terms of sines and cosines for the function
$g(x)$ which has value $x$ for $x \in [0,\pi]$. We don't really care if this
representation does not work well outside the interval $[0,\pi]$, but
inside the interval, the representation better work perfectly.  This seems
to be crying out for a Fourier series approach, and so that's what we will
use.  The standard approach begins with functions that are periodic with
period $2\pi$ and so let us extend the definition of $g(x)$ to cover
one period $[0,2\pi]$.

If we choose $g(x)$ to have value $2\pi-x$ for
$x \in [\pi, 2\pi]$, then, when repeated periodically,
this extended $g(x)$ is an even function of $x$ with average value $\pi/2$.
Furthermore, as is the case with all even functions, the sine terms all
have $0$ coefficients.  Thus, the Fourier series works out to
have only cosines in it, and the DC term is $\pi/2$.
If we choose $g(x)$ to have value $x$ for $x \in [\pi, 2\pi]$, then
when repeated periodically,
this extended $g(x)$ is not an odd function of $x$, but nonetheless
it is true that all the cosine terms in the Fourier series
turn out to be $0$ except the DC term which has value $\pi$.


The Goertzel algorithm allows you to compute individual terms of the Discrete Fourier Transform, and is more efficient than the FFT. However, if you wish to later compute the spectrum of other bands besides alpha, Goertzel will turn out to be more computationally expensive than the FFT.

Unless you are seriously limited by computational power, in principle, you can compute the Fast Fourier Transform in quasi-real-time, using a sliding window (also known as the Short Time Fourier Transform, or STFT). The Discrete Wavelet Transform is also sometimes used for real-time EEG power estimation. To find the power in the alpha band, you can look at the corresponding frequency bins of the STFT for either the original or the filtered signal.

The output of the STFT will be delayed by the length of your window (as you need all the samples inside the window to arrive before computing the power in that time bin). However, with a window length of 512 samples under 8kHz sampling rate, for example, this amounts to 0.064 seconds delay, which is fairly close to 0.0 seconds!

The STFT (Short-Time Fourier Transform) is effectively a computation of the FFT one 'chunk' or 'time-frame' of the signal at a time. Often these chunks are overlapped by some amount, such as 50%, meaning that if the window length is 512 samples, then halfway through the first frame (at the 256th sample), you already begin to process samples that contribute to the second time-frame.

The purpose of using overlapping frames is to avoid artefacts caused by taking square slices of the signal. (Computing the power of one sample at a time, like you suggest, amounts to convolving your signal with a square window, which is equivalent to convolving the frequency representation of your signal with a sinc function. This causes spectral leakage, or 'ringing' artfeacts. Thus, it is common to use a sliding window with smooth edges, such as a Hamming window, to reduce these artefacts.)

Note that the choice of window length will affect your spectro-temporal resolution. A longer window will have better frequency resolution, but worse temporal resolution, and vice-versa. This phenomenon is known as the Gabor limit, or Heisenberg-Gabor limit.