A PassBand is the Band that can be passed through a BandPass filter, or for that matter, any filter. In other words, the filter is a Bandpass, what signal passes through is a Passband.
The classic is a high pass and low pass together. For instance, a capacitor, inductor, and resistor in series is an example of a bandpass filter. Digitally, you combine a high pass and a low pass filter together.
Essentially, if you are limited to a small region, you can infer what frequency you are limited to. The classic example is Nf->2*Nf. The frequencies reported at 0 is 2*Nf, the frequency reported at Nf is actually NF, and continue the pattern for the rest. This only works when you are guaranteed to only have frequencies in the particular band.


Yes, that conclusion would be valid. However, as you suspected, there are some limitations. Consider the case where the system's frequency response $Y[k]$ contains one or more zeros. In that case, corresponding frequency bin in your estimate of the input signal $X[k]$ would diverge to infinity (because of the division by zero).

What you're really trying to perform by doing this sort of calculation is deconvolution, which is a complex topic in itself, but suffice it to say that deconvolution is often difficult; you illustrated one of the reasons why. This naive approach of just dividing in the frequency domain is known as a zero-forcing type of algorithm.


Are Convolution and Deconvolution Kernels the Same?

Question - When Are the Convolution Operator (Kernel) and the Deconvolution Operator (Kernel) the Same?

Discrete Convolution (Cyclic) is described by Circulant Matrix.
Circulant Matrices are diagonalized by the DFT Matrix (Will be denoted as $ F $).
So given a convolution kernel $ g $ with its matrix form given by $ G $ one could have its diagonalization by:

$$ G = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F $$

Where $ \mathcal{F} \left( g \right) $ is the Discrete Fourier Transform (DFT) of the kernel $ g $.

In the general case the deconvolution is given by:

$$ x = G h \Rightarrow {G}^{-1} x = {G}^{-1} G h = h $$

Where $ {G}^{-1} $ is the Deconvolution Operator (Inverse of $ G $ the convolution operator).
In our case we're looking for the case $ G = {G}^{-1} \Rightarrow G G = I $.

Utilizing the diagonalization of Discrete (Cyclic) Convolution Operators by the DFT Matrix (Which is Unitary Matrix):

$$ G G = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F $$

Namely we need $ \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) = I $ which means $ \forall i, \, {\mathcal{F} \left( g \right)}_{i} \in \left\{ -1, 1 \right\} $, namely the Fourier Transform of $ g $ has the values of -1 or 1.
So there is a set of kernels $ \mathcal{S} $ which is defined by:

$$ \mathcal{S} = \left\{ g \mid \forall i, \, {\mathcal{F} \left( g \right)}_{i} \in \left\{ -1, 1 \right\} \right\} $$

One should notice that the identity Kernel ($ I $) is indeed part of this set of kernels.

The above kernels obey $ \delta \left[ n \right] = \left( g \circ g \right) \left[ n \right] $ where $ \circ $ is the Circular Convolution (Which is multiplication of the DFT).

How could one generate such a kernel?
Well, just generate a vector of numbers which each is composed of $ 1 $ or $ -1 $ and use its ifft().
Pay attention that if you are after real operator the vector must obey the symmetry of the DFT for Real Numbers.

You may find a code to replicate the above in my StackExchange Signal Processing Q19646 GitHub Repository.


If your goal is to plot a magnitude spectrum with the frequency axis on a log scale, but with a roughly even visual resolution along that axis, then a single FFT might provide too low a density of plot points (without interpolation) at low frequencies, and more plot points than can be plotted on a line (without averaging) at high frequencies for a given print or pixel display resolution or smoothness.  If you use different FFT lengths for different octaves or sub-octaves in frequency, and select subsets of the results from each FFT, then you can maintain a lower delta in density of FFT result points when plotted on such a log scale.  How many FFTs you might want to use depends on the maximum variance you want in log frequency resolution of the final joined plot.

However, since the FFTs are of different lengths, then they are for different sets of data.  A sequence of longer FFT windows can be done with larger offsets if more time resolution isn't needed at low frequencies.  You will have the problem of how to join or blend all these different FFT results (usually subset segments of the results) so that the magnitude response might be uniform across FFT boundaries.  In total, you also end up with an overdetermined set of FFT results (but so will highly overlapped windows).

Even better results (in terms of even log frequency plot resolution) might be obtained by using some form of wavelet transform (a Morlet or Gabor wavelet or constant-Q transform, for example) instead of a bunch of semi-redundant FFTs.

But libraries for optimized FFTs might be more available on some platforms.  Each basis vector of a windowed FFT in a "multi-time" set can be considered a non-optimally-sized wavelet.  So, in some ways, using these so-called "multi-time window" FFTs is a "poor man's" wavelet transform.  But some of the many FFT results might be useful for other forms of more traditional DFT analysis of the signal input (in conjunction with the log plot), thus serving a dual purpose.

A partial psycho acoustic justification for increasing FFT size at lower frequencies (or using wavelets) is that, over a certain mid-frequency range, the human ear/brain takes a length of time to determine frequency or pitch of a sound roughly proportional to the period of the frequency.


Well, we are summing over $k$ in the range $n_0$ to $n$. The number of elements in the sum is therefore $n-n_0$. We also know that the sum must be positive and that all terms in the sum ($x[n]$) have absolute value less than $M$.

To account for $n-n_0$ being negative, but our sum MUST be positive, we take its absolute value.


Extended Kalman Filter in mechanics, electronics and hydraulics?

It's the application which dictates whether you should use a linearized, extended, unscented Kalman filter or you can go easy with the simpler linear Kalman filter.

In other words, the extended Kalman filter is not preferred over the linear one because it provides better performance than the latter. Just the opposite; its estimation error is larger than what could have been provided by an equivalent linear Kalman filter if the problem could be solved by the linear.

When the physical nature of the problem does not permit a linear mathematical model, then you should either abondon the Kalman filter at all, or you should find a way around (such as the linearization of the nonliner model) to continue using the Kalman filter and still expect its benefits unless you have some alternative algorithm to go with.

Most fundamentally the accuracy of the extended Kalman filter depends on the the sampling interval (step size) used to update the linearized matrices and other parameters. Hence you should use dense enough sampling so that the nonlinearity does not change too much to add excess linear apprximation error. However smaller sampling interval also means more computations per second.

Also subtle issues such as noise de-coupling, independence assumptions, dimension reduction have thier own solution strategies which can also affect the accuracy of the filter estimates per given sampling rate.

When your problems are inherently nonlinear, then you cannot use the linear Kalman filter, therefore you can either use the extended Kalman filter or think about other algorithms that could be useful.


The first part of the response may be a result of the transient response of the Butterworth filter.  Also note that the signal is phase shifted and the basic shape is altered.  This looks like the result of phase distortion.  You may want to try again using a Bessel filter which has better transient and phase characteristics.  If you are using DSP, you can also try a true linear phase FIR approach.  This will yield the best transient and phase characteristics.


Autoregressive Exogenous model on multiple datasets in MATLAB

It sounds like you need a model on top of your models, which sounds similar to a Hierarchical Bayesian Model:

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

I used this technique on a set of linear regressions, essentially following

http://seor.vse.gmu.edu/~klaskey/SYST664/Bayes_Unit7.pdf.

This is not AR modeling but perhaps might be a fruitful Google set of search words.

You might want to crosspost to Cross Validated

https://stats.stackexchange.com/


For a signal that rides at VSS/2, goes down to 0V and up to VSS, just use a pair of GPIO pins, each with a $10\mathrm{k\Omega}$ resistor, with the far ends connected together.  When both pins are high you'll get VSS, when both are low you'll get 0V, when one is high and the other low you'll get VSS/2.
Then sequence the pins to get the pattern you want.  If you need something that's faster and has less jitter than you can accomplish by bit-banging the pins, use a pair of timer outputs and sequence them in hardware (which you may not be able to do with the ATMega found in most Arduinos, and which the Arduino hand-holding software may not allow you to do easily -- it's your job to figure that part out).
To get your 2V, 3V, 4V pattern, you'll need up to four resistors, connected together in the middle, and to VSS, GPIO_n, GPIO_{n+1}, and ground on the other ends.  Finding their value takes some algebra, but it's straightforward if you know basic circuits.


From pressure to SPL : microphone calibration using sound level meter calibrator and sound level meter

First, I would check the positioning of sound level meter and microphone. If you are putting them close to the oscilloscope a few centimeters of displacement may change the measured value dramatically. (6dB/doubling of distance)

Second, which microphone do you use? The link you provided leads to a preamplifier. The actual capsule can either be a free field or a pressure field one. A free field capsule should be facing the sound source directly, while a pressure field capsule should face it sideways.


In MATLAB you can create a .wav file waveFileName.wav from a waveform y with sampling frequency Fs using MATLAB's audiowrite as follows:
filename = 'waveFileName.wav';
audiowrite(filename,y,Fs);

You can also read and output a waveform y and its sampling frequency Fs from the .wav file waveFileName.wav using MATLAB's audioread as follows.
[y,Fs] = audioread('waveFileName.wav');

I am sure there must be corresponding functions in Python.


MUltiple SIgnal Classification (MUSIC) Implementation

I'm going to group number 1 and 3 as related.

For a high level description of MUSIC, you can take a look at MATLAB's overview here. One of the main steps in the algorithm is to find the eigenvectors of a correlation matrix, which can be done via singular value decomposition or other methods. MATLAB has functions for this, so you may want to find equivalent functions in C/C++, or implement your own. The other steps in the algorithm are more implementation agnostic: mostly data collection and simple vector operations.

Using an N-element array, MUSIC can resolve at most N - 1 sources. These can be real sources or their images and must be dealt with accordingly.

For some more specific implementation details, I have an answer here where we go through the process of collecting data and performing MUSIC in MATLAB.


So the issue is that your filter order is too high.  There are 2 problems with this:

SciPy has a bug that generates inaccurate filters at high orders.
On any platform, higher-order filters cannot be done in a single stage.

SciPy bug:
SciPy was previously generating prototype filters as a list of poles and zeros, then converting them to transfer functions, then transforming them to their new frequency and band-type, then converting them back into poles and zeros.  This is bad because transfer function representation has numerical problems and loses precision.  So your filter currently has these poles and zeros:

but it's supposed to look like this:

This probably isn't your problem, since these poles and zeros are still approximately the bandpass filter you asked for.  At even higher orders the poles can be so wrong they jump outside the unit circle and make the filter explode.  This has been fixed in SciPy 0.14.0, as long as you use output='zpk'.
Second-order sections:
When you actually do the filtering, you have to use a transfer function, not a list of poles and zeros, so your high-order filter would suffer the same numerical corruption if you tried to do the entire filter in a single stage.  What we do instead is (starting from the list of poles and zeros, not from the inaccurate transfer function) break up the filter into second-order sections, which each have minimal numerical problems, and cascade the sections together to get the final equivalent output.
This is true for any platform, hardware or software, but Matlab and Octave have functions like sosfilt to do the work for you.  SciPy does not, but it will in the future.
For now, you can break up the sections manually and use lfilter repeatedly with each section to get the desired output or copy the code from that pull request.
Update
SciPy has SOS now.  So the best general-purpose filtering method is to generate your filter in SOS format, something like:
sos = butter(order, [low, high], btype='band', output='sos')

and then filter using sosfiltfilt:
y = sosfiltfilt(sos, data)


Having a look into the sources of librosa.display.specshow reveals how bins are converted into frequencies internally: The plotting uses librosa.core.fft_frequencies, which shows that it is basically the same as numpy.fft.fftfreq:
>>> librosa.fft_frequencies(sr=22050, n_fft=16)
array([     0.   ,   1378.125,   2756.25 ,   4134.375,
         5512.5  ,   6890.625,   8268.75 ,   9646.875,  11025.   ])


Can channel capcity be explained in simple terms?

There is a limit where simplifying a scientific concept starts becoming inaccurate.

However, the capacity of a channel is a really simple idea that is completely illustrated by the following equation:

$C = B\log_2(1+\frac{S}{N})$

The equation says that the maximum amount of information (in bits per second) that can be transmitted in a channel error-free is proportional to the Bandwidth $B$ in Hertz, and is also proportional to the logarithm of base-2 of 1+ the signal to noise ratio.

So let us break this down to two points, the bandwidth and signal to noise ratio.

The bandwidth is the range of frequencies that you use in your channel.  Let us take an example where your channel has a really small bandwidth ($B=0$) where you only have one frequency that you can operate on. This translates to a single sinusoid of a certain frequency. Note that you cannot turn it on or off, or even change its phase because that will result in a range of other frequencies. Therefore communication is impossible on this channel.

The signal to noise ratio, is a measure of how strong a signal is to how strong the noise is and can take values between 0 and $\infty$. A signal power equal to 0 or a noise power equal to $\inf$ will make $\log_2(1+S/N)$ go to zero where communication is impossible.

On the other hand, 0 noise power, or an infinite signal power (given the bandwidth is not 0) will mean that you can send as much data as you like through the channel without any errors. Which is also the case for an infinite bandwidth and a non-zero value of $\log_2(1+S/N)$.


You can easily create measurement for Kalman Filter model.
Just run the model itself on some noise with starting point and you'll have endless number of data points.

Namely, if your model is given by:

$$ x \left[ k \right] = F x \left[ k - 1 \right] + v \left[ k \right] $$

Then all you need is to generate $ x \left[ 0 \right] $ randomly according to $ P \left[ 0 \right] $ and start moving it according to the above model while generating $ v \left[ k \right] $.


Discrete filters are divided in two major classes: Finite Impulse Response (FIR) filters and Infinite Impulse Response Filters (IIR).

Both types of filters are structured around a fundamental unit. A fundamental operation that is responsible for the filtering. That is, the delay unit. Of course, besides the delay unit you also have the sumation points. But sumation, on its own, without delay doesn't produce filtering.

So, both filters, produce outputs that depend on the past of the inputs (and / or outputs).

The concept of "filter order" is defined for both types of filters but in a slightly different way (it's still the same thing really).

Let's look at an FIR filter:

$$y_n = \sum_{m=0}^{Nh-1} x_{n-m} \cdot h_m$$

Here, to produce one value for $y_n$ we sum $m$ products of the past of the signal and some impulse response $h$. $Nh$ is the total length of $h$.

The order of the filter here is $N_h$, directly proportional to the number of delay units in the filter.

So, by increasing the order, we do get a sharper filter but also a longer filter and the longer the filter is, the longer it takes for one of the $x_n$ to propagate through the series of delays.

And that is where the increased group delay is coming from! Furthermore, in an FIR filter all frequencies are delayed proportionally. So, bass frequencies are delayed less than treble frequencies in a predictable (and linear) way because the time it takes for the sinusoids to propagate through the delay line is fixed.

Now, let's look at an IIR filter:

$$y_n = \sum_{m=0}^{Nh-1} x_{n-m} \cdot h_m - \sum_{k=1}^{Ng-1}y_{n-k} \cdot g_k$$

Which, if you take a closer look, is like two FIR filters back-to-back. The left part, we have seen before and it is the FIR filter from the first equation. The right part is an FIR but this time it is on past values of the output. We therefore have introduced feedback...and all the "joys" that come with it.

The order of an IIR filter is now the length of the feedback delay line (i.e. the $Ng$ here) and because of this interaction between the input and output, the group delay of an IIR filter varies with frequency. So, treble frequencies might be going through much faster than bass frequencies in an IIR filter.

OK, so now let us look at your experiment. You design a filter and you pass through it a simple sinusoid at some fixed frequency. Then you increase the order and you observe a longer delay.

Which is perfectly reasonable, because, irrespectively of the type of filter (FIR / IIR), increasing the order, means adding more delay units.

Hope this helps.


ZP-OFDM vs. CP-OFDM channel equalization

Hint: you could convert zero padded ofdm into cyclic ofdm at the receiver, and use the same equalization. Though there will be more noise on this scheme.
I had asked a puzzle some time back on this topic:
Modelling zero padded OFDM as circular convolution between channel and FFT frame


The a priori SNR is the ratio of the power of the clean signal and of the noise power. The a posteriori SNR is the ratio of the squared magnitude of the observed noisy signal and the noise power. Both SNRs are computed for each frequency bin.

Of course, the only signal we have is the observed noisy signal. The noise power as well as the power of the clean signal must be estimated. The noise power can be estimated during absence of the desired signal, but in practice it can be difficult to detect the presence or absence of the desired signal, at least when the SNR is bad. There are so-called noise tracking algorithms that estimate the noise power without relying on detecting the presence or absence of the desired signal. They are based on the assumption that the desired signal is "more non-stationary" than the noise. The power of the desired clean signal can be estimated using a decision-directed approach, as described in Ephraim and Malah's paper. In that approach, the estimated clean amplitude, i.e. the output or 'decision' of the algorithm, and the a posteriori SNR are used to estimate the a priori SNR.


The frequency spectrum of a static non-linearity driven by colored noise

If your nonlinearity can be expressed as a polynomial (i.e., in terms of addition and multiplication), you can make use of:

The linearity of the Fourier transform, i.e., if $f$ and $g$ are (benign) functions, $a$ and $b$ are numbers and $ℱ$ denotes the Fourier transform, then:

$$ℱ(a·f+b·g) = a·ℱ(f) + b·ℱ(g)$$
The convolution theorem, which states that multiplication in the time domain corresponds to convolution in the frequency domain (and vice versa):

$$ℱ(f·g) = ℱ(f) ∗ ℱ(g)$$

So, in your specific example, the spectrum of $y$ is the spectrum of $v$ convoluted with itself:

$$ℱ(y) = ℱ(v·v) = ℱ(v) ∗ ℱ(v)$$


Expression for Frequency-shifted Gaussian pulse

Well the expression for the $\frac {d^n}{dt^n} p(t) $ is $V_n(t) p(t)$, where $V_n(t)$ is a polynomial of degree $n$, and $V_n(t)$ may be determined recursivelyy
$$  V_{n+1}(t)e^{-t^2/(2\tau^2)} = \frac{d}{dt}\left(V_{n}(t)e^{-t^2/(2\tau^2)}\right) $$
Solving for $V_{n+2}(t)$ we have
$$ V_{n+1} = \frac{d}{dt} V_{n}(t) - \frac{t}{2\tau^2} V_n(t)$$
In order to link the two forms of the equation we need
$$\begin{eqnarray}
 g(t) &=&  p(t)cos(2\pi f_{c}t) \\
      &=& V_{TX}exp(-\frac{t^{2}}{2\tau^2})cos(2\pi f_{c}t) \\
      &=& \frac{d^{n}}{dt^{n}} \left( V_{TX }exp(-\frac{t^{2}}{2\tau^2}) \right) \\
      &=& p(t)V_n(t)
\end{eqnarray}$$
Thus $cos(2\pi f_c t) = V_n(t)$, the right hand side is a polynomial. As $n \to \infty$, the are the Hermite polynomials have some asymptotic expansions in terms of cossines. So the two may assume approximate values for certain values of $f_c$ and $\tau$.
Numerically evaluating the derivaties we see some resemblance with gausiand windowed cosine.
function hermite;
  tau=pi;
  t = tau * linspace(-5, 5, 1000);
  for order = 4:2:20
    subplot(3,3,order/2-1)
    a = -1/(2*tau^2);
    P = V_coefs(order, a);
    plot(t, polyval(P(end:-1:1), t).*exp(a*t.^2));
    title(num2str(order, 'n=%d'))
  end

end

function coefs = V_coefs(order, a)
  coefs = 1;
  for ii = 1:order
    coefs = [(1:(ii-1)).*coefs(2:end), 0, 0] + [0, coefs .* (2 * a)];
  end
end


What you understand is right. The second peak is the negative frequency component. I guess what your professor meant is that the Fourier transform of a sine wave is theoretically perfect Dirac $\delta$ function, and the discrete Fourier transform of a sine wave is a Kronecker $\delta$ function, but in your implementation, it's two peaks instead of two $\delta$ impulse.
The problem is the relationship of the frequency of sine wave and FFT size. As we know the frequency resolution of FFT is
$$
\Delta f = \frac{f_s}{N}
$$
where $f_s$ is sampling rate and $N$ is the length of FFT.
The FFT results in perfect $\delta$ only when the sinusoid frequency is exactly integer multiples of $\Delta f$.
In your case $\Delta f=15.625$ Hz, try something like 125 Hz and see what happens.


How many psychoacoustical critical bands are there?

Neither. The basilar membrane in your ear performs a frequency-to-place transformation which is then picked up by an array of thousands of hair cells (cilia). Therefore there are thousands of heavily overlapped bands. The processing that follows is very complex and not completely understood. The shape of the effective band pass filter at each location on the membrane can be measured in animals, and the masking curves can be measured by listening tests. Masking is measured by increasing the amplitude of a test tone relative to the masking tone until it is audible, so the situation is a bit different that what you assumed in your question (where you assumed that if the masker and the test tone were in the same band, you only hear one tone). The actual result is that as the test tone approaches the masker frequency, it must have a higher amplitude before it can be detected.

Bob


Draw a structure for IIR Filter given H(z)

It would help to apply some specificity to the general equation:

$H(z) = \frac{b_0 + b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}$

Rearrange into a difference equation, skipping the steps for brevity:

$y[n] = b_0x[n] +b_1x[n-1]+b_2x[n-1]-a_1y[n-1]-a_2y[n-2]$

Does the difference equation form show now that the Direct Form structure is derived purely by inspection?


What would the plot of an two sided digital signal look like?

A function of $x$ is right-sided if and only if there is some
$x_0$ such that the value of the function is zero for all $x <
x_0$, and left-sided if and only if there is
some $x_0$ such that the value of the function is zero for all $x >
x_0$.

Likewise for signals with time represented by $x$. A double-sided signal has neither restriction. For example a Gaussian function and $\sin(x)$ (both pictured below) are double-sided.

$f(x) = \frac{e^{-\frac{x^2}{2}}}{\sqrt{2 \pi }}$

$f(x) = \sin(x)$


Removing high pitch noise embedded in audio spectrum

Can any body tell how to remove it completely so that it is not hearable?

As suggested already, the signal has harmonics.

You can diminish it further by keep filtering the rest of the peaks of the spectrum with a notch filter, or using the noise reduction option in Audacity.

The latter option is better when the interference is over a useful part of the recording. To apply that filter, select a part of the recording that is representative of the noise profile you are dealing with. Audacity will use this to extract the noise spectrum and then create a filter that cancels it. You can then take that profile and apply it to other points in the recording to reduce its impact.

  It will be bonus for me if someone can tell why this high pitch noise occur in audio?

It sounds like interference and it is likely that it is picked up as a square wave. It definitely doesn't sound like a hum, more like a high pitch pipe organ sound with two or three components.

It might be coming from unshielded cables, strong sources nearby an amplifier (e.g. WiFi router / dongle, LED lightstrip, TV, other) or if it is an audio USB dongle, of dubious built quality, it might be picking up interference from the data exchange going on on the bus. It's hard to tell without knowing what is this a recording of (?).

Hope this helps.


Selecting gradient orientation for Generalized Hough Transform while scaling and rotation is present

If the shape is rotated by $\theta$, then the gradient orientation ($\phi$) for a given edge point changes. So, shouldn't we do either one of the following:

  Rotate all the ϕ values in R table by θ?

  OR

  Rotate gradient vector by ($-\theta$) and then calculate $\phi$ for the edge point?

The generalised Hough transform does #1 implicitly by the $S, \theta$ parameters when accumulating.

Note here that usually, you don't really know that the extracted contour has been rotated by some angle (let's call it $\alpha$). You usually are trying to find that angle.

During the preparation stage, you construct the "template", $\phi$. This quantises the circumference of a circle that encloses the shape to be recognised into $k$ levels ($k$ distinct angles). Then, for each $k^{th}$ level, you look at the contours that make up your template image at the $\phi_k$ bearing and note down their $(r, \beta)$. The $(r,\beta)$ is nothing more than the polar coordinates of an element of your template image's contour that is oriented at angle $\phi_k$.

If you wanted to recognise a simple circle, then you would only have one column of $(r, \beta)$. If that circle enclosed another, smaller circle (imagine for example a road sign), then at some orientation $\phi_k$ you would have more than one "edges" and therefore more than one columns of $(r, \beta)$ data.

Now, during the recognition phase, you use the $\angle G$ to lookup an entry in $\phi$. Once you have found that entry, your take that entry's $(r, \beta)$ and use them to find the "right Hough counter cell" $H(x_c,y_c)$ and increment it.

The point here is that the $(r, \beta)$ are calculated with respect to the centroid and it is that centroid $x_c,y_c$ that accumulates quickly in $H$. When the shape is matching, it means that your test image contains a lot of $\angle G$ entries that have resulted in the $H(x_c,y_c)$ counter getting a very high value.

The next step is to be able to recognise the template $\phi$ irrespectively of scale $S$ and orientation $\theta$.

Notice here that if you were to take the Hough transform of the template with its centroid translated at $(0,0)$, then different rotations $\theta_1, \theta_2, \theta_3, \ldots$ are basically cyclic lateral shifts in the Hough space (where the $x$ dimension now corresponds to bearing $\theta$) and different (uniform) scalings $S_1, S_2, S_3, \ldots$ correspond to different "arc amplitudes" (in the $y$ direction that corresponds to distance from origin).

So, the template doesn't really change, all that you do is consider different scalings and rotations at the same time. From another point of view it is as if you were running $N$ different Hough transforms where you scale / rotate the original image by $(S,\theta)_n$. You would produce one $H$ for every $(S, \theta)_n$ such as $\left(1,0\right), \left(1, \frac{\pi}{4}\right), \left(1, \frac{\pi}{2}\right), \left(1, \frac{3 \pi}{4}\right), \left(1, \pi\right), \left(0.5,0\right), \left(0.5, \frac{\pi}{4} \right), \ldots, \left(0.25, 0 \right), \ldots$. In this case here we are looking for 4 orientations at scale 1, 0.5, 0.25 and so on.

But, rotating and scaling the original image can be costly because there are definitely more pixels in an image than elements in the $k$ level quantised periphery of the curve you are trying to find.

For this reason, instead of rotating the original image and running parallel Hough transforms, you simply scale and rotate $\phi$ but you do it on the fly. Instead of modifying the $\phi$ you just introduce $S, \theta$. $S$ modulates your $r$ (the radius of the template) and $\theta$ offsets your $\beta$ (the bearing of your template).

In this way, you produce $N$ matrices $H(x,y)$, one for each $(S, \theta)_n$ you would like to consider. Again, the $x_c, y_c$ counter that crosses a certain threshold marks the location of the recognised contour but in addition now you can also infer at which scale and rotation by fixing the $(S, \theta)$ pair.

Hope this helps.


How does this represent a low-pass filter?

This is one tortured piece of code. I strongly recommend to toss it and write whatever you need from scratch.
Arguments

N filter order or some design parameter
fft_size: size of the final vector
low_cut: the actual filter impulse Response padded to FFT size

Code

First loop creates a cosine window of length N and bias it with 0.5
2nd loop zero pads to FFT length
3rd loop calcuates the sum of the amplitudes
4th loop divides by the sum and multiplies with -1.
5th loop seems to time reverse the filter (pointless since it's symmetrical) AND copies towards the end of FFT buffer AND also writes past the end of the FFT buffer in the process.
6th loop shift N samples forward by N/2-1 overwriting what was previously in there

What happens exactly here depends on values of N and fft_size and there is some minor low cut behavior in there, but overall it seems  like and extremely convoluted, imprecise and dangerous way of doing a low cut


Is it appropriate to use PSD on the analysis of a plucked guitar string?

I'd really recommend going for a short-time fourier transform such as that produced by: matplotlib.mlab.specgram. It gives you the change in the spectrum over time and can be calculated using PSD or the magnitude spectrum.

The reason I'd recommend this is that your guitar signal is not spectrally stationary; your graphs are only a snapshot at one point in time. What gives your guitar string their tonal differences is the way all the harmonics shown above change over the length of the note.

There's a very good summary of PSD vs FFT magnitude here and here.


Find camera pose by knowing triangle coordinates

No.  No.  You can't.
Not completely at least.  With only three points, the math has two solutions.  Your triangle could be tilted in one of two positions.  Like picking a square root.
From a single frame, you might be able to tell and measure motion from a blur gradient, but that is difficult.
From multiple frames, calculating the position of the camera for each frame, the velocity is simply a matter of differentiation techniques on a finite sequence.


2nd Order Type 2 PLL for Tracking Frequency Error

Notation: Incoming carrier frequency : $f_c$ and Local Oscillators synthesized frequency : $f_o$.   A First Order PLL to track a constant phase offset $(\phi)$ in the received carrier sinusoidal is based on the following PLL daigram. In working of the single PLL Loop it is essential that $f_c = f_o$:  (image courtesy : Software Receiver Design by Johnson, Sethares & Klein) Here, $r_p(kT_s) = BPF_{2f_c} \{ r(kT_s)^2 \}$, meaning received signal is squared and then bandpass filtered with BPF's center frequency at twice carrier frequency $2f_c$. This preprocessing is done to emphasize the carrier, since this will boost carrier SNR. So, input to the PLL is basically, $$r_p(kT_s) = cos(4\pi f_c k T_s + 2\phi)$$Now if instead of constant phase offset of $\phi$ if we are having a constant frequency offset of $\Delta f = f_c - f_{o}$ and a constant phase offset of $\phi$ in carrier frequency then, we can extent this single PLL Loop to a combination of 2 PLL loops as follows: (image courtesy : Software Receiver Design by Johnson, Sethares & Klein)  Here, input to this frequency Tracking PLL is as follows: $$r_p(kT_s) = cos(4\pi (f_c - \Delta f)k T_s + 2\phi)$$  Now, if the first Loop processing is worked out, it can be shown that it tracks a ramp:$$\theta_1[k] \rightarrow 2\pi (f_c - f_{o})k + \beta$$, where $\beta$ is the y-intercept of the ramp being tracked by first PLL.This is added to $\theta_2[k]$, which makes it possible for the $2^{nd}$ PLL to track the frequency accurately :$$\theta_2[k] \rightarrow \phi - \beta$$ , which mean, $$\theta_1[k] + \theta_2[k] \rightarrow 2\pi (f_c - f_{o}) + \beta + \phi -\beta$$ $$\theta_1[k] + \theta_2[k] \rightarrow 2\pi (f_c - f_{o}) + \phi$$, And this sum is being used to drive the sinusoid of Local Oscillator, hence the local oscillator is now generating a sinusoid which is synchronized with $r_p(kT_s)$ in both frequency and phase. Basically, top PLL Loop is tracking a ramp and hence generating the correct frequency for the bottom loop, and bottom loop is further correcting the residual phase error.(image courtesy : Software Receiver Design by Johnson, Sethares & Klein)


Help or suggestions with Line detection in microscopy images

I used Median Blur instead of Gaussian Blur, and used Sobel Filter to detect the edges, the results are not perfect, but you can improve them if you tuned the parameters well(and they are many!)

#==================
# Import Libraries
#==================
import numpy as np
import cv2
import matplotlib.pyplot as plt

# Read Image
img = cv2.imread('1.png', 0)

# Median Blur
img = cv2.medianBlur(img, 3)

# Sobel filter for Edge detection
sobelx = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=5)
sobely = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=5)
sobel = np.sqrt(sobelx**2+sobely**2)

# Threshold Sobel to locate the lines
maxi = 0.30*sobel.max()
sobel[sobel<maxi] = 0

# Normalize
sobel = (255*sobel) / (sobel.max()-sobel.min())
sobel = sobel.astype('uint8')

# Hough Line transform.
rho = 1
theta = np.pi / 180
threshold = 15
minLineLenght = 18
maxLineGap = 13
myBlankLines = img.copy()
lines = cv2.HoughLinesP(sobel, rho, theta, threshold, np.array([]),
                    minLineLenght, maxLineGap)
for line in lines:
    for x1,y1,x2,y2 in line:
        cv2.line(myBlankLines,(x1,y1),(x2,y2),(255,0,0),1)

# Visualize Results
plt.imshow(myBlankLines, cmap='gray')
plt.axis('off')
plt.show()


Relation between the matrix trace and the amplitude of each element

The trace of a matrix just scales with the scaling of the matrix:
$$\textrm{tr}\{c\mathbf{A}\}=c\,\textrm{tr}\{\mathbf{A}\}\tag{1}$$
So if you scale the vector $\mathbf{p}$, the trace scales with the squared magnitude of the scaling constant. With $\mathbf{p}=c\mathbf{p}_0$ you have
$$\begin{align}\textrm{tr}\{\mathbf{Ap}\mathbf{p}^H\mathbf{A}^H\}&=\textrm{tr}\{\mathbf{A}c\mathbf{p}_0c^*\mathbf{p}_0^H\mathbf{A}^H\}\\&=\textrm{tr}\{|c|^2\mathbf{A}\mathbf{p}_0\mathbf{p}_0^H\mathbf{A}^H\}\\&=|c|^2\textrm{tr}\{\mathbf{A}\mathbf{p}_0\mathbf{p}_0^H\mathbf{A}^H\}\end{align}$$


how to add two images in one, next to each other

combinedImage = [leftImage, rightImage];


Determining the mean of of a sampled signal (improvements to time averaging method)

Edit: After thinking about it some more I came up with something much simpler than the phase-locked loop.

The problem you are having is because you are filtering with a boxcar.  The boxcar filter has a lot of ripples in the frequency domain, so if you choose the wrong width you don't get good attenuation of your approximately 10Hz signal.

If you use a Butterworth filter you will get a frequency response that has no ripples.  Everything over (say) 5Hz will be attenuated by more than -3db.  A Butterworth filter is also cheap to calculate.

I went to http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html and asked for a 3rd order Butterworth low-pass filter with a cutoff frequency of 5Hz on a sample rate of 1000Hz and got the following recurrence relation:

y[n] = (  1 * x[n- 3])
     + (  3 * x[n- 2])
     + (  3 * x[n- 1])
     + (  1 * x[n- 0])

     + (  0.9390989403 * y[n- 3])
     + ( -2.8762997235 * y[n- 2])
     + (  2.9371707284 * y[n- 1])

You'd probably get a better result by using Matlab to design an elliptic filter (and probably one with higher order than 3).  You'll get sharper attenuation beyond the cutoff frequency.

Here is my original answer about phase-locked loops:

I would try a discrete time phase-locked loop.

So something like this:

The idea is to multiply your input signal by a sinusoid of the same frequency, followed by a low pass filter.  This shifts the sinusoidal part of your input signal to (nearly) zero frequency.  By observing the changes to the output of the low pass filter you are getting an estimate of how far off your estimated frequency is from the actual frequency.  So you feedback the adjustment.

The output of the multiplier in the picture above is:
$$
Be^{2\pi i \theta_n} + Ae^{2 \pi i (\theta_n + \omega n)} + Ae^{2 \pi i (\theta_n - \omega n)} + \varepsilon'(n).
$$
$\theta_n$ is the running sum of the estimated frequency, $\hat{\omega}_n$.  When $\hat{\omega}_n = \omega$ then $\theta_n-\omega n$ is constant, so $A e^{2 \pi i (\theta_n - \omega n)}$ is constant.

The derivative you want here is the angular derivative.  Given samples from the low pass filter, $r_ne^{2\pi i \alpha_n}$ and $r_{n+1}e^{2\pi i \alpha_{n+1}}$ you really want $\alpha_{n+1}-\alpha_{n}$.  Dividing low-passed sample $n+1$ by low-passed sample $n$ and taking the complex part gets you essentially the right thing.

The "adjust $\hat{\omega}_n$" box is also non-trivial.  The estimated phase error is likely to be a little noisy so you might want to low-pass filter the phase error and you might want to only adjust $\hat{\omega}_n$ by a fraction of the estimated phase error, rather than by the entire phase error.

$\hat{\omega}_n$ is the estimated frequency, and $1/\hat{\omega}_n$ is the low-pass filter width you are looking for in your question.  I think you can apply the same trick to the low-pass filter in the phase-locked loop.  By dynamically adjusting the phase-locked loop low-pass filter width to $1/\hat{\omega}_n$ you completely eliminate $Be^{2\pi i \theta_n}$.


number of possible component in sinusoidal model

I guess that your problem can be modeled as a harmonic retrieval problem.
You can check the Kung's method.


About Discrete Fourier Transform vs. Discrete Fourier Series

okay, i'm gonna expound a little.
quoting (except for any typos that may result) from the 1989 text of O&S (Introduction to Chapter 8, The Discrete Fourier Transform, p 514):

Although several points of view can be taken toward the derivation and interpretation of the DFT representation of a finite-duration sequence, we have chosen to base our presentation on the relationship between periodic sequences and finite-length sequences. We will begin by considering the Fourier series representation of periodic sequences.  While this representation is important in its own right, we are most often interested in the application of Fourier series results to the representation of finite-length sequences. We accomplish this by constructing a periodic sequence for which each period is identical to the finite-length sequence. As we will see, the Fourier series representation of the periodic sequence corresponds to the DFT of the finite-length sequence. Thus our approach is to define the Fourier series representation for periodic sequences and to study the properties of such representations. Then we repeat essentially the same derivations assuming that the sequence to be represented is a finite-length sequence. This approach to the DFT emphasizes the fundamental inherent periodicity of the DFT representation and ensures that this periodicity is not overlooked in applications of the DFT.

section 8.1, p 516 on the DFS:

Eq. (8.11) $\quad \tilde{X}[k] = \sum\limits^{N-1}_{n=0} \tilde{x}[n] \ e^{-j2\pi n k/N} $

Eq. (8.12) $\quad \tilde{x}[n] = \frac{1}{N} \sum\limits^{N-1}_{k=0} \tilde{X}[k] \ e^{+j2\pi n k/N} $

regarding the DFS, $\tilde{x}[n]$ (with the tilde) is defined to be periodic with period $N$ such that $$ \tilde{x}[n+N] = \tilde{x}[n] \quad \forall n $$ and $\tilde{X}[k]$ turns out to also be periodic with period $N$ (so $ \tilde{X}[k+N] = \tilde{X}[k] \quad \forall k $)
later, in section 8.6, p 532 on the DFT:

Eq. (8.59) $\quad X[k] = \begin{cases} \sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\ 0, & \text{otherwise} \end{cases} $

Eq. (8.60) $\quad x[n] = \begin{cases} \frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 0, & \text{otherwise} \end{cases} $

Generally the DFT analysis and synthesis equations are written as

Eq. (8.61) $\quad X[k] = \sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N} $

Eq. (8.62) $\quad x[n] = \frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N} $

In recasting Eqs. (8.11) and (8.12) in the form of Eqs. (8.61) and (8.62) for the finite-duration sequences, we have not eliminated the inherent periodicity. As with the DFS, the DFT $X[k]$ is equal to samples of the periodic Fourier transform $X(e^{j\omega})$, and if Eq. (8.62) is evaluated for values of $n$ outside the interval $0 \le n \le N-1$, the result will not be zero but rather a periodic extension of $x[n]$. The inherent periodicity is always present. Sometimes it causes us difficulty and sometimes we can exploit it, but to totally ignore it is to invite trouble.

so the first obvious thing i would say is that the tildes used for the DFS (to explicitly depict a periodic sequence) are symbols and still do not change any mathematical fact.  the direct relationship between the periodic $\tilde{x}[n]$ and the "finite-length" $x[n]$ is
$$ \tilde{x}[n] = x[n \bmod N] \qquad \forall n \in \mathbb{Z}, \ N \in \mathbb{Z}>0$$
where $ \qquad\qquad\qquad n \bmod N = n - N \left\lfloor \frac{n}{N} \right\rfloor  $
now i know some folks will point to the Eqs. (8.59) and (8.60) definition of the DFT that has truncated (to $0$) values outside of the interval $0 \le n,k \le N-1$.
however, that definition is contrived.  it could just as well be expressed as
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
5, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 5, & \text{otherwise}
\end{cases} $
or
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
5000, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 5000, & \text{otherwise}
\end{cases} $
or
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
\text{the man on the moon}, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ \text{and his hot girlfriend}, & \text{otherwise}
\end{cases} $
because that $0$ in that contrived DFT definition will never ever be used in any theorems regarding the DFT.  when that contrived definition is used for the DFT, then when using any DFT theorems to do any real work (other than the linearity and scaling by constant theorems), then one must use modulo arithmetic in the arguments of $x[n]$ or $X[k]$.  and using that modulo arithmetic is explicitly periodically extending the sequence.
so (sorta responding to hotpaw) there are two or three processes that you should think about when using the DFT on a real signal.

the sampling process. what happens to the spectrum of $x(t)$ when you sample it with a "dirac comb" or whatever you want to call the sampling function?

windowing to finite length. what happens when you window either $x(t)$ or the sampled version, $x[n]$, with a rectangular window of length $N$?

periodic extension.  what happens when you periodically extend it by repeatedly shifting the windowed $x[n]$ by $N$ samples and overlap and add it?

deal with each step by itself.


FIR design for arbitrary magnitude and phase: Why can I separate real and imaginary parts like this?

Neither of you is wrong. Your final expression is just the same as $H(\omega)\cdot e^{-j\omega(N-1)/2}$, which is (an approximation of) the desired complex frequency response with an additional linear phase term that depends on the chosen filter length $N$. This means that you need to subtract that linear phase term from the desired phase response, because it will be added in the design process.

Let the desired complex frequency response be given by

$$D(\omega)=|D(\omega)|\cdot e^{j\phi(\omega)}\tag{1}$$

Now you have to compute a new complex response by subtracting the linear phase term $-\omega (N-1)/2$ from the phase of $(1)$:

$$\tilde{D}(\omega)=|D(\omega)|\cdot e^{j[\phi(\omega)+\omega(N-1)/2]}\tag{2}$$

You have to use the real and imaginary parts of $\tilde{D}(\omega)$ to design the filter.

What is meant by the option 'hilbert', is that you tell the design routine that the desired response (the imaginary part of $\tilde{D}(\omega)$) should actually be multiplied by $j$ and it is an odd function of frequency, because otherwise it will be interpreted as a purely real and even function. This will result in an asymmetric impulse response. In principle this can be done, but I don't know if it can actually be done for arbitrary prescribed responses with the current versions of Matlab's filter design routines.

Note that such an approach is actually not really necessary for the least squares design of FIR filters, because solving a complex least squares problem (for non-linear phase designs) is just as simple as solving a real least squares problem (for linear phase designs). As far as I know it's just not implemented in Matlab's Signal Processing Toolbox. I've written a few Matlab/Octave functions that solve the complex least squares approximation problem for FIR filter design: lslevin.m for complex least squares design of real-valued filters, and cfirls.m for filters with complex coefficients.

For equi-ripple designs, the complex (non-linear phase) approximation problem is indeed much harder than the real (linear phase) problem. For the latter we have the Remez exchange algorithm, which was adapted to the linear phase FIR filter design problem by Parks and McClellan and which is implemented in Matlab (firpm.m). There are algorithms for solving the complex equi-ripple design problem, but they are much less efficient than the Parks-McClellan algorithm. One of them is implemented in Matlab (cfirpm.m).

Note that for the least squares approximation, solving two real-valued approximation problems for the real and imaginary parts gives the same solution as solving one complex least squares approximation problem. This is easily seen as follows. The complex approximation error is given by

$$E(\omega)=D(\omega)-H(\omega)\tag{3}$$

where $D(\omega)$ and $H(\omega)$ are the desired and actual complex frequency responses. The complex least squares solution minimizes the integral over $|E(\omega)|^2$, whereas the individual real-valued least squares solutions minimize the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$, where $E_R(\omega)$ and $E_I(\omega)$ are the real and imaginary parts of $E(\omega)$, respectively. Since

$$|E(\omega)|^2=E_R^2(\omega)+E_I^2(\omega)\tag{4}$$

minimizing the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$ is equivalent to minimizing the integral over $|E(\omega)|^2$.

For the complex Chebyshev approximation problem the situation is different. The quantity that is minimized is the maximum of the magnitude of the complex error over the bands of interest $\max_{\omega}|E(\omega)|$, but since

$$\max_{\omega}|E(\omega)|\neq \max_{\omega}|E_R(\omega)|+\max_{\omega}|E_I(\omega)|\tag{5}$$

independent Chebyshev approximation of the real and imaginary parts is not equivalent to complex Chebyshev approximation. If for a given complex approximation problem the maximum error of the optimal solution equals $\delta$, then in the worst case the errors of both real-valued approximations will also be equal to $\delta$. This means that in the worst case, the maximum complex error of the solution obtained from independent optimization of the real and imaginary parts becomes

$$\max_{\omega}|E(\omega)|=\max_{\omega}\sqrt{E_R^2(\omega)+E_I^2(\omega)}\le\sqrt{\max_{\omega}E_R^2(\omega)+\max_{\omega}E_I^2(\omega)}\le\sqrt{2}\delta\tag{6}$$

This is equivalent to a maximum of $3$ dB increase in error compared to the optimal solution.

I'll show an example design to illustrate a few points that I made. I designed a bandpass filter with a desired phase response that is linear but with a smaller group delay than an FIR filter with perfectly linear phase. The filter length is $N=61$, so the delay of a linear phase FIR filter would be $(N-1)/2=30$ samples. I chose the desired delay to be $20$ samples instead. I used a complex least squares and a complex Chebyshev criterion to approximate the given complex frequency response. For each criterion I designed the filter in two different ways. First, directly in the complex domain, and second by using the real and imaginary parts of the desired frequency response to design two linear phase filters, which, after combining them, should approximate the original complex frequency response. The figure below shows the design results. As expected, for the least squares design, both solutions are equivalent (up to numerical accuracy) and only one of them is shown in the plots (in blue). For the Chebyshev criterion both solutions are slightly different, as predicted. The solution from splitting the response into real and imaginary parts (in red) has a larger error by about $1.2$dB (as can be seen at the first sidelobe in the stopband). In the passband all three solutions are very similar. Note that all three filters approximate the desired group delay value of $20$ samples relatively well over a large portion of the passband.


I like this paper a lot: The philosophy of PCM. It explains PCM very clearly, includes references to some of the earliest papers on it, and you can hardly argue with the qualifications of the authors :)
PCM was conceived as a three-step process. First, sampling (at frequency $f_s$ sufficiently above Nyquist). Second, quantizing (with enough resolution $\Delta$ to obtain the required quantization SNR). Third, coding, or mapping the quantized samples to $M=2^k$ pulse amplitudes, with $M$ chosen to obtain the required line pulse rate.
PCM, then, has parameters $f_s$, $\Delta$ and $M$. Each of these affect system performance in different ways.
Why use the term "modulation"? In communications, "modulation" is very general, and is used almost as a synonym to "mapping" or "converting elements from an input set to an output set". In fact, what you understand as "modulation" (multiplying by a carrier) is often called "upconversion" and "downconversion", to avoid confusion.


Multi Scale Shape and Detail Enhancement from Multi Light Image Collections

What the authors meant is to create a matrix of Weights, $ {U}^{\left( i, j \right)} $.

It is a matrix of the size of the image.
The given calculation is by the exponent of two terms (Each of them is a matrix, the calculation is element wise).
The final step is to blur it using a Gaussian Blue (2D Blur, like on an image).

In equation 10 you can see they take care of the normalization of the weights.

It would be nice if you come back and show us your results.


It's something like "finite energy signals"

$$ \int_{-\infty}^{\infty}|x(t)|^2 dt<\infty $$

$$ \sum_{n=-\infty}^{\infty}|x[n]|^2 <\infty $$

except that the "finiteness" is a tougher constraint with

$$ \int_{-\infty}^{\infty}|x(t)| dt<\infty $$

$$ \sum_{n=-\infty}^{\infty}|x[n]| <\infty $$

for example, the signal

$$ x(t) = \frac{1}{\sqrt{1+t^2}} $$

is finite energy, but is not finite action.


Initial rest condition applied on $x(t)$ vs $h(t)$

It should be clear that a property of a system, such as causality, cannot be determined by looking at its input signals. For a linear time-invariant system, it is its impulse response $h(t)$ from which properties such as causality or stability can be determined.
Only the second definition in the question is correct: a causal LTI system has an impulse response $h(t)$ which equals zero for $t<0$. This implies that the system's output at time $t_0$ only depends on values of the input signal $x(t)$ for $t\le t_0$, and not on future values $x(t)$, $t>t_0$.
The output signal of a causal LTI system can be written as
$$y(t)=\int_{0}^{\infty}h(\tau)x(t-\tau)d\tau=\int_{-\infty}^tx(\tau)h(t-\tau)d\tau$$
An initial-rest condition just means that if the input $x(t)$ is zero for $t<t_0$ then the output $y(t)$ must also be zero for $t<t_0$.


Your first equation is indeed valid for all $t$, but this has nothing to do with causality. You have

$$\begin{align}v(t)&=\frac{1}{C}\int_{-\infty}^ti(\tau)d\tau\\&=\frac{1}{C}\int_{-\infty}^{t_0}i(\tau)d\tau+\frac{1}{C}\int_{t_0}^ti(\tau)d\tau\\&=\frac{1}{C}\int_{-\infty}^{t_0}i(\tau)d\tau-\frac{1}{C}\int_{t}^{t_0}i(\tau)d\tau\end{align}$$

If $t_0>t$ you simply subtract from the first integral the contribution from the interval $[t,t_0]$. But this is still causal, you don't need the future input, you just need to subtract it in case you added it in the first integral. The output could of course be computed without it.

Concerning your second question, you choose the time $t_0$ as the beginning of your processing. The value of $t_0$ doesn't matter but once you choose it, it's fixed. The system is linear if the initial conditions at $t=t_0$ are zero.


Equivalence of Convolution Through Polynomial and Matrix Products

Convolution and polynomial multiplication are equivalent by definition. Offsets are usually introduced through indexing. The treatment of ends is important.


I saw same issue, and I found the answer here.
for short, you can set imaginary part of filtering kernel zeros before you call mulSpectrums.

Mat planes[] = {Mat::zeros(dftImage1.size(), CV_32F), Mat::zeros(dftImage1.size(), CV_32F)};
Mat kernel_spec;
planes[0] = your_filter; // real
// planes[1] = your_filter;
// imaginar should be zero! Don't touch it.
merge(planes, 2, kernel_spec);
mulSpectrums(...

Hope this helps.


Is concept of "bit" in computer programming similar to the concept of "bit" in information theory?

They are not the same, but they're related. In particular, if you look at a computer memory holding $M$ "computer" bits, where each bit can be considered random and independent of all other bits, and there are roughly 50% of zeros, then the memory also holds roughly $M$ "information theory" bits.

Of course, this is often not the case: computer bits are usually correlated, and not uniformly random. This is why they can be compressed. Compressor programs such as LZW ("source coders" in information theory parlance) work, in a sense, by making each computer bit hold one information bit.

Edited to add: This example may make the distinction clearer. Consider a memoryless source with two outputs, $m_1=000$ and $m_2=001$, with probability 0.5 for each. Clearly, the information in each message is one (information) bit, but its length is three (computer) bits. A source coder, such as the Huffman algorithm, will readily code the messages to $c_1=0$ and $c_2=1$, compressing the source output. You can easily extrapolate this example to a source that produces ASCII-encoded text.

Note that, in the case of written languages in general and English in particular, nobody knows what the actual source entropy is, because there is no model for it. That is why there are contests for the best compression of large bodies of text; nobody is really sure what the optimum compression algorithm for English is.


use prony method for damped exponentials

you have C and U, and C can be written as C_i = (A_i/2)*exp(j*phase)
U_i can be written as exp(damp_i + j*2*pi*f_i)*T
its should be easy to extract the info from that. T is sampling period.

%to get damping factor i and frequency i use something like this
%assume damping factor at sample 1 is 0.3, frequency at 1 is 20
%and out sampling frequency is 500 Hz
damp1 = [0.3 0.4 0.5];
f1 = [20 30 40];
T = 1/500;
%assume this is your ui vector from th ealgorithm
ui = exp((damp1 + j*2*pi*f1)*T);
%get the damping and frequencies
lv1 = log(ui);
r_damp1 = real(lv1)/T;
r_f1 =imag(lv1)/(T*2*pi);

%assume this is your Ci vector from algorithm
Ai = [0.1 0.2 0.3];
phasei=[pi/8 pi/4 pi/2];
Ci = (Ai/2).*exp(j*phasei);
%extract phase and amplitude
r_Amp = abs(Ci)*2;
r_phase = angle(Ci);


Checking the "cleanliness" of an image for Optical Character Recognition

Your images are gray-scaled. In order to perform the OCR you might want to first threshold it to apply OCR on the binary image.

The effectiveness of the thresholding will definitely affect the performance of the algorithm. A simple way is say Otsu's thresholding. In your case, you can think of this as character vs. non-character pixels. And the measure of the cleanliness can be treated as the amount of energy (such as class variance) of the non-character pixels.

Unfortunately the measure I described is post applying your (first step) of OCR. Hence, I am not sure if this helps you.


The mistake could be in assuming you know which angle $a$ to use when you don't know if the source is in front of or behind the array yet. In other words, you need to solve the same problem for $a$ that you're solving for $b$.

For your approach to work, you need to make sure side $a$ is the closest side to the source. You could find the closest side by finding the pair of mics where the maximum arrival time of the pair is less than the minimum arrival time of the opposite pair. This beam points outward, and the beam on the opposite side points inward. Only one pair of mics will meet those conditions at a time.

Also note that your approach for choosing which angle to use for $b$ only works if you consider the whole array as a single point (the distance between the mics relative to the distance to the source is vanishingly small). Otherwise, if the source is between the middle of side $a$ and one of $a$'s mics on the same axis as side $a$, angle $b$ will point the wrong way.

Instead of doing it like that, you might look at solutions using more than two mics at once.

An approach recently suggested here uses three equidistant mics and takes the cross product of the vectors pointing from one mic to another, with vectors having an additional component for the time difference.

Here's a really simple solution based on that, using axis-aligned diagonal mic pairs. Set your mics up like this:

      4
    1   2
      3

The time delta between mics 1 and 2 is the X component of a vector towards the source, and the time delta between mics 3 and 4 is the Y component. Check it out here.


Spectrum analyser and resolution bandwidth

This document is seriously awesome.  It gave me a good understanding of how spectrum analysers work in an easily understandable way.

http://literature.agilent.com/litweb/pdf/5965-7920E.pdf


Finding amplitudes from FFT data of a 2D image

If you can adjust the orientation of your image so that the sinusoid pattern runs vertically to the image, then you can obtain a vertical slice, perhaps at the middle of the image and use that as a 1D waveform for further analysis.

If the waves can come out of your experiment in any orientation you basically have two options:

1) Estimate the angle of rotation ($a$), apply a rotation by $-a$ and then select a slice of the image (just as in the introduction). To estimate the angle of rotation you can look at the output of the Hough Transform which will give you strong "spots" at a given angle for each line (or line-looking) form within your image. That would be the elongated folds of the wrinkles in your image. For more information please see: http://en.wikipedia.org/wiki/Hough_transform

2) Do a straightforward 2D Fourier Transform to recover both orientation and amplitude spectra (at least for the example you provide). The sinusoidal-like form will appear as a linear feature of some orientation in your FFT spectrum which you could pick up with a simple thresholding operation. For more information please see http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf

Hope this helps.


GNU Radio | Packet Encoder/Decoder | Unwanted bytes added to the transmitted stream

Packet Encoder and Decoder are broken; they drop data.

That's why they are in the deprecated category (for years now!). We've removed them, because as a project, GNU Radio has not been able to fix them (and also, they were terrible from an architecture point of view).

So there's exactly one solution: don't use packet encoder / decoder.


Statistical significance of coherence values

Great question! I've been grappling to determine how to calculate statistical significance of coherence for several days. This question helped me find some key references I hadn't previously located!

One detail regards whether you're talking about "coherence" or "magnitude-squared coherence". The first resource uses the former and the second uses the latter. I generally prefer the latter:

$$C_{x,y}(\omega)= \frac{|\langle{X^*(\omega)Y(\omega)}\rangle|^2}{\langle|X(\omega)|\rangle^2\langle|Y(\omega)|\rangle^2}$$

Where $X(\omega)$ is the fourier transform of $x(t)$, and $\langle \rangle$ denotes an ensemble average of $n$ independent ensembles.  Obviously the plain "coherence" is the square root of what is defined above.

The significance level of this magnitude-squared coherence as defined by the two resources is:

$$L_1(\alpha, q) = 1-\alpha^{1/q}$$
$$L_2(\alpha, q) = F_{2,2q}(\alpha)/(F_{2,2q}(\alpha)+q) $$

Where $\mathrm{DOF} = 2n = 2 ( q + 1 )$, and $F$ is the inverse-CDF (percent point function) of the F-statistic. Also, just FYI, the second formulation can also be found in Priestley 1981 page 706.

The following Python code implements the two approaches for $\alpha=0.05$ (i.e. a 95% significance level):

import numpy as np
import scipy.stats as st

alpha = 0.05
p = 1 - alpha

dof = np.arange(5, 100)
n = dof / 2.

C_1 = 1 - alpha ** (1. / (n - 1))

fval = st.f.ppf(p, 2, dof - 2)
C_2 = fval / (n - 1. + fval)

Plotting the two functions shows that they are apparently identical:

import matplotlib.pyplot as plt
ax = plt.gca()

ax.plot(dof, C_1, 'k-', lw=10, label=r'$1-\alpha^{1/q}$')
ax.plot(dof, C_2, 'r-', lw=5, label=r'$F_{2,2q}(\alpha)/(F_{2,2q}(\alpha)+q)$')
ax.plot(dof, 6. / dof, 'b-', lw=2, label=r'$6/d$')

ax.set_xlabel('Degrees of Freedom')
ax.set_ylabel('{:0.0f}% Significance Level'.format(100 * p))
ax.legend()
ax.set_title(r"Coherence level for $\alpha=0.05$")

I compared $L_1$ (black) and $L_2$ (red) for several different values of $\alpha$ and they always agree (their difference was always $O(10^{-16})$). I didn't bother doing the algebra to show that $L_1=L_2=L$, but perhaps someone else will find that to be an interesting undertaking?

Lastly, I've included the approximation: $L(\alpha=0.05)\approx6/\mathrm{DOF}$ (blue). I got this from U. Washington Applied Physics Lab professor Jim Thomson, which he apparently got from legendary oceanographer Carl Wunsch. It appears to be a very good approximation when $\alpha=0.05$ and $\mathrm{DOF} > 20$ or so. I'd be interested to see the algebra that shows that:

$$\lim_{\mathrm{q}\rightarrow\infty}L(1/20,q) = \frac{6}{\mathrm{DOF}}$$


How Butterworth low-pass filter can be applied on a digital signal (i.e. discrete representation of a signal)?

What you need to do is find a digital filter $H(z)$ that behaves similar to the $H(s)$ of your Butterworth filter. That is, you need to read about IIR implementation.

There are several methods described in signal processing books about this subject. I recommend you to check for example "Essentials of Digital Signal Processing" from Lathi and Green.

Two of the most common methods are the "Impulse invariance" and the "Bilinear transform".

The impulse invariance requires you to find the impulse response $h(t)$ of your $H(s)$. After doing this, you sample it replacing $t$ by $nT$, to finally Z-transform your $h[n]$ into $H(z)$. Bare in mind that this method only works for band-limited filters (lowpass and bandpass) because of aliasing issues. You should determine an appropriate small enough T to avoid this.
The bilinear transform, on the other hand, consists on applying the following mapping of your variable $s$ into $z$:

$s=\frac{2}{T}(\frac{z-1}{z+1})$

This method has a drawback as well, which is the frequency warping, you can read about it here https://en.wikipedia.org/wiki/Bilinear_transform#Frequency_warping, but is in general easy to counter by pre-warping your original $H(s)$.


Element-wise multiplication in the spatial/temporal domain with a windowing function serves to reduce the effect of the potentially large jump you get at the signal edge, when making the sampled signal periodic. This jump would otherwise introduce lots of frequencies that might not be present in the signal.

Element-wise multiplication in the frequency domain with a windowing function applies a low-pass filter. Multiplying by a Gaussian window is equivalent to applying a Gaussian filter in the spatial/temporal domain.

If the code you found online were to compute the inverse transform after applying the window, then whoever wrote it would be low-pass filtering. With the code as you posted, I tend to think that person doesn't know what they're doing (or maybe it's just a typo).


Pattern matching methods: Based on Area and Based on features

Yes, if applied properly, HOG is a good feature extractor even for OCR. I could point you towards an OpenCV sample and a MATLAB digit classification tutorial which do exactly that. However, if you like to get state-of-the-art performance, I would suggest deep learning methods, more specifically convolutional neural networks. Many works demonstrated the power of these approaches on MNIST dataset. You might like to check this one or this one for some beginner's tutorials.


How we can quantize a sampled signal in MATLAB?

First you will need to determine the number of quantization levels. I am going to assume a power of two for digital convenience's sake.

 nbits = 8 % 256 qantization levels
 qLevels = 2^nbits

The next step will be to scale your signal to have the same magnitude as your number of bits.

 signalMin = -1
 signalMax = 1
 scalingFactor = (signalMax-signalMin)/qLevels
 signal = signal / scalingFactor

This gives you a signal ranging from -128 to 128.
You now have the choice of four functions to use:

floor()  %round down to the nearest integer
round()  %round to the nearest integer
ceil()   %round up to the nearest integer
fix()    %round towards zero

I will use round(), then scale the signal back to its original magnitude

signal = round(signal)
signal = signal * scalingFactor

You already seem to have the ability to graph down. I will leave graphing the function up to you.

Edit: To respond to abtj's comment about quant():
I wasn't familiar with quant, but it seems it would work just fine.
The second argument needs to be the value of the least significant bit. This is the same as the scalingFactor as calculated in the code above.

The scalingFactor is simply a way to scale the original signal to the size of the quantization. i.e. scale a signal from -1 to 1 volts to ±8 for a 4 bit quantization. This is to make a function like round() useful, which will only round to integer values.
This is done by taking the range of the original signal (signalMax-signalMin) and dividing it by the number of quantization levels desired.

Please point out any bugs in my code; I don't have access to MatLab to test right now.

Edited: Swapped the * and / when scaling the signal.
Edited: Added Parentheses per A_A's comment.


This paper does a good job of explaining it, based on the energy of a mass-spring system.

R.C. Guido, "Enhancing Teager Energy Operator based on a novel and appealing concept: signal mass", J. Franklin Institute, 2019. https://doi.org/10.1016/j.jfranklin.2018.12.007


Can filter "depth" be adjusted by mixing dry and wet signals?

adding the input to the output of a scaled 2nd-order bandpass IIR will get you the classic peak/cut EQ curve.

adding the input to the output of a scaled 1st-order LPF or HPF will get you a 1st-order low-shelf or high-shelf EQ.

adding the input to the output of a scaled 2nd-order LPF or HPF will get you a particular form of a 2nd-order low or high shelving EQ, but there are many others.  doing it this way will often leave an unintended null or lump in the frequency response.  if you want a symmetric (in log frequency) and perfectly monotonic gain vs. frequency curve, i might suggest referring to the Audio EQ Cookbook.


Validity of differentiation property of Fourier transform

You're right about your intuition, but note that the differentiation property is always valid. This can be seen as follows:

If $x(t)$ is a real-valued function which is constant everywhere except at $t=0$ where it has a discontinuity such that $x(0^+)-x(0^-)=1$, then you have

$$\frac{dx(t)}{dt}=\delta(t)\tag{1}$$

In the frequency domain, this is equivalent to

$$j\omega X(j\omega)=1\tag{2}$$

From $(2)$ we can conclude that

$$X(j\omega)=\frac{1}{j\omega}+c\delta(\omega),\quad c\in\mathbb{R}\tag{3}$$

because $j\omega\delta(\omega)=0$.

Note that the real part of $X(j\omega)$ corresponds to the even part of $x(t)$, and the imaginary part corresponds to the odd part of $x(t)$. For the given $x(t)$, the even part is just its DC value. So if $x_0$ is the DC value of $x(t)$, i.e.,

$$x_0=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}x(t)dt\tag{4}$$

then we must have

$$\mathcal{F}\{x_0\}=2\pi x_0\delta(\omega)=c\delta(\omega)\tag{5}$$

from which we obtain

$$c=2\pi x_0\tag{6}$$

Consequently, the Fourier transform of $x(t)$ is given by

$$X(j\omega)=\frac{1}{j\omega}+2\pi x_0\delta(\omega)\tag{7}$$

For $x(t)=\frac12\textrm{sgn}(t)$ we have $x_0=0$ and consequently,

$$\mathcal{F}\left\{\frac12\textrm{sgn}(t)\right\}=\frac{1}{j\omega}\tag{8}$$

and for $x(t)=u(t)$ we have $x_0=\frac12$ and

$$\mathcal{F}\left\{u(t)\right\}=\frac{1}{j\omega}+\pi\delta(\omega)\tag{9}$$

Also take a look at this and this question, which deal with the same problem in discrete time.


By convention they mean the following:

$x(t)$ = Input in the time domain

$h(t)$ = Filter impulse response (time domain)

$H(w)$ = Filter frequency response

$X(w)$ = Input in the frequency domain

In the equation that you wrote the two "*" symbols mean different things.  In the time domain it is convolution, in the frequency domain it is multiplication.  The equation is just saying that convolution in the time domain (which is how you do time domain filtering) is equivalent to multiplication in the frequency domain.


HI: I haven't read it but this discusses the non-linear case so it may be helpful atleast to get an idea of what the issues are. http://sites.utexas.edu/renato/files/2017/04/recursiveUpdate_v3.pdf


How do I compute the impulse response of a robot?

Impulse response makes sense for linear systems. If your system is nonlinear forget it. If you assume that your system is linear time invariant, you may want to measure the impulse response by exerting a short duration high energy pulse to your robot. It must have sufficiently high energy so that most of the exerted energy is transferred to kinematic energy of the robot to enable motion.

If this is not possible you can apply longer duration pulses to your robot and then measure its response (motion). You must deconvolve then the response and the input to get the impulse response of your robot. You may try this for a number of different shape pulses and if you get impulse responses similar to each other every time, you can be fairly sure that your system is linear.

Your comment on having an infinite impulse response is not in line with correct systems modelling, so forget about it.


Compute Similarity between Fourier Transforms

Most simply, you can compute simple statistical measures on each FFT frame and compare those directly. The signal energy, its distribution, variance, etc across the frequency spectrum are interesting measures. You can also directly compute the correlation coefficient between accelerometer signal pairs (normalized to their length).

It might be useful to look into research on “gait recognition”, the literature points to “low frequency-domain entropy” as feature that is useful in discriminating accel/gyro data that differ in “complexity”

Have you seen this thread on comparing 3d data from accelerometers. It is worth considering if using the raw 3d data, doing a dimensionality reduction (PCA, etc) and then using a simple distance metric may get you further than analyzing raw FFT output. Doing an FFT as you've proposed doing is only one way to simplify the data, but I think the Fourier transform is too coarse & noisy to give you directly comparable results, every time.

As an alternative, it might be easier/faster to use machine learning, i.e. training a "people walking" classifier using any modern supervised learning technique.


You should use the synthesis equation of an impulse train with period $T$ (which is easy to derive):

  $$x(t)=\sum_{k=-\infty}^{\infty}\delta(t-kT)=\sum_{k=-\infty}^{\infty}\frac{1}{T}e^{jk\frac{2\pi}{T} t}\tag{1}$$
  That is: the Fourier coefficients for all terms is a constant ($\frac{1}{T}$).

Now assume that there are two impulses with different amplitudes $a$ and $b$ per period, or $$x(t)=a\delta(t)+b\delta(t+2),\ -4< t \le0$$ In such case we have

$$\begin{align}
x(t)&=\sum_{k=-\infty}^{\infty}a\cdot\delta(t-kT)+b\cdot\delta(t+\frac{T}{2}-kT)\\&=\sum_{k=-\infty}^{\infty}a\cdot\delta(t-4k)+b\cdot\delta(t+2-4k)\\
&=\sum_{k=-\infty}^{\infty}\frac{a}{4}e^{jk\frac{2\pi}{4} t}+\sum_{k=-\infty}^{\infty}\frac{b}{4}e^{jk\frac{2\pi}{4} (t+2)}\\
&=\sum_{k=-\infty}^{\infty}\frac{1}{4}e^{jk\frac{\pi}{2} t}\left(a+be^{jk\pi}\right)\\
&=\sum_{k=-\infty}^{\infty}\frac{1}{4}e^{jk\frac{\pi}{2} t}\left(a+b(-1)^k\right)\\
&=\begin{cases}
\displaystyle\sum_{k=-\infty}^{\infty}\frac{a+b}{4}e^{jk\frac{\pi}{2} t},& k\text{ even}\\[10pt]
\displaystyle\sum_{k=-\infty}^{\infty}\frac{a-b}{4}e^{jk\frac{\pi}{2} t},& k\text{ odd}
\end{cases}
\end{align}$$
Now refering again to $(1)$ and comparing it with the question, we should have
$$c_k=\begin{cases}
\frac{a+b}{4}=1\\[10pt]
\frac{a-b}{4}=2
\end{cases}$$


It appears that your understanding of harmonics is not entirely correct yet. Your example of two sinusoids with frequencies $f_1$ and $f_2$ can actually be explained in terms of harmonics IF
$$f_1=n_1f_0\quad\text{and}\quad f_2=n_2f_0,\qquad n_1,n_2\in\mathbb{Z^+}\tag{1}$$
is satisfied for some positive $f_0$. Note that this means that the ratio $f_1/f_2$ is rational. In that case the resulting sum is periodic. Choosing the smallest possible integer values of $n_1$ and $n_2$ such that $(1)$ is satisfied will allow you to compute the fundamental frequency $f_0$ of the resulting waveform.
In the case of your example we have $f_1=100$ and $f_2=201$ (units don't matter here). So we get $n_1=100$, $n_2=201$, and $f_0=1$ (same units as $f_1$ and $f_2$). So the sum is periodic with fundamental frequency $f_0=1$, even though the Fourier series coefficient of the fundamental is zero! Hence, we have a periodic signal with two harmonics and no fundamental. Nevertheless, the (fundamental) frequency of that signal equals $f_0=1$ (check that!).
With the period $T=1/f_0=1$ we have
$$x(t+T)=\sin\big[2\pi\cdot 100(t+T)\big]+\sin\big[2\pi\cdot 201(t+T)\big]=x(t)\tag{2}$$
where the last equality follows from the fact that we add an integer multiple of $2\pi$ to the argument of $\sin(x)$, which of course has period $2\pi$. Note that this doesn't show that $T$ is the smallest possible period, i.e., that $f_0$ is the largest possible fundamental frequency, but I trust that you can easily verify that.


I explained in "Why is time inversion?"

Red is a reference delta impulse (of height 1), whereas green is a typical response to that impulse, denoted by function h(t). In LTI, output is proportional to the input, that is if we have a delta impulses at the input at time 0, output at time t will be a*h(t). Now, instead of single impulse at the origin, you apply a series of input impulses at various times. What will be the output? Say, there was input impulse of height $a_1$ at $T_1$. Since current time is $t$, impulse occurred $t-T_1$ seconds ago and its contribution to the current output y(t) is $a_1 h(t - T_1)$. There is another contribution from another impulse, occurred at $T_2$. Its contribution is $a_2 h(t-T_2)$. So, $y(t) = a_1 h(t-T_1) + a_2 h(t-T_2)$. You simply add up the contributions beacuse of LTI linearity.

In general, you have $y_t = \sum_0^t {a_i h(t-i)}$. That is a convolution formula.

It also appears when you multiply to polynomials $(a_0 + a_1 z + a_2 z^2 + ...)(b_0 + b_1 z + b_2 z^2 + ...) =\sum_0^\infty {c_n z^n} $ where $c_n = a_i b_{n-i}$. That is why you tend to represent series as z-transforms. In this case you can simply multiply them and have convolution on the background.


How do you do signal averaging on a realtime data?

@Greyfrog. Here are the descriptions of four different kinds of averaging operations:


How to recognise different symbols from their slope ? like (,^) etc

My way of doing it would be:

Extract LBP feature vectors for the reference symbols are store it.
Now extract LBP feature vector of the test symbol.
Compare it with the list of available reference feature vectors and choose the one with minimum difference.

This method can also be extended to other symbols.


Average Energy of modified QAM Constellation

The procedure is always the same. You need to compute the expectation $E\{|A_k|^2\}$, where $A_k$ are the complex symbols of the constellation:

$$E\{|A_k|^2\}=\sum_kP_k|A_k|^2\tag{1}$$

$P_k$ is the probability that the $k^{th}$ symbol is chosen. Usually you can assume that all symbols are equally likely, i.e., $P_k=1/M$, where $M$ is the number of symbols.


How to detect quadrilaterals in a point cloud?

There is one possible solution which is roughly $O(n\cdot (n-1)\cdot log(n))$ where $n$ is the number of points.

Build 2-dimensional k-d tree from all the points
For each pair of unmarked (undiscovered) points

2.1 Determine two areas (based on rectangle criteria) where the other two points can possibly lay

2.2 Find nearest point to center of each of the areas (using the k-d tree)

2.3 If each nearest point is unmarked and lays within its respective area, store new rectangle, mark all four points as discovered

The algorithm can be optimized by keeping undiscovered points in a set, picking new ones from the, progressively reducing its size.

Another speedup would be removing discovered points from the KD-tree so the algorithm does not have to check them in step 2.3.


The ripple you see is due to Gibbs phenomenon and is caused by your filter having a step response. You will have to pass the frequency response through a window function first.


If you are referring to a buffer as in a logic gate, this question may be more appropriate on the electronics stack exchange site but here are some suggestions.

If it is a logic gate, you must be referring to a fixed frequency clock, and shifting the phase for this is quite trivial. Here are some approaches to doing this:

For 180° simply invert with a not gate (inverter).

For 90° and 270° you can double the frequency with a frequency doubler and then invert the 2x clock and then divide by 2 with a flip flop (since it's output changes state on each rising edge). You need to gate the reset on the divide by twos with your original signal to ensure a consistent phase state (either 90° or 270°)

For 90° you can phase-lock a like frequency VCXO using an XOR gate as the phase detector which will have a lock point at 90°, invert this to get 270°.

For 90° you can implement a delay lock loop (similar to above but with a voltage variable delay line and compare one clock output to the delayed one with the xor gate, integrate it's output and use that to control the delay line for a first order DLL which will lock to 90°)

For non precise (temperature sensitive) approximate phase shifts, cascade buffers with RC filters in between.

For variable phase shifts there are plenty of delay lines on the market: https://www.digikey.com/products/en/integrated-circuits-ics/clock-timing-delay-lines/688


length of window and overlap rate in STFT

In general, shorter windows contain less information, and will thus provide less frequency resolution.  But short windows will be more localized in time, and thus a sequence of short windows might allow greater resolution in locating or separating time domain behavior.

Longer windows will contain more information, but this information will be averaged across the entire window duration, thus potentially smearing out any transients or other non-stationary behavior that might change more rapidly than the window length or duration.  But the greater information content might allow a higher resolution (both for separating peaks and locating or estimating isolated stationary frequency peaks).

Overlapped windows will provide some partially redundant information, but perhaps better localized data regarding rate of phase change for some types of signals that are stationary across multiple frames.

Added:  Overlapping is also common when a non-rectangular window is used.  Since the windowing process (combined with finite arithmetic) can be quite lossy at the edges, overlapping helps provide some of the information that was thrown away by the windowing.


Creating Periodic Impulse Train In MATLAB

When you load a signal, it is already sampled. So continuous-time does not make much sense here. Creating an impulse train and convolving it cannot help you either. You should know the sampling rate of the signal you are going to load. Assume it is originally sampled at the rate $f_{s,\text{old}}$. This means the samples are apart in time by $1/f_{s,\text{old}}$ seconds.

Now you want to sample this signal with a sampling rate of $f_{s,\text{new}}$. That is, in the sampled signal you will have $1/f_{s,\text{new}}$ seconds time difference between every two consecutive samples. This problem is called resampling, or sampling rate conversion. There are different techniques to do so which you can learn in the given link or similar sources.

A quick approach to do it in Matlab is as follows:

old_fs = 11200; % or any other value it is
old_time_axis = (0:length(input_signal)-1)/old_fs;
new_fs = 8000; % new sampling rate
[p,q] = rat(new_fs/old_fs); %find the ratio of the sampling rates
resampled_signal = resample(input_signal,p,q);
new_time_axis = (0:length(resampled_signal)-1)/new_fs;

Notice that at sharp signal edges distortion can appear. You can minimize such effects, depending on your signal, by adjusting the resampling filter. Documentation can be found here.


Is applying a linear operation to the input and output of a continuous linear operator OK?

To be a bit picky, the claim does not hold for all linear systems. Here's a counter-example:

Let the system $Tr$ be defined by the input-output relation $$y[n]=Tr\{x[n]\}=nx[n].$$
This system is clearly linear, since

$$Tr\{ax_1[n]+bx_2[n]\}=n(ax_1[n]+bx_2[n])=aTr\{x_1[n]\}+bTr\{x_2[n]\}$$

Now, let's apply the following linear operation $H$ to an input signal: $y[n]=H\{x[n]\}=x[n-1]$. This operation is also linear. Let's see what happens to the claim:

$$Tr\{H\{x[n]\}\}=nx[n-1]\neq H\{Tr\{x[n]\}=(n-1)x[n-1].$$

So, in general the claim does not hold for a linear system.

However: Assuming the claim is stated for linear time-invariant systems (LTI), then it is correct. (The claim was missing the constraint of time-invariance). The proof is simple: Any LTI system is described by its impulse response, and the output of a system is the convolution of the input signal with its impulse response. Given a system $L$ with impulse response $l[n]$, and a system $M$ with impulse response $m[n]$, we get

$$g=L\{f\}=l*f$$
$$M\{g\}=m*g = L\{M\{f\}\}=l*m*f=m*l*f=M\{L\{f\}\}$$

The reason why this is true is the commutativity of the convolution operation (i.e. $x*y=y*x$). Since the commutativity holds for both discrete and continuous time, the claim holds for both discrete and continuous time.


Fast Hartley Transform Implementation in MATLAB

You can resort to the properties of the DHT with the DFT, given in the wikipedia article. For a real-valued function $f$, its DHT in terms of the DFT is given by

$$
\mathcal{H}f = \Re\{\mathcal{F}f\}-\Im\{\mathcal{F}f\}.
$$

So, a function for fast Hartley transform can be given by

function Hf = fht(f)
F = fft(f);
Hf = real(F) - imag(F);
end

Though, you need to take care of the normalization factors, and tdhis is up to you and how you define the scale factors of the DHT.


Channel coder: Does encoder help? Input signal distortion:

BUT: Can encoder help to eliminate effects such as Doppler, offset...?

No, since the encoder doesn't "see" them. The overall system, including en- and decoder can be optimized to deal with such effects.
But honestly, something that corrects a Doppler offset would be called "frequency synchronization", and work on the digitized signal, not on codewords, so it's not much to do with channel coding.
There is techniques that you find both in decoders and for example equalizers, but those are really more shared methods, not "channel coding". Rule of thumb would be for me: if it's something in a receiver that has to do with optimizing the signal prior to mapping it to information, it's digital signal processing, afterwards, it might be channel or source coding.

when we are talking about improving encoder in practice, what does it mean? What should be improved?

This is a  bit to broad, because nothing is ever optimal for all use cases. Example: One man's block code is "optimal" if it's as large as possible to come as close as possible to theoretical limits in code gain. Another man needs to decode something within a very short time, and can't wait for a gigantic code word to be completely received. These two people will have very different ideas what an "optimal" thing is.

Why has an interleaver complicated implementation in practice? How to reduce complexity?

Interleavers are typically among the least complex things in a receiver or transmitter, so not really sure what you're referring to.
NB: there's often things that can be tackled with knowledge that the decoder has; any decision-guided synchronization algorithm is testament to that. In some cases, simple types of ISI might be corrected by a convolutional decoder - but that's really, like, a special case that you can't generally rely on.


Speech segmentation for speaker recognition

I'm going to answer my own question, for those who are looking for the same information as I was. Adaptive Thresholding is the most common way to segment audio, and it gives you the most bang for the buck. While there are more accurate methods out there for segmenting, adaptive thresholding is fast and very accurate, and uses simple features such as zero-crossings, energy, and entropy. I'm going to have to look into this more myself but this should be enough to get started.


What is meant by "Cubic Difference Frequency"?

After skimming through the paper, I can see more clearly now. The measure $D_3$ quantifies the relative strength of the 3rd order intermodulation product. If two sinusoidal signals with frequencies $f_1$ and $f_2$ are input to a non-linear device (such as a microwave amplifier, as referred to in the paper), there will be intermodulation products at the output. Third-order intermodulation products occur at frequencies $2f_1+f_2$, $2f_1-f_2$, $2f_2-f_1$, and $2f_2+f_1$. Problematic are the ones at frequencies $2f_1-f_2$ and $2f_2-f_1$ because they can end up in the frequency band of interest.

The quantity $D_3$ is defined as the ratio of the power of one of these undesired frequency components and the power of one of the desired frequencies (in dB) at the output of the nonlinear device. There are four possible definitions:

$$D_{3,1}=-10\log\left(\frac{P_2(2f_1-f_2)}{P_2(f_1)}\right)\\
D_{3,2}=-10\log\left(\frac{P_2(2f_1-f_2)}{P_2(f_2)}\right)\\
D_{3,3}=-10\log\left(\frac{P_2(2f_2-f_1)}{P_2(f_1)}\right)\\
D_{3,4}=-10\log\left(\frac{P_2(2f_2-f_1)}{P_2(f_2)}\right)$$

where $P_2(f_x)$ is the output power at frequency $f_x$. These four definitions do not necessarily result in the same value for $D_3$. It's just important to specify which of them is used, and to consistently stick to that definition.

To answer your questions: no, $D_3$ is not the power of the 3rd order intermodulation product, but it is the ratio (in dB) of that power to the power at one of the signal frequencies. This quantity is defined based on the assumption that the input signal is the sum of two sinusoidal signals, but it is used to quantify the strength of 3rd order nonlinearities at the output of a nonlinear system regardless of the input signal.

Concerning the power compression factor, you can view the changes $\Delta P_1$ and $\Delta P_2$ as infinitesimal changes, and so the question if $P_1$ and $P_2$ are the values before or after the change becomes irrelevant. At a given frequency, an ideal linear amplifier would satisfy

$$P_2=g\cdot P_1$$

where $g$ is the gain factor at the given frequency. So you have

$$g=\frac{P_2}{P_1}=\frac{dP_2}{dP_1}\Longrightarrow\frac{dP_2}{dP_1}\frac{P_1}{P_2}=1$$

That's why the power compression factor

$$c=1-\frac{dP_2}{dP_1}\frac{P_1}{P_2}$$

is a measure for the nonlinearity of power amplification.


Is this waterfall characteristic of QAM or OFDM (or something other) modulation?

The sharp rectangular shape of the spectrum, with the immensely flat top definitely says "OFDM, with suitable whitening/PAPR reduction"; what constellation is used on the individual carriers is a bit hard to tell.

Can you, with a no-samples-left-out spectrogram (e.g. gr-fosphor sink, or by playing back your recorded file slower and using a suitably high update rate and suitable low number of FFT points in your QT waterfall sink, determine how long the OFDM "burst" is?

Take a bunch of samples from the end, and a bunch of samples from the start (like, one fourth of the burst duration), and calculate a cross-correlation. The maximum correlation should give you the time when what is used as the cyclic prefix at the end of the burst starts; that way, you'll be able to make an educated guess at the symbol number that makes up the CP, and by that the OFDM symbol length/IFFT size. Take frame, and try to manually resample it, so that both CP starts are exactly that (guessed) IFFT size apart. FFT the hell out of the frame; because your resampling had the purpose of time-stretching your signal so that the length of it actually (hopefully) matches the transmitter's sampling, you'd have a timing-synced demodulator.

Look at the spectrum and the constellation plot of the resulting signal; maybe you can already guess the constellation used to modulate the OFDM subcarriers from that.


If the complex frequency bin outputs of the two Goertzel filters are $a$ and $b$, you can calculate the phases as $\alpha = \text{atan2}\left(\text{Im}(a),\ \text{Re}(a)\right)$ and $\beta = \text{atan2}\left(\text{Im}(b),\ \text{Re}(b)\right)$, where $\text{Re}$ returns the real component and $\text{Im}$ returns the imaginary component of the argument. Most programming languages will have the two-argument $\text{atan2}$ function. The phase difference can be calculated as $\alpha - \beta = \text{atan2}\left(\text{Im}\left(a\ \text{conj}(b)\right),\ \text{Re}\left(a\ \text{conj}(b)\right)\right)$, where $\text{conj}$ denotes taking the complex conjugate, saving you one $\text{atan2}$ evaluation and wrapping of the phase, at the cost of one complex multiplication.

There's one more trick available to calculate $a\ \text{conj}(b)$. Multiply the Goertzel output $a$ by the second input signal and lowpass filter the result. In frequency domain it looks like this:

The lowpass filtering isolates the circled spectral peak that corresponds to $a\ \text{conj}(b)$. This saves you one Goertzel filter at the cost of a lowpass filter. If you use this trick you must also take into account the phase shift due to the lowpass filter.

For fast $\text{atan2}$ see Methods of computing fixed point atan2 on FPGA.


Notch filtering line noise outside of frequency band of interest

A general rule of thumb is to filter as little as necessary, because every filter distorts your desired signal, even if just a little bit. Also: a notch filter is advised in a situation where narrow-band interference is present within your signal's band. The notch filter will reject the interference while having minimal effect (ideally) in your signal.

In the problem you describe, my preferred approach would be to go with a low-pass filter only, and increase its order to achieve sufficient rejection of the 60 Hz noise.


Frequency translation of an ideal low pass filter

Multiplying the impulse response by $e^{jTn}$ would give a complex impulse response and a complex-valued filter output. Seems that they were aiming for a real-valued output. They did this by summing two filters, one with the pass band at $T - W \le \omega \le T + W$ and the other with the pass band at $-T - W \le \omega \le -T + W$. If there is overlap between the two pass bands then the frequency response of the composite filter will have a weird shape like:

Seems that they did not intend there to be overlap, but that $T \gt W$.


What do the L, M and S pixel values in the LMS color space stand for exactly?

Colors specified in CIE LMS don't have more or less meaning than colors specified in other CIE based color spaces. CIE LMS coordinates can be readily transformed between CIE XYZ, CIE LAB etc using standard matrices.

However some applications are more suited to some color spaces than others, and LMS has its uses.

In LMS color space, the color matching functions are very close to being the sensitivity functions of the human LMS cones (cone fundamentals). Since each dimension of LMS color space relates to a single cone, it makes it easy to model some of the per-cone gain adjustments that happen in the retina.

cone sensitivity functions:

Basically the retina adjusts the sensitivity of different cones when trying to white balance a scene and some digital imaging white balancing algorithms mimic these operations in LMS space.

This is the paper on CIE LMS color, also known as the '2006 Physiologically Relevant Color Matching Functions' by Andrew Stockman
http://cvrl.ucl.ac.uk/people/Stockman/pubs/2006%20Physiological%20CMFs%20SS.pdf Though it mostly talks about the derivation not the use cases of LMS.

This chapter is also by Stockman and has a lot of discussion on the cone fundamentals http://color.psych.upenn.edu/brainard/papers/Brainard_Stockman_Colorimetry.pdf

--

  To what extent is my understanding correct?

The biggest problem with your understanding is that the values in a LMS triplet do not correspond 1:1 with the # of photons collected in some sensor. They are only relative values. They indicate how the L cone was stimulated relative to the M and S cone, for example.

Every other color space besides LMS also indicates the relative retinal response, just less directly. So LMS is not unique in its capabilities but it is certainly more suitable for chromatic adaptation algorithms.


what exactly does scale mean in scale-space theory?

There's an image processing algorithm called Retinex that uses scales to perform local contrast enhancement. Here's an OK presentation on the Multi-Scale Retinex algorithm for more detail.

The image size doesn't change because you are performing convolution of a Guassian kernel with the image. In practice this is generally accomplished in the frequency domain and the image is padded to avoid the effects of circular convolution. This is standard stuff in image processing where you convolve a 'kernel' against an image. It's a just averaging surround pixels with a 2-D Gaussian function instead of something like a block of ones.

I've seen the choice of sigma be a somewhat adhoc "This looks good to me" selection process. There are ranges published in a few papers that seem to work well for most images. Essentially small values preserve the detail in the image, medium values provide a mix of general detail and global energy, and large values provide more global representations of the image energy. I think something like small is 1-25, medium 25-75, and large is +75. When you convolve and sum all of the scales, with a bit more light math, you get a local contrast enhanced image. The results are pretty good. You can try it out on some images to see the result of each to get a better feel for it.