I have two methods to suggest:
1. Use the erosion (morphological operator) (you can also use other operators like closing/opening). But i think it will destroy the steam data also...
2. A better solution is to detect the lines using hough transform. Google it and you'll find tons of info about it (if your'e using MATLAB, there is a guide how to do it in their website)

Have fun...


What is the equivalent smoothing function to running the same Gaussian 8 times?

From Wikipedia:

  Applying multiple, successive gaussian blurs to an image has the same effect as applying a single, larger gaussian blur, whose radius is the square root of the sum of the squares of the blur radii that were actually applied.

So in your case, that would be equivalent to blurring the image with a single Gaussian with $\hat{\sigma}$:
$$
\hat{\sigma} = \sqrt{8\sigma^2} = 2\sqrt{2} \sigma
$$


Why are MFCCs of two equivalent signals completely different?

You say:

  I stream this file via VLC-Player and Icecast2-Server, receive it in Java (with the Player of Javazoom-Lib)

When you receive the stream in Java, what is the stream format?

You say that it is the "same" because the waveform and/or spectrogram look similar, but the MFCCs will come out different if the stream format (sample rate, bit depth, etc) are different.

Can you confirm by printing out the stream formats

Before sending to Java
Upon receiving from Java?

I recommend using Soxi to print out the stream format (if you cannot otherwise dump it from R, or Java)

If you can eliminate stream format issues, mp3<->wav issues (deal only in wav), if you analysis windows are identical (as jojek says), then your MFCCs should come out identical.


extend parts in an existing audio file to make it sound longer without stretching the file

Sounds for me, like non-linear time stretching can accomplish what you want. In your example, you would apply a stretching function with a slope of one most of the time and for the single note (you want to be played very long) the slope then needs to be strongly different from one.

Overview paper: Jonathan Driedger and Meinard Müller. A Review on Time-Scale Modification of Music Signals, Applied Sciences, 6(2): 57–82, 2016. [DOI]

Code: https://www.audiolabs-erlangen.de/resources/MIR/TSMtoolbox


Well, not really. The wavelet transform (the discrete to be precise) ist just a base transformation of a signal. However, it can be used to compress signals (dimension reduction in some sense), due to the fact that the energy of most natural signals is compactified under the wavelet transform.
The continuous wavelet transform does even the opposite. It blows up your signal and introduces a lot of redundancy.


Relating transfer functions with step responses

In the given example you have 3 types of systems:

the ideal integrator (D) with a step response that grows infinitely (fig. 1); that's the obvious one, as you've found out by yourself.
underdamped systems with complex conjugate poles (A, B, and C): their step response oscillates
critically damped systems with a double real pole (E and F): no oscillation in the step response

So systems A, B, and C must correspond to figures 2, 4, and 5, whereas systems E and F must correspond to figures 3 and 6. You had already arrived at that conclusion.

But we have more cues. Clearly, system C has less damping than systems A and B, which have the same damping. Smaller damping corresponds to more oscillation (i.e. higher overshoot and undershoot). So we can conclude that system C corresponds to figure 4. The difference between systems A and B can be seen by using the Initial Value Theorem:

$$f(0^+)=\lim_{s\rightarrow \infty}sF(s)\tag{1}$$

For system A we get

$$g_A(0^+)=1$$

whereas for system B we have

$$g_B(0^+)=-1$$

These are the values of the impulse responses, i.e., the derivatives of the step response. We see that the step response in figure 5 has a negative derivative at $t=0^+$, from which we can conclude that it corresponds to system B. Then we're left with figure 2 for system A.

The difference between systems E and F is the distance of the double pole from the imaginary axis. The further away from the imaginary axis, the faster will the impulse response decay, and consequently, the faster will the step response rise. So system E must correspond to figure 6, and system F corresponds to figure 3.


Conceptually, I would go for the following steps:

Perform thresholding of the image to create black/white image.
Detect the grid in the image. How you do this depends on the image properties. For example, will the grid lines always be horizontally and vertically or can the grid also be rotated? If it's always horizontal, you can simply add up the pixel values in all rows or in all columns, and you will see the grid lines in the sums. Also, you can detect the grid lines using Hough Transform.
Remove the grid lines from image, e.g. by setting all pixels that belong to the grid to 0.
Use dilation to close the gaps in the curve, where the grid lines were removed.
check, where the curve crosses the grid lines by checking, where the grid image and the dilated image are both 1.


Moving Average Filter's Magnitude Response

A true moving average would have unity gain coefficients but also be divided by the total number of samples (as per the definition of average). This is trivial since it is just a scaling factor, but if the OP divides the output of the moving average by the total number of samples in the filter then the magnitude response will be normalized to 0 on a dB scale, with all other magnitude values given in negative dB.


Calculating distance between two blobs in a binary image in microns

This is relatively straight-forward with OpenCV and numpy.

Step One:
Using OpenCV's Blob Detector you can do something like this:

import cv2

im = cv2.imread("blobs.png", cv2.IMREAD_GRAYSCALE)
params = cv2.SimpleBlobDetector_Params()
params.blobColor = 255
detector = cv2.SimpleBlobDetector_create(params)

keypoints = detector.detect(im)

kpt_image = cv2.drawKeypoints(im, keypoints, numpy.array([]),(0,0,255),cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

cv2.imshow("Keypoints", kpt_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

Step Two:
Write a function to calculate the distance between two keypoints:

import numpy

def distance(kpt1, kpt2):
    #create numpy array with keypoint positions
    arr = numpy.array([kpt1.pt, kpt2.pt])
    #scale array to mm
    arr = arr*40/1280
    #return distance, calculted by pythagoras
    return numpy.sqrt(numpy.sum((arr[0]-arr[1])**2))

Putting it all together:

import cv2
import numpy

def distance(kpt1, kpt2):
    #create numpy array with keypoint positions
    arr = numpy.array([kpt1.pt, kpt2.pt])
    #scale array to mm
    arr = arr*40/1280
    #return distance, calculted by pythagoras
    return numpy.sqrt(numpy.sum((arr[0]-arr[1])**2))

im = cv2.imread("blobs.png", cv2.IMREAD_GRAYSCALE)
params = cv2.SimpleBlobDetector_Params()
params.blobColor = 255
detector = cv2.SimpleBlobDetector_create(params)

keypoints = detector.detect(im)

kpt_image = cv2.drawKeypoints(im, keypoints, numpy.array([]),(0,0,255),cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

for i,keypoint in enumerate(keypoints[1:]):
    print("Distance: {0:6.3f}mm".format(distance(keypoints[0], keypoint)))

cv2.imshow("Keypoints", kpt_image)
cv2.waitKey(0)
cv2.destroyAllWindows()

I assume that when you say 'real time' you may want to look at annotating the image with the distance.  But I'll leave you to work with that yourself.  I've not actually done an error analysis, but you may need to tweak the Params above to increase accuracy.

The output image:


Power Spectral Density as a single number confusion

Question 1: Yes you are correct, the power spectral density is the power distribution per unit frequency so is a continuous function of frequency.
Question 2: The single number as given can be an estimate of total power.
What they gave is completely incorrect starting with the formula as given:
$$P = lim_{T \rightarrow \infty}\frac{1}{T}\int_0^T|x(k)|^2dt$$
This is an attempt to provide a formula for the total power in the signal, as the time average of the energy of the signal (but where is t in the function for x?). I assume $k$ is a typo, but if it refers to an FFT bin, then as $T \rightarrow \infty$ discrete $x(k)$ would become a continuous function of time as $x(t)$ so that the above would properly represent the total power as:
$$P = lim_{T \rightarrow \infty}\frac{1}{T}\int_0^T|x(t)|^2dt$$
Note that the absolute value squared of the FFT bins, $|x(k)|^2$ can be used as an estimate of the power spectral density (see links below for conditions on this), with proper scaling to any units desired as the the relative power in each bin.
For further details on estimating PSD from the DFT, see:
Power spectral density vs. FFT bin magnitude
What is the difference between PSD and squared magnitude of frequency spectrum?


Impact of the sample rate of the emitting signal : communication between HackRF One and USRP NI 2921

The HackRF doesn't support those sampling rates. From reading the firmware over reading libhackrf to reading the gr-osmocom hackrf_sink_c:

The set of supported sample rates is very limited:

osmosdr::meta_range_t hackrf_sink_c::get_sample_rates()
{
  osmosdr::meta_range_t range;

  /* we only add integer rates here because of better phase noise performance.
   * the user is allowed to request arbitrary (fractional) rates within these
   * boundaries. */

  range += osmosdr::range_t( 8e6 );
  range += osmosdr::range_t( 10e6 );
  range += osmosdr::range_t( 12.5e6 );
  range += osmosdr::range_t( 16e6 );
  range += osmosdr::range_t( 20e6 ); /* confirmed to work on fast machines */

  return range;
}

So, use 8 MS/s.

The hackrf is instructed to use a different sampling rate than you think it uses, which obviously "scales" the spectrum.


What should be the dimensions of a Fourier Transformed image?

Here is a third different output. So I have tried to perform an FFT with Matlab, on four versions of your image:

the grayscale version of the colored image,
the red, green, and blue channels.

I do get  the following, with $512\times256$ size, in a $\log$ scale.

If as suggested by @Peter K. you zero-pad below, with $512\times512$ size, one gets:

Looks like  the FFT of the Red component, padded.
If you have the opportunity to do other tests, we can settle this riddle


How to find a certain frequency in a signal that varies in time?

If I understood it right, you intend to figure out how some frequency content of a given signal evolves in time.

Next to short time Fast Fourier Transform, you have the opportunity to make use of a time-frequency analysis. I would suggest you to use either the Hilbert Transform (https://en.wikipedia.org/wiki/Hilbert_transform) or the Teager-Kaiser energy operator (http://folk.uio.no/eivindkv/ek-thesis-2003-05-12-final-2.pdf). Both methods are able to make estimation of istantaneous magnitude and frequency of a signal over time. I assume that your signal is multicomponent, therefore you should first isolate the frequency components you want to analyze by means e.g. of a filter-bank and then apply one of the two mentioned methods.
Otherwise, you could directly go for an Empirical Mode Decomposition (https://en.wikipedia.org/wiki/Hilbert%E2%80%93Huang_transform). Such a method will decompose the signal in terms of frequency components: after that you can again apply Hilbert or Teager-Kaiser methods to check how your target-frequency develops over time.


How can I translate what I have learnt in Transfer Functions to differential equations?

You assign variables for each node that branches off or where signals join (eg the feedback adder). Then you write down the transfer functions between those variables. And as a last step you convert the functions back from Laplace domain into the time domain. Et voilà, you have differential equations.

You actually don't need that many variables. Only input, output and one per integrator and per differentiator is required. But starting off with more variables makes the functions simpler and thus the back transformation. At the end you can always get rid of the excess variables.


(How to ask a Homework Question): Define poles by using proportional controller

Nice work and your conclusion is correct.

Further confirmation is the poles in the closed loop system are the roots of $$1+K G(s)=0$$

Which results in the polynomial with a numerator

$$s^2+(3+K)s+(2-K) $$

In order to have complex roots, using the quadratic formula $\sqrt{b^2-4ac}$ must be negative.

So we want to test if $4(2-K)>(3+K)^2$ for any K

Since $4(2-K)=(3+K)^2$ has roots at approximately -.101, -9.898:

$$K^2+10K+1=0$$

The range of K where the poles can be complex is therefore from -.101 to -9.898.

Showing all poles for the closed loop system as a function of K (root locus) over the range of K=-10.5 to 0 as plotted below confirms that $-4 \pm i$ is not a closed loop pole.

We could mathematically prove this by using the relationship of "Break-Outs" on the root locus:

(See Dr. Cheever's excellent pages at http://lpsa.swarthmore.edu/Root_Locus/Example2/Example2.html)

As Dr. Cheever summarized, break-out or break-in points occur where N(s)D'(s)-N'(s)D(s)=0 where N(s) and D(s) are the numerator and denominator polynomial of the open loop transfer function, and N'(s) and D'(s) are the differentiation.

$N(s)= s-1$

$N'(s) = 1$

$D(s) = s^2+s3+2$

$D'(s) = 2s+3$

Solving this results in roots at approximately -1.4495 and 3.4495 indicating the break out and break in locations. The root locus will break out at -1.4495 when K=-0.101 and follow a trajectory as shown in the diagram above to the break-in at 3.4495 when K=-9.898. Therefore $-4 \pm i$ is not a closed loop pole.

Note: Root locus was plotted using Octave with the rlocus command as part of the control toolbox

pkg load control %if not already loaded
sys=tf([1 -1],[1 3 2]);
rlocus(sys,.1,-10.5,0);


One way to think about this:

Everything that can be represented as a rational function in the Z-domain can also be represented as a linear difference equation in the time domain. As such at can be interpreted (and implemented) as a filter. It may be non-causal or unstable, but it's still a filter.

You can certainly represent things in the z-Domain that are NOT rational functions, however, that would be a system transfer function anymore, because the transfer function is only defined for LTI (Linear Time Invariant) Systems


Effect in Fourier-domain when the signal is "zero-stuffed" in time-domain

Let's be more general and consider the insertion of $L-1$ zeroes between each sample. Mathematically, this can be written as
$$
y[m]=
\begin{cases}
x[n] & \text{if } m = nL\\
0  & \text{otherwise}
\end{cases}.
$$
The DTFT of $y[n]$ is then given by
$$Y(e^{j\Omega}) = \sum_{m=-\infty}^{\infty} y[m]e^{-jm\Omega}.$$
Because $y[m]$ is non-zero only when $m$ is a multiple of $L$, that is $m = nL$, this becomes
$$Y(e^{j\Omega}) = \sum_{n=-\infty}^{\infty} x[n]e^{-jnL\Omega} = X(e^{jL\Omega}).$$
From this, you can see that the DTFT of $y$ is the DTFT of $x$ with a period $L$ times smaller (i.e. $2\pi/L$ instead of $2\pi$). Intuitively, this makes sense: you know that "expanding" a signal in the time-domain is equivalent to "squeezing" its spectrum in the frequency domain (see the time scaling property of the Fourier transform).


Nyquist Maximum data rate formula for PCM

The first document (from mason.gmu.edu) gives a confusing explanation of the formulas. The second document is what's the standard way of referring to the channel capacity

Now the two formulas are:

$C = 2B ~\log_2(M) ~~~~~~~~~~~~,~~~\text{Nyquist}$
$C = B ~\log_2(1 + \text{SNR}) ~~~,~~~\text{Shannon-Hartley}$

Eventhough the first formula, (referred to as Nyquist in the first document), is assumed to yield channel capacity (of a noiseless! channel which is infinite) it's actually giving the necessary minimum data bit-rate to represent an analog message signal of bandwidth $B$ Hz, and quantized to $M$ levels ($\log_2(M)$ bits) using a PCM representation. According to the Nyquist sampling theorem, the minimum allowed sampling rate would be $2B$ Hz, hence there will be $2B$ samples per second each quantized to $M$ levels, yielding a total of $ D = 2B ~ \log_2(M)$ bits per second data rate to represent the analog message signal. Hence it's not channel capacity, it's the minimum necessary bitrate to represent an anlog signal using $M$-bit PCM technique.

The second formula is what we call: maximum channel (reliable) transmission capacity in bits per second rates, or Channel-Capacity in short, over a more realistic noisy channel (such as telephone wires, all sorts of RF channels, and also the fiber optic or satellite links by modeling them with the inherent noise floors). It gives you the maximum achievable rate of information (or data) transmission through a channel of bandwith $B$ Hertz, by using a modulation scheme with the given SNR (ratio of transmitted signal power to the channel noise power).

Note that despite being treated in communications sections, PCM is not a data transmission format but a data representation format. Analog data is converted to digital by samlping and quantizaiton through an ADC of $N$-bit quantization at a sampling rate of $F_s$ samples per second. The resulting sampled-quantized data is practically represented in binary digital system for processing and transmission. In such a case every $N$-bit quantized sample of the message signal is laid out as $N$ consecutive binary digits (bits) to represent a single PCM sample of message data. When it comes to transmission, those bits can be transmitted one by one as is done by using BPSK or in groups of M-bits as is done by 8-QAM (where 8 PCM bits are transmitted at once by a single 8-QAM transmission symbol).


Triangle inequality is your friend: if $y = a_1x_1+  a_2x_2+ a_3x_3 $, then $$|y| \le |a_1||x_1|+  |a_2||x_2|+ |a_3||x_3|$$
Hence if $x_{.}$ is bounded, the sum is bounded as well.


What proportion of a padded FFT should be actual values

Be careful with this in thinking that you would increase your frequency resolution- you won't! Zero padding is very effective in iterpolating more samples between the samples you have, but it does not add any more information about the frequency content of your signal between your original bins.  You do not lose anything either.

The frequency resolution in all cases is $1/T$ where T is the time length of your data.

Please see my answer at this post specifically for more details on zero-padding and how it works: What happens when N increases in N-point DFT


Real time signal processing use cases for eigenvalues of symmetric matrices

OK, you know that complex sinusoids are eigenfunctions (eigenvectors) of linear time-invariant (LTI) systems (which are linear operators). Guess what you see when you look at a spectrum: it's the eigenvalues corresponding to the eigenvector of any given frequency. So, basic signal theory uses eigenvectors to understand the fact that the frequency- and time-domain representation of signals are equivalent.

In MUSIC especially, the eigenvalues are of great importance – they signify the power of that eigenvector of the observed system (i.e. the power of the frequency component in case of that eigenvector spanning the signal subspace). So, MUSIC depends on the eigenvalues.

Basically, all symmetrical MIMO schemes that involve precoding decompose a matrix using the EVD into especially a diagonal matrix with the Eigenvalues as entries – and the question whether, and with what power, a decomposed "virtual" channel is used depends on the magnitude of the respective eigenvalue. So, MIMO, Precoding, and in the context of power assignment waterfilling for these channels depend on eigenvalues.

Now, for things that are "borderline" signal processing:

Graph theory has a relatively rich interaction with linear algebra. The Spectrum of a Graph is basically equivalently used to Spectra of Matrices.

Especially for random graphs, statements on bounds and distribution of eigenvalues of the (often unboundedly large) adjacency matrix are useful to describe degree distributions, and how reliably one can remove random edges from the graph without cutting the graph.

You can imagine the relevance of that to the theory of routing in communication networks; but in the case of random parity-check codes, the distribution of node degrees and cycle lengths is an important thing that can be at lest approximated through this.


Is there a way to evaluate block size when doing image processing or video processing

I assume the block size is related to image quality no matter you are doing the segmentation or compression. Generally speaking, smaller blocks are better as far as correlation (the pixels in the block are more correlated), but if they're too small, you won't be able to eliminate many coefficients from each block without causing severe distortion, especially when there are high variance and more details in the block. Some compression algorithms use different block sizes within the same image so that high-detail regions are divided into smaller blocks, while smooth regions are divided into larger blocks (Chong et al's Compression of medical images through adaptive block-size DCT coding, and The Influence of DCT Block Size on Coding Efficiency: An Adaptive Variable Block Size Transform Coding System with Lagrangian Cost Function Decision Criterion).

Besides, Karunasekera et al's A distortion measure for blocking artifacts in images based on human visual sensitivity provides you a reference to estimate the distortion induced by blocking artifacts. Pramod el al's Unsupervised Segmentation of Medical Images using DCT
Coeﬃcients applied DCT coefficients to the image segmentation.


A transfer function is only defined for linear time variant systems. Since your proposed structure is non-linear, the concept of a transfer function doesn't apply any more. You could make approximations for certain special cases and assumptions, but overall the answer to "what happens to the frequency response" is "undefined".


I just cropped the original gray image:

and converted it to JPEG using irfanview, settings to 15% (my first guess, good feeling isn't?):

You can see see the same yellowish shades on the right, due to heavy chroma quantization. Specifying that the image should remain grayscale gives you this:


How we can find ROC(Region of convergence) given a signal in z domain?

You need to determine the values of the complex variable $z$ for which the series

$$\sum_{n=1}^{\infty}n^2z^{-n}\tag{1}$$

converges. So you need to know a few things about infinite series.

For this case a simple test is the ratio test. You take the ratio of two successive elements of the series and compute the limit:

$$L=\lim_{n\to\infty}\left|\frac{(n+1)^2z^{-(n+1)}}{n^2z^{-n}}\right|=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|z^{-1}\right|=\frac{1}{|z|}\tag{2}$$

The series converges absolutely for $L<1$, i.e., for $1/|z|<1$, or $|z|>1$. So the ROC is the region $|z|>1$.

Note that the ROC would not change if we used any other power of $n$ in $(1)$.

EDIT:

A step for step explanation of Eq. $(2)$:

$$\begin{align}L&=\lim_{n\to\infty}\left|\frac{(n+1)^2z^{-(n+1)}}{n^2z^{-n}}\right|\\&=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|\frac{z^{-(n+1)}}{z^{-n}}\right|\\&=\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\left|z^{-1}\right|\\&=\frac{1}{|z|}\lim_{n\to\infty}\left|\frac{(n+1)^2}{n^2}\right|\\&=\frac{1}{|z|}\lim_{n\to\infty}\left(\frac{n^2+2n+1}{n^2}\right)\\&=\frac{1}{|z|}\lim_{n\to\infty}\left(1+\frac{2}{n}+\frac{1}{n^2}\right)\\&=\frac{1}{|z|}\left(1+\underbrace{\lim_{n\to\infty}\frac{2}{n}}_{0}+\underbrace{\lim_{n\to\infty}\frac{1}{n^2}}_{0}\right)\\&=\frac{1}{|z|}\end{align}$$


The CPU's datasheet will provide information about how many cycles each command requires (e.g. in here look for "cycles", and with something like this you can also work out what else might a command require to determine if additional cycles will be required). Sometimes they might differ (for the same command), depending on the operands.

What you can do is disassemble the processing part, look at the set (i.e. unique) of commands that are used, look up how many cycles they take on average and from that, work-out a best-case / worst-case scenario for your algorithm (taking also into account how many times each command is used, potentially even, in a loop).

If you are already at a low-enough level, so that you are sure about what commands have been used, then your estimate will be closer to the real value.

Hope this helps.


The "symbols per second" (or, more strictly, the baud) unit is a bit weird. It is really a unit that is used for human convenience and not for mathematical rigor. Just like "samples per second", bauds can be interpreted as hertz and still make sense, since they measure the "number of times something happens per second".

So, in your calculation, the units of $B/R_S$ are $\text{Hz}/\text{Hz}$, making it unitless.

I know this is somewhat unsatisfying (at least, it is to me), but this is the best explanation I can give.

(All your unit definitions at the top of the question are correct, BTW).


The difference in implementation between a FLL and PLL is completely dictated by what device is used as the error detector: if the error detector produces an error that is directly proportional to phase (phase detector) that is then passed into the loop filter, then the implementation will be a PLL and we would treat the VCO (Or NCO) as a phase integrator (since that device producing a frequency in proportion to its input, so produces a phase in proportion to the integral of the input). If the error detector is directly proportional to frequency (such as a frequency discriminator) that is then passed into the loop filter, then that device would be an FLL and we would treat the VCO or NCO as a proportional device producing a frequency in direct proportion to its input. Note that we can have a phase detector that we then take the derivative of its output which would then be a frequency discriminator (and this approach can be used for a combined FLL/PLL) so the actual block diagram and functionality needs to be carefully reviewed for each case.
Frequency is simply the derivative of phase with respect to time. So in a phase estimation algorithm you estimate absolute phase (relative to a reference which is typically the phase of a symbol location on the complex IQ plane), while in a frequency estimation algorithm you estimate the rate of change of phase with respect to frequency which is independent of any offset.
Picture a wheel that is spinning at a constant frequency and we use a control loop to stop the spin; if we measure error using a frequency estimation algorithm to then reverse and ultimately stop the spin, it will eventually stop but at an arbitrary offset like we would typically expect from a roulette wheel. If we use a phase estimation algorithm to measure error which is then used to stop the spin (and assuming we can keep up with and overcome the rate of change which implies sufficient loop order) it will stop at the exact place you desire (and you’ll be a winner every time!)


In case the four vectors are row vectors

r2=r(size(c,2)+size(z,2)+1:end-size(z,2))

In case the four vectors are column vectors

r2 = r(:,3)

PS: thanks @ThP


QPSK Data Rate vs Bandwidth at Passband: Why can it only transmit 1 bit per Hz and not 2 bits per Hz?

Let us assume you have a bandwidth $W = 1\,\text{Hz}$ in passband. This means that you have a bandwidth $B = W/2 = 0.5\,\text{Hz}$ in baseband.

In baseband, using sinc pulses, we can transmit at a symbol rate $R_p = 2B = 1\,\text{Bd}$, where $\text{Bd}$ (baud) is the SI unit for "symbols (or pulses) per second".

Using BPSK modulation, we can transmit one bit per symbol, for a bit rate $R_b = 1\,\text{b/s}$.

Now, QPSK is esentially two BPSK signals in parallel. In other words, we have two BPSK baseband signals $s_I(t)$ and $s_Q(t)$, each having a bandwidth of 0.5 herz and transmitting one bit per second.

By using quadrature, we transmit both $s_I(t)$ and $s_Q(t)$ at the same time in passband, using the signal $$ s(t) = s_I(t)\cos(2\pi f_c t) - s_Q(t)\sin(2\pi f_c t). $$

The QPSK signal $s(t)$ has a bandwdith of 1 herz and transmits 2 bits per second (one over the in-phase component and one over the quadrature component).

However, note that the relationship $R_p = 2B$ is valid only for sinc pulses, which are almost never used in practice. If you use other pulses, then the answer will be different. Many introductory textbooks are not explicit about this.

Let us say that you wish to use a pulse $p(t)$. Furthermore, $p(t)$ is a Nyquist pulse. Then, the spectrum of the pulse determines the pulse rate and subsequently the bit rate.

As an example, let us repeat the exercise above for pulses shaped like a half-sine wave. In this case, the pulse rate is $R_p = B/2.5$. Then, if $B=0.5\,\text{Hz}$, $R_p = 0.2\,\text{Bd}$ and, using QPSK, you will only be able to transmit at a bit rate $R_b = 0.4\,\text{b/s}$ per herz in passband.

Note that the $2.5$ figure I used above is my (conservative) personal measurement and reflects what I'm comfortable with. Many introductory textbooks are extremely optimistic about the actual bandwidth of different pulses. In particular, it is common for textbooks to claim that $R_p = B$ for rectangular pulses, which I believe borders on the dishonest (Stallings and Tomasi do this, among others). This is what is behind the claim that "QPSK transmits 1 b/s/Hz" that you (very sensibly) suspect is false.

Note also that the answer you linked to is correct, for the specific pulse shape that is assumed. My answer is more general and allows you to calculate the bit rate for any Nyquist pulse you choose to use.


Fast Vocal Tract Length Normalization(VTLN) implementation?

VTLN is pretty old technology. There is a fast linear VTLN or basis fMLLR which is about the same thing. These days everyone is doing neural networks though, they provide more advanced speaker adaptation, at the same time they could be pretty efficient with quantization and proper weight pruning.

For example you can check DOMAIN AND SPEAKER ADAPTATION FOR CORTANA SPEECH RECOGNITION by Microsoft.

Best adaptation method depend more on the task you are trying to solve, on the resource constrains you have, on the model you are using.


Why do we use white noise to estimate echoes?

The mathematical justification would be that the Cramer Rao Bound on the estimate of time delay is inversely proportional to the time bandwidth product of the pulse, so in cases where there are many distinct echoes, the shorter the pulse is better to resolve each arrival, so one would increase the bandwidth to minimize the bound on time delay.

In the case where the signal is long, as illustrated by the matlab code below for a single delay:

  clear all
  x=randn(1,16*32768);
  x1=[x zeros(1,.5*128)];
  x2=[zeros(1,.5*128) x];
  figure(1)
  pwelch(x1+x2)

The psd has periodic nulls, and you can get an estimate from the null spacing, and the wider the band (white noise), the more nulls that can be resolved.

$$
H(\omega)=1+\alpha e^{\jmath \omega \tau}
$$

For 3 delays it will be a bit more of a bother, but doable.  Actually if you don't know the waveform, like a passive SONAR problem, this is an essential approach.


Is filtering necessary if I'm not downsampling?

Filtering is not necessary as it would only effect those frequency bins outside of the passband.

However, depending on how small your passband is, you could potentially save yourelf a lot of FFT processing by low-pass filtering and then downsampling before your perform your FFT.  For example, if you downsample by 4, the resulting FFT could be 4 times shorter and still have the same frequency resolution.


How about answer your question in an ultimate way? What will happens to if you apply a low pass filter infinite times or practically a large number of times? I believe it is more meaningful for you to understand this one.

I want to first emphasize the fact that apply a filter $F$ to $d$ for $k$ times, is equivalent to apply a filter $F'$ to $d$ for once, if we have
$$F' = \underbrace{F*F*\cdots*F}_{k \,{\rm times}} $$.

The answer to the 'ultimate' version of your question is that: no matter where you start with your lowpass filter, say $F=[f_0,f_1,\cdots,f_n]$, with $\sum_{i=0}^{n} f_i = 1$ and $f_i \geq 0$, you will end up to apply a Gaussian filter $G$ when you apply filter $F$ to some discrete signal $d=[d_0,d_1,\cdots,d_m]$ for $k\rightarrow\infty$ times.

Note: this conclusion does not rely on
1) what is the signal to be processed $d$,
2) which $F$ you use, as along as it is lowpass.

Here is some ipython code, you might play with.

import numpy as np
# generate a random low pass filter
F = np.random.randint(0,1000,(100,1) ).astype('float').ravel();
# ensures sum to 1
F = F / F.sum();
print "sum(F) = ", F.sum();
# create copy to buffer the last version of convolving F for k times
Fk = F.copy();
Fk_list = [];
for k in range( 5 ) :
    Fk = np.convolve( Fk, F, 'full' );
    Fk = Fk / Fk.sum();
    # donot worry about this line, see comments below
    Fk = np.asarray( filter( lambda x : x > 1e-6, Fk ) );
    Fk = Fk / Fk.sum();
    Fk_list.append( Fk );
# plot results
f, axarr = pyplot.subplots( 3, 2 );
# plot data
axarr[0,0].plot( F ),
axarr[0,0].set_title( 'k = 0' );
axarr[0,1].plot( Fk_list[0] ),
axarr[0,1].set_title( 'k = 1' );
axarr[1,0].plot( Fk_list[1] ),
axarr[1,0].set_title( 'k = 2' );
axarr[1,1].plot( Fk_list[2] ),
axarr[1,1].set_title( 'k = 3' );
axarr[2,0].plot( Fk_list[3] ),
axarr[2,0].set_title( 'k = 4' );
axarr[2,1].plot( Fk_list[4] ),
axarr[2,1].set_title( 'k = 5' );
pyplot.show();

the line that with filter keyword is simply to discard very small components, because I want to keep you see beautiful bell shape Gaussian. Otherwise, the results are still Gassuian, but since we might include many things outside of 4 times of standard deviation, you will see something like a peak instead.

Here are the results of the above ipython code. As you can see, the more times of convolve to $F$ self, the resulting lowpass filter is more similar to a Gaussian.

Now let me give you another visual demonstration by using a very different $F$, which only contains two nonzero values. Though this time you cannot see $F$ to be Gaussian like in 5 times of convolution, the long term behaviour namely, $k\rightarrow\infty$ does not change. It will be Gaussian like after a large number of self-convolution.

Finally, you might ask why this thing happens? Because God loves Gaussian. This conclusion is so-called the 'Central Limit Theorem' ( see details in http://en.wikipedia.org/wiki/Central_limit_theorem ). I am not gonna to explain everything for you, because I wish you could learn this on your own. But here is a hint: probability / discrete distribution, arithmetic mean, low-pass filter coefficients, linear operators, and convolution, they are all related. If you have any difficulty to understand something new, try to understand it via something you already knew with close relationships.


How Can I Remove Flickering Background in a Video?

For those classic Video Processing operations there is nothing better than the Plug In's of AviSynth.

Specifically for De Flickering look at:

LMFlicker.
ReduceFlicker.
DeFlicker.

The source code of each is available so you can get inspiration from them.

All of them basically works on smoothing the Luma Data along the temporal dimension.


What would the behavior be of an FFT without the phase (imaginary part)

The real part isn't the amplitude.

The imaginary part is the sine wave part of the FFT.  Thus exactly periodic- in-aperture sine waves (with a phase of zero at the window start) won't show up on the real part, and the FFT real part will see a zero signal.  Only sinusoids with a periodic cosine wave (with a phase of zero at the window start) component will show up at all.

The "beating" you are seeing is your signal alternating between being in and out of phase compared to a cosine basis vector.


$\mathcal Z$-transform, ROC of a system from dividing to others

I agree with your doubts. The $\mathcal{Z}$-transform

$$\hat{G}(z)=1/G(z)=1-3z^{-1}$$

has the ROC $|z|>0$. The ROC of the multiplication of $H(z)$ and $\hat{G}(z)$ must be the intersection of the ROCs of $H(z)$ and $\hat{G}(z)$ (unless there occur pole-zero cancellations, which is not the case here). Since the ROC of $H(z)$ is $|z|>2$, the ROC of $F(z)=H(z)\hat{G}(z)$ must also be $|z|>2$.

Another way to think about it is that since $H(z)$ and $\hat{G}(z)$ are both causal, also $F(z)$ must be causal, i.e., its ROC must be outside a circle in the complex plane.


Yeah, I answered my own question.
So tried the same filter adding zeros to a total of 32 taps.
Plays great now and spectrum is better in the highs.
So this is the filter:

1, -0.15, 0.026, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
(Total taps: 32)


Can Normalized Cut be used to find connected components in binary images?

Lamar, first welcome to StackExchange. Laplacian of a graph has very strong links to the topological properties such as connected components. For simple scenarios, the graph Laplacian, $L=I-D^{\frac{-1}{2}}AD^{\frac{-1}{2}}$ constructed using the degree matrix $D$ and the adjacency matrix $A$, can be processed to result in connected components. First, the number of connected components in the graph is the dimension of the nullspace of the Laplacian and the algebraic multiplicity of the $0$-eigenvalue. Then, for a graph with multiple connected components, $L$ is a block diagonal matrix, where each block is the respective Laplacian matrix for each component, possibly after reordering the vertices.

Simply look at the eigenvalues of the Laplacian matrix and identify clusters (e.g. Spectral clustering). This will give you the individual components if you form the graph with correct adjacency matrix and degrees.

I'm not sure, but if you directly like to use normalized cuts for finding the CC: Simply form the graph ($1$s are adjacent and $0$s are adjacent in themselves). Run the NC and select the segments which do not belong to the background. These should be the connected components.


How to get the values for C1 and C2 in Chan-Vese Model?

The actual formulas are given on the publication, which I linked by editing your question. Besides, there are many publicly available implementations of Chan-Vese algorithm, including, but not limited to:

https://www.mathworks.com/matlabcentral/fileexchange/23445-chan-vese-active-contours-without-edges

https://github.com/kevin-keraudren/chanvese

http://www.ipol.im/pub/art/2012/g-cv/

https://www.mathworks.com/matlabcentral/fileexchange/19567-active-contour-segmentation

http://math.arizona.edu/~rcrandall/code.html

I think you could find all the details in those generously shared source codes.


What does the visual representation of matched filtering of images look like?

There's really nothing special about image processing here: the matched filter is still the hermetian, spatially-inverse filter (=convolution kernel) to the original filter (=convolution kernel).
Since images tend to be real-valued, most convolutional kernels are, too. In which case you just flip the kernel along both axes and are done.
If the original kernel is symmetrical along both axes, it's its own matched filter.
The explanations used for the one-dimensional matched filter are the same as for the 2D case:
Cauchy-Schwarz inequality says that, for any given filter energy, a filter that fulfills the above criteria leads to the highest SNR in the additive uncorrelated noise case.


Questions on Need of filtering before decimation?

If you think for a moment why you filter, it is to remove frequencies that would fold to the remaining bandwidth when decimating.

So if these frequencies that would fold do not exist in the audio then it will have same bandwidth whether you filter it or not. So in the first case, there is really no need to apply LPF first (unless there is noise that would fold and cause problems).

In the second case, it needs to be filtered, and you calculated the cut-off correctly. However, in practice it won't be a brick-wall filter, so depending on needed filter performance and allowable transition bandwith, you might want to select for example cut-off of 4.5 kHz and allow the filter to have about 1kHz transition band to have some reasonable attenuation level at 5.515 kHz.


How are these 2 frequencies aliases of each other?

Note that the cosine has a positive and a negative frequency component. Consequently, in the range $[0,f_s/2]$ you get a component at $-50\textrm{ Hz}+75\textrm{ Hz}=25\textrm{ Hz}$.


This is an expression for the computation of kurtosis.

For a stationary process $x$, moments and cumulants, when they exist, depend on relative time lags, hence $n-1$ variables when the moment or cumulant is of order $n$ (see for instance Basic Signal Processing and Statistics, p. 446 sq.):

$ \mathcal{M} \{ x(k),x(k+\tau_1),x(k+\tau_2),x(k+\tau_3) \} = \mathrm{E} \{ x(k)x(k+\tau_1)x(k+\tau_2)x(k+\tau_3) \}\,. $

A standard notation is $m_{4x}(\tau_1,\tau_2,\tau_3)$, where $x$ denotes the  process, $4$ the $n$-th order, $\tau_1,\tau_2,\tau_3$ the lags, with a triple $0$ for the $0$-lag cumulants. The first $0$-lag cumulants are classicaly named variance, skewness and kurtosis. Your version seems a normalized kurtosis, which could be written, with the above notations:
$$ \kappa_{4x}(0,0,0) = \frac{m_{4x}(0,0,0)}{m_{2x}^2(0)}\,. $$
I am quite surprised thant you only have two zeroes in your parenthesis.

On the  origin side, kurtosis comes for a greek word meaning arched, curved, measuring the tailedness of a distribution.


If we take Z-transform of a signal & find its ROC. How to use this ROC? There are ∞ signals in ROC, suppose we choose any signal what to do with it?

Of particular interest is the unit circle ($z=e^{j\omega}$) when it falls within the region of convergence, since the Z-transform $H(z)$ with $z=e^{j\omega}$ is the Fourier Transform. If we compute the Z-Transform of the system’s impulse response (unit sample response), this will give us the frequency response of the system. We can thus use the ROC to determine if the Fourier Transform even exists.
The Z-transform of a discrete function is only unique if the ROC is also provided. Thus, and importantly, by knowing the ROC we can uniquely determine the original discrete function by using the inverse Z-transform. If the system is known to be causal then the ROC is known to be the outside of a circle which resides on the inner-most pole in the z-plane. This property of the relationship of the ROC allows us to mathematically work with non-causal signals by specifying the ROC as a ring.
The actual locations of the singularities (poles and zeros) provides a lot of additional information about the system not provided with the Fourier Transform alone, and these singularities would be the primary locations of interest within the ROC.
Here is another perhaps more intuitive view that expands on the unit circle of the Z-transform being the Fourier Transform and the concept of convergence with regards to the Z-transform result: Every circle centered on the origin within the ROC is the result of the Fourier Transform of the original function multiplied by a real decaying or growing function given by $|z|^{-n}$ where $n$ is a sample number. Thus if the unit circle is not in the ROC, we can determine from ROC a weighting function that we can use in order to determine the Fourier Transform. We see this by comparing the formulas for the z-transform with the Discrete-Time Fourier Transform (DTFT) for an arbitrary function $x[n]$:
DTFT
$$X(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$
Z-Transform
$$X(z) = \sum_{n=0}^{\infty}x[n]z^{-n}$$
$z$ is a complex variable that has a magnitude and phase given as $|z|e^{j\omega}$ (in general the form $Ae^{j\phi}$ is $A\angle \phi$ for real $A$ and real angle $\phi$), and thus the z-transform is related to the DTFT as:
$$X(z) = \sum_{n=0}^{\infty}x[n]z^{-n} = \sum_{n=-\infty}^{\infty}u(t)x[n]|z|^{-n}e^{-j\omega n}$$
Where $u(t)$ is the unit step function.
Where we see that the z-transform is the Discrete Time Fourier Transform of the causal function $u(t)x[n]$ weighted by the real function $|z|^{-n}$ which is  decaying or growing with $n$ (or constant if on the unit circle).
Following an example provided by the OP in the comments:

Consider a signal having Z-transform X(z) = 3+(2/z) & its ROC is whole
z-plane except 0. Suppose z=1, so X(z) =5

With the explanations above we can see intuitively what is occurring and what this means. First, consider how the z-transform on the unit circle is a correlation of our original function $x[n]$ with the set of complex rotating phasors (individual frequency tones) given by $z^n = e^{j\omega n}$, given that correlation is simply a sum of complex conjugate products. When $z=1$ we are correlating our function $x(t)$ with $z^{n}= 1$ for all $n$, basically a DC signal. So $5$ is the DC response of our system, assuming $x[n]$ was the impulse response (unit sample response) of the system, as then the FT is the frequency response and this would imply that a DC signal at the input to the system would have a gain of $5$. If $x[n]$ wasn't an impulse response but an arbitrary signal, $X(z)=5$ for $z=1$ would imply that the DC content of the signal $x[n]$ (which is the mean value) would be $5$ in proportion to the levels everywhere else for $z$ on the unit circle (basically the FT at $f=0$ or DC, which is a correlation of our signal to all possible frequencies).


Quantized error in scatter plots for a hamming windowed LTE signal

That's the picture you'd get usually if you had a single-carrier system with a slight IQ imbalance. (but that would usually be kind of slanted; but we don't know what your OFDM receiver does, so we can't tell you.)

In a multi-carrier system like OFDM, I'd look into whether I accidentally added inter-carrier-interference (ICI), for example by having an uncorrected fine frequency offset. Or, and that's not that uncommon for simulations, by having a control loop trying to achieve frequency synchronization, but goes wrong. (But that would usually be less "nice and orderly looking". However, we don't know what your OFDM receiver does, so we can't tell you.)

Since this is in simulation: Test. Send a stream of constant symbols, from only one subcarrier. That way, you can see where the effect comes from. Iterate. Debug.

You're the engineer sitting in front of that problem – go and investigate!


Difference between the use of 'bandwidth' in DSP and computer networking

They both define the amount of data you can get through a channel -- the channel bandwidth, together with the SNR, defines the channel capacity $C$, which has the unit of $\frac{\mathrm{bit}}{\mathrm{s}}$ (which, by the way, is simply $1$ binary decision per second), giving you the maximum amount of information you can get through that:

$ C_\mathit{Shannon} = f_\mathrm{bandwidth} \log_2 ( 1+\mathrm{SNR} )$

You see, bandwidth in the signal sense just says "we can see so and so many changes per second", and you can directly translate that to bits you get through that bandwidth, if you know how well your signal gets through.


Slicing and Recognition of book title and author from bookshelf image

Here is the link to a research paper that tries to do the same thing as you wanted. It might help you.using image features
Also a cool video on the youtube


Some common approaches:

Simple post-processing

Apply a mode filter on your resulting sequence of pitches. This will increase temporal inaccuracies though!

Pre-segmentation

Instead of analyzing your signal according to a grid of fixed-size windows, try to detect first the boundaries between notes using an onset detection algorithm. Once you have performed the onset detection, your signal is divided into segments in which you can make the assumption that the pitch is constant. You then apply your pitch detection algorithm to each segment. This works best for instruments with a sharp attack - like guitar or piano (playing monophonic melodies) - but not very well in which the pitch of a held note is modulated without sharp attacks between notes (flute, whistling, singing voice).

Dynamic programming

You define a function $C_i(f)$ which indicates how unlikely it is that frequency $f$ is the fundamental frequency of the $i$-th signal window. You can use, for example, the reciprocal of the autocorrelation function; or sort your pitch hypotheses. For example, if the 4 highest peaks in the autocorrelation function of the 6th signal window are at 220 Hz, 440 Hz, 500 Hz and 250 Hz, $C_6(220) = 0$, $C_6(440) = 1$, $C_6(500) = 2$ and so on. The higher the function, the most unlikely the argument is to be true fundamental frequency.
You define a transition cost $T(f_1, f_2)$ which indicates how unlikely it is that a window with fundamental frequency $f_1$ will be followed by a window with fundamental frequency $f_2$. This transition cost can embed basic musicological knowledge (which notes are more likely to follow other notes) but will mostly serve as temporal smoothing.
You use dynamic programming to find the sequence $f^*$ that minimizes:

$$\sum_i C_i(f^*_i) + T(f^*_{i - 1}, f^*_i)$$

Intuitively, this searches for a sequence of pitches that is both compatible with your observations, but also doesn't change too much. Adjusting the relative weight between $C$ and $T$ will allow you to control the amount of smoothing.

HMM

Manually annotate a large collection of files. Train a HMM using autocorrelation vectors as your feature vectors, and musical notes as your states.

To transcribe a melody, find the most likely sequence of states that explain your autocorrelation vectors. This is done by the Viterbi algorithm, which will take exactly the same form as the dynamic programming method given above. The only difference is that there is this time a rigorous statistical framework for defining and learning the scoring and transition functions.


noise remains after applying median filter

Convert the image to gray scale before applying the median filter.


As of 2019, which discrete nonlinear, time-invariant systems with memory are considered "easy" to model and identify?

It is very strange phenomena that one object is completely dropped out of attention of researchers. It is Urysohn operator. First of all Urysohn is equivalent to multiple parallel Hammersteins and Urysohn followed by static nonlinearity is a model of any deterministic dynamic object, it maps any given input to any provided output. I obtained Ph.D. in modelling dynamic systems as Urysohn operators in 1991. Since that time I showed how to do that to many people, I published articles, I maintain the site, explaining it, and for this almost 20 years noone started using it, there is no single reference on my articles. People write 10,000 lines to identify Hammerstein, I show how to identify Urysohn by 10 lines of code, write here. Look

for (int i = T - 1; i < N; ++i) {
   double predicted = 0.0;
   for (int j = 0; j < T; ++j) {
      predicted += U[(int)((x[i-j] - xmin) / deltaX), j];
   }
   double error = (y[i] - predicted) / T * learning_rate;
   for (int j = 0; j < T; ++j) {
      U[(int)((x[i-j] - xmin) / deltaX), j] += error;
   }
}

That is it. It replace 10,000 line Hammerstein identification. It published, it showed, why noone use it I don't know. Google show my articles on "Urysohn identification" at the top. People don't search.

Explanation to code. x[i], y[i] elements of input output arrays. U[k,m] elements of model matrix. T is one size of matrix U, another one is determined by deltaX, discretization interval, learning_rate is similar to LMS.


Why do we need to normalize a WAV file before calculating the FFT?

It seems that the WAV file used in that example has samples with values between 0 and 255 (likely stored as unsigned chars). The normalization step just changes the samples to floating point values in the range [-1,1). This is, strictly speaking, not necessary to perform the FFT, but it is a good idea. For instance, it removes the large DC component that results from interpreting the unsigned char values as actual signal amplitudes.


What is different between Euclidean distance and the Geodesic distance?

Simple: the Euclidean distance completely ignores the shape when finding a path from the start point to the end point while, for the geodesic distance, the path is constrained to be within the given shape.

That's why the distances at the bottom left of the figure are so different.


Is gradient vector of pole zero carry usful information?

Does gradient vector of pole zero carry useful information?

Yes. The partial derivatives can be used in creating iterative search algorithms for fitting IIR filters to arbitrary targets. Examples of algorithms that use the derivatives are Steepest Descent or Conjugate Gradient.

It's not a trivial process though and there are lot of details to be worked out.

Poles, zeros and transfer function targets are all complex, so you either need to work with complex derivatives are use a suitable real-valued representation (real/imag, magnitude/phase, etc.).
It's important to formulate a suitable error criteria that properly
reflects the requirements and trade-offs for your specific
application,
The search algorithms often get stuck in local minima
Log spaced frequency grids may have to be pre-warped to straighten out the error surface a bit.
Many search algorithm have "heuristic" parameters that are hard to get right


Simplest way to implement FIR filter in C/C++

Have you considered using the Parks McClellan algorithm to generate your FIR filters. The source code for it is available on several sites in Fortran or C.

The original Fortran code is available on Wikipedia.

Here are two sources for C code

Github/scipy
Iowa Hills

You said you were trying to understand filters, so let me explain that the Parks McClellan algorithm is a classic. It was developed in the early 70's, and is still considered to be one of the best FIR design algorithms available.


When you upsample by a factor of two by inserting a zero between each sample you create an alias of your signal centered at the new Nyquist frequency (half the new sample rate).  You get rid of the alias by low-pass filtering.

The specifications of the filter that you need depend on your signal.  Say that your signal has one-sided bandwidth $B$ (by one-sided I mean that we are only talking about the positive frequencies, not the negative ones).  $B$ must have been less than the Nyquist frequency of the old sample rate, $f_{sOld}$, to avoid aliasing.

Once you upsample the alias appears.

Hopefully these pictures will help to gain an intuitive grasp of how to characterize the required low-pass filter.  The passband must be flat from 0 Hz to $B$ Hz, and it must cutoff at $\frac{f_{sNew}}{2} - B = f_{sOld} - B$.  The amount of attenuation that is required in the cutoff region depends on your application.  Usually 50 dB is more than sufficient.


Why does this parallel IIR filter have such a weird form?

The main difference between the form the OP is showing and the referenced solution is that the referenced answer isolates the pole and moves it further away from the unit circle resulting in a more stable design applicable to fixed point solutions. In particular the solution the OP gave has delay accumulation from each element that is within the feedback loop (four multiplier and adder elements) which would limit the maximum clock rate that it could operate in a stable condition for in a realizable implementation.

This is depicted in the diagram below showing the synchronous and combinatorial elements. In order for this diagram to accurately match the first order IIR filter (a single feedback tap), the internal functionality shown must be combinatorial. This will run at 1/4 of the clock rate, but will take significantly longer for the state to settle between clock cycles than the implementation shown in the OP's Wikipedia link copied below this image.

OP's Diagram

Wikipedia Diagram
(note this is only showing two of 4 outputs ($y[4k]$,$y[4k+4]$), in Ben's answer he links a pdf that shows the expanded version of this on p. 35 that shows all four outputs that would be needed after each clock cycle. In either diagraml this partial one or the complete one, the feedback loop is the same, which is the main point: it consists of only 1 adder and 1 multiplier running at 1/4 the clock rate, rather than 4 adders and 4 multipliers in the critical feedback loop).


Why would the audio from .mp4 files be audible on some devices but not others?

When such is the case, it's most probably a broken flash or incompatible browser for those who cannot hear the sound of your videos. However if those two can hear all other internet videos but only yours, then this is really a complex issue. Considering the fact that there are so many different platforms, browsers, devices, drivers and media formats it will be really hard to find out why...

for example I can not even run the videos you have uploaded on my windows 8.1 laptop with firefox browser. It's probably because I'm using a flash blocker to prevent unwanted streaming of flash videos whenever a page is loaded. Also firefox is known to casue erratic behaviour for flash videos through its latest versions...


OpenCV detectin skin colour invariant with illumination changes

In my experience best method for this is converting it to Lab color space.
L represents the light, and a and b are light independent.
OpenCV supports Lab color scale conversion.


How to add pitch-shifted version of signal to itself?

okay, try googling "pitch shifting AMDF" or "pitch shifting TDHS" and see what you hit there.  also check out DSP Dimension and maybe music-dsp and the mailing list.

you need three concurrent processes:

a good pitch detector that gives you a good estimate of the period (to a precision of a fraction of a sample).  this period estimate should be updated as least as often as every 15 or 20 ms.
a real-time resampling and interpolation alg with a pointer that has precision at the fractional sample.  this pointer will have both integer and fractional parts, the integer part tells you where in your delay buffer you are interpolating, the fractional part tells you how you're interpolating (like the interpolation coefficients).  try linear interpolation at first, it's easier and you can improve it later.  now this pointer will advance faster than the write pointer for upshifting and slower for downshifting.  either way, that pointer will run into the edge of the delay buffer (and glitch or click or pop when it wraps around) unless you do something about it.
a splicing and crossfading alg that allows you to splice or seamlessly "jump" from one instantaneous pointer position in the delay buffer to another pointer position.  the difference between those two pointer positions should be one period as determined by your pitch detector.

if i hint more than that, i might be spilling beans (internal trade secrets) from one manufacturer or another.  but this is the main idea.  think of splicing as necessary to time-scale the input (say time-expand if you're upshifting), and then step 2 is playing that time-expanded audio back faster which both returns it to the original tempo, but bumps it up in pitch.


Comparison between Fourier transform, short-time Fourier transform and wavelets

All three transforms are inner product transforms, meaning the output is the inner product of a family of basis functions with a signal. The parametrization and form of the basis functions determine the properties of the transforms.The number of basis functions for a complete picture (i.e. a result that contains enough information to reconstruct the original signal) is  the same for all three cases. That is a result of the Heisenberg uncertainty or a bit more general, conservation of information.

For the standard Fourier transform the basis functions are simply the complex oscillations
$$b_\omega(t) := \exp(i \omega t)$$
where $t$ is the time axis of the signal and $\omega$ is the single frequency parameter that determines the basis function in the family. There is one basis function for every $\omega$. The fourier transform of the signal $s(t)$ is then simply the inner product written as an integral:
$$\mathcal{F}\left\{s(t)\right\}(\omega)=\langle b(\omega,t),s(t)\rangle=\int_{-\infty}^\infty \exp(-i \omega \tau) s(\tau)d\tau$$
The negative sign in the exponential comes from the complex conjugation in the general inner product definition $$\langle a(t),b(t)\rangle=\int_{-\infty}^\infty a(\tau)^*b(\tau) d\tau$$

The short time fourier transform adds a time dimension to the base function parameters by multiplying the infinitely long complex exponential with a window to localize it. The base functions are then something like $$b_{(\omega,t_0)}(t):= w(t-t_0)\exp(i \omega t)$$ where $w(t)$ is the window functions that vanishes outside some interval and $(\omega,t_0)$ are the time-frequency coordinates of the base function in the family. The inner product becomes then
$$\mathcal{S}\left\{s(t)\right\}(\omega,t_0) = \int_{-\infty}^\infty w(\tau)^* \exp(-i\omega \tau) s(\tau)d\tau$$
and gives you time and frequency information. The actual tradeoff between time and frequency is determined by the choice of the window function. The result of this transform can also be regarded as a filterbank with bandpass filters that have the fourier transform as the window $w(t)$ as frequency response, but shifted to the center frequency $\omega$. All the filters therefore have the same bandwidth.

The wavelet transform does away with the constant bandwidth constraint and adapts the window size to the frequency. And that happens in a specific scale invariant way that doesn't even need the complex modulation anymore. You start with a generic base function that is localized and oscillates, and it also has zero mean (i.e. the integral over the complete space is 0). That is your wavelet. Now if you scale the wavelet in time direction the oscillation frequency changes as well, so you can control localization and oscillation with a single parameter that links the two. The family of base functions is therefore something like $$b_{(\sigma,t_0)}(t) = w\left(\frac{t-t_0}{\sigma}\right)$$ where $w$ is the original wavelet and $\sigma$ the scale parameter. The inner product becomes then $$\mathcal{W}\left\{s(t)\right\}(\sigma,t_0)=\int_{-\infty}^\infty w\left(\frac{t-t_0}{\sigma}\right)^* s(t)$$
Again this gives you a 2-dimensionally parametrized result that can also be seen as a filterbank, but this time the bandwidth of the filters is proportional to the center frequency of the bands. The time-frequency plane is therefore partitioned non-uniformly and more suitable to many real-world signals where the typical time constants of processes is naturally scale dependent.


I think you are trying to mix orthogonality for $GF(2)$ and $R^n$. For the sequence $010101000001$, if station 2 had tried to demodulate it, he will have to take groups of 2 bits, multiple (AND) and add(XOR) them, to decode whether bit 0 or 1. Station 2 would have got all 0s because $10$ and $01$ is 'orthogonal' for station 1 and 2. So yes, #3 is indeed a valid CDMA scheme in $GF(2)$. I stress, it is a valid CDMA scheme only in $GF(2)$. It is not related to #1 and #2 (see below why).

Setups #1 and #2 are valid CDMA schemes in $R^n$ (real number vectors) because orthogonality is validated by $x_1^Tx_2$. For #2, you cannot validate orthogonality by dot product of vectors containing 0 and 1. Because for $W_4$, rows 3 and 4 are not orthogonal. So, $0$ has to be represented by $-1$. Even though logically we switch $0$ and $1$ as NOT operator while forming Walsh matrices, orthogonality is not valid in $GF(2)$.

  On the other hand, I cannot make setup #3 work with codes that resemble normal Walsh matrices

$[1,1;0,1]$ would not correspond to CDMA because whether in $GF(2)$ or $R^2$, they are not orthogonal. $[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1]$ is orthogonal in $GF(2)$ as well as in $R^2$ which is same as your multiplexing in time but it is still technically a CDMA scheme (if you transmit each bit at 4x higher rate).


Have a look at these notes.

In particular, this slide:

shows that using a backward difference model sets up the approximation:
$$
\frac{d}{dt} \approx \frac{z - 1}{Tz} = \frac{1 - z^{-1}}{T}
$$

There are many ways to approximate the continuous-time derivative with discrete-time difference operations. That's all that substitution is doing: choosing one way of making the approximation.


The width of the Sinc (both the main lobe and the window) needs to correspond to a low pass filter with a cutoff frequency suitable for the bandlimit of the data and the target sample rate.  Depending on the number of coefficients, a non-rectangular window may be needed.


Are you trying to detect speech vs non-speech, or are there classes of non-speech sounds you are trying to discriminate from? I'm not clear from your question.

I think a decent first approach would be to block your signal into frames and compute the Mel-Frequency Cepstral Coefficients (MFCCs), as well as the delta-MFCCs (differences between adjacent frames' MFCCs) and delta-delta MFCCs (differences between MFCCs in frames that are two frames apart). This isn't the only way to do it, but without more specific knowledge of the problem domain, this is probably a good place to start.

Just googling should give you some good reference on how to compute the MFCCs if you are not already familiar with them. Basically you take the DFT, take the magnitudes, compute the energies inside triangular windows corresponding to human hearing, take the DCT of these coefficients, essentially as a compression step, and then discard the high order coefficients, usually only taking about the first twelve coefficients. I have an explanation of the meaning of the DCT step in this post: How do I interpret the DCT step in the MFCC extraction process?

You could then, say, use these coefficients as features for an SVM.


Thanks to the plot in Olli Niemitalo's answer I got convinced that the formula given in the book has a sign error. The non-linearity used for fuzz or distortion is always some type of smoothed clipping function, which compresses the input signal. So small input amplitudes experience little change whereas high input amplitudes are (more or less) softly clipped. And the figure shown in Olli's answer does exactly the opposite.

So I'm convinced that the correct formula should be

$$f(x)=\frac{x}{|x|}\left(1-e^{-x^2/|x|}\right)=
\text{sgn}(x)\left(1-e^{-|x|}\right)\tag{1}$$

For small values of $x$ we have $f(x)\approx\text{sgn}(x)|x|=x$, and for large (magnitude) values we get $f(x)\approx\text{sgn}(x)$, i.e., clipping.

This is a plot of the corrected non-linearity $f(x)$ (WolframAlpha):

The formula should also be simplified like the right-most expression in $(1)$, because a beginner might be inclined to literally implement the other formula and try to evaluate the terms $x/|x|$ and $x^2/|x|$, which is unnecessarily complex and also gives trouble when $x$ is close to zero. A typical implementation would look like this:

if (x > 0)
   y = 1 - exp(-x);
else
   y = -1 + exp(x);
end


I assume here that your device is not in the feedback chain.

If you can't afford a FFT or filter-bank decomposition (and then detect over successive frames the FFT bins in which the amplitude gets almost exactly multiplied by the same complex number over successive frames), I would suggest looking at these few parameters:

Fit a line to the log of the RMS envelope (computed on 50 or 100ms windows) over the past 2 seconds. Compute the slope, compute the coefficient of correlation. First should be positive, the second should be close to 1 for a feedback build-up, which is exponential.
Compute the standard deviation of the zero-crossing intervals. This is a cheap harmonicity measure, and it should go close to zero as a pure tone settles in and overwhelms the signal.
% of clipped samples. Seeing clipping is a bad omen - but maybe it's too late for you?

You should get by with a few decision rules based on these criteria.

Now, if your device is in the feedback loop (eg: audio DSP box, hearing aid), things are easier because you can "challenge" what you think is a feedback frequency by inducing a delay in the processing chain, or more subtly an all-pass filter (to delay only the suspicious frequency), and if you notice a change in the frequency of the suspect peak, you get a confirmation that it was induced by feedback rather than coming from the input signal.


Correlation Using FFT / IFFT (Convolution in Frequency Domain) in Java

A real (non-complex) time-domain input will always give a FFT that is duplicated (mirrored)

Some implementations will perform a shift to put the DC (frequency=0Hz) term at index 0. Otherwise you will find the DC term halfway through the array.


FFT of square wave: maths and actual numbers do not coincide

The continuous-time Fourier transform of a single rectangular pulse $p(t)$ of duration $[-d,d]$ is :

$$ P(\Omega) = \frac{ 2 \sin(\Omega d) }{\Omega } = 2 d \cdot \text{sinc}( \frac{\Omega d}{\pi} ) \tag{1} $$

where $\Omega$ is the frequency in radians per second.

You can not represent $p(t)$ or $P(\Omega)$ using a sampled-data computer system because $p(t)$ is not bandlimited. However, there is an analogous definition of a rectangular pulse in discrete-time, defined as a sequence $x[n]$ of duration $[-d,d]$, with $d$ being integer, whose discrete-time Fourier transform is :

$$X(e^{j\omega}) = \frac{ \sin( \omega (2d+1) /2 ) }{ \sin( \omega /2 )} \tag{2}  $$

which is called a periodic-sinc or a Dirichlet kernel. Note that this is not a simulation of the continuous case.

Using the computer with N-point DFT/FFT, implementation of discrete-time case will be

$$X[k] = X(e^{j\omega})|_{w_k = \frac{2\pi}{N}k} = \frac{ \sin( \frac{\pi k}{N} (2d+1)) }{ \sin( \frac{\pi k}{N} )} \tag{3}  $$

A Matlab/Octave implementation is as follows

d = 3;       % rectangular pulse in [-d,d]
N = 32;      % DFT/FFT length

k = 0:N-1;   % DFT index

X = sin( (pi/N)*(2*d+1)*k)./sin( (pi/N)*k));
X(1) = 2*d+1;

x = real(ifft(X,N));   % compute inverse DFT on Xk

figure,stem(k,xp);
title('rect-pulse from inverse DFT of periodic-sinc');

Note that the circularity of DFT implies that negative indices of the pulse $x[n]$ are at the right edge.


I'm not sure what your intent is, but the organization of the computation can be done using straightforward overlap-and-add or overlap-and-save methods.

If you're attempting to look at data statistics in the frequency domain, a good starting might be to look at Welch's method and its descendents.


Matlab : Help in getting the correct constellation diagram for QAM modulation

You are getting the correct constellation. You need to remember that you added noise to your symbols in the channel, So, you won't get a perfect constellation as you have in the first figure. Try increasing SNR (EbNo) and you will see better the constellation.

This is what I got:


Energy definition for Autocorrelation lag 0 and lag 1 for complex signals

I had some problems understanding your questions, but from what you wrote I'm assuming you meant $\mathbf x(k) = [x(k), x(k+1)]^T$ so that  $\mathbf x(k-1) = [x(k-1), x(k)]^T$ (using bold faced letters to distinguish scalars from vectors). Is that right?

In this case, expanding your equations, we get the following: $$ R(0) = \begin{bmatrix} r(k,k) & r(k,k+1) \\ r(k+1,k) & r(k+1,k+1) \end{bmatrix}, \quad R(1) = \begin{bmatrix} r(k,k-1) & r(k,k) \\ r(k+1,k-1) & r(k+1,k) \end{bmatrix}.$$

What we can see is that $R(0)$ is Hermitian symmetric (since $r(k,k+1) = r(k+1,k)^*$) and positive definite. Otherwise, unless further assumptions are made, we cannot argue much.

If you assume that your sequence is white sense stationary, we have $r(k,\ell) = r_{\ell-k}$ and we obtain $$ R(0) = \begin{bmatrix} r_0 & r_1 \\ r_1^* & r_0 \end{bmatrix}, \quad R(1) = \begin{bmatrix} r_1^* & r_0 \\ r_2^* & r_1^* \end{bmatrix}.$$ Moreover $r_0$ represents the signal's variance ("energy") and you have $|r_1|\leq r_0$ and $|r_2|\leq r_0$ from Schwartz' inequality. In particular, you can say that $r_1 = r_0 \rho_1$ where $\rho_1$ is the normalized (Pearson) correlation coefficient between adjacent lags of your sequence. Its magnitude is in $[0,1]$ (0 for an uncorrelated sequence, 1 for a constant [or alternating] sequence) and its phase can be anything. This answers why $r_1$ is complex and how to interpret real and imaginary part. It's only close to $r_0$ if the temporal correlation is high. For uncorrelated sequences its zero and hence nowhere close to $r_0$.

Now, $R(0)$ will have eigenvalues $r_0 \pm | r_1|$, so the max eigenvalue is between $r_0$ (for an uncorrelated sequence) and $2r_0$ (for a constant sequence). The eigenvalues of $R(1)$ will be $\left(r_1 \pm \sqrt{r_0 r_2}\right)^*$. From this, one can see that the dominant eigenvalue of $R(0)$ will be greater than or equal to the dominant eigenvalue of $R(1)$ with equality only for a constant sequence.

Finally, the eigenvalues of $R(0)+R(1)$ can be shown to be $ r_0 + r_1^* \pm \sqrt{(r_1^* + r_2^*) (r_1 + r_0)}$, so it's not just the dominant eigenvalue of $R(0)$ plus something.


Classifying 2D shapes by the smoothness of their boundaries

Fourier shape descriptors are quite easy to use and can do well to differentiate smooth objects from jagged ones. Imagine a polar coordinate system with the origin at the centroid of the 2D object. Store a vector of $r$ values as $\theta$ varies in $[0, 2\pi)$ where $r$ is the distance of the boundary from the centroid at each fixed angle $\theta$. Next, take the Fourier transform (FT) of this vector of $r$ values. If the boundary of the 2D object is jagged, the FT will have lot of non-zeros even at high frequencies. On the other hand, imagine doing this with a perfectly smooth circular object---you will end up with a vector of constant $r$ values whose FT is simply a non-zero DC term.

The DFT vector can be used as a feature vector.

See slides 16-20 in this pdf:
http://homepages.cae.wisc.edu/~ece533/project/f06/karrels_ppt.pdf


What is the difference between a'mixer' and a ' multiplier' used in modulation process?

For most purposes a mixer is a multiplier.

The issue regarding sum / difference frequencies is that when you multiply two signals of different frequencies, $\omega_1$ and $\omega_2$ say, then you get one component at the sum of the two frequencies and one component at the difference of the two frequencies:
$$
\cos(\omega_1 t) \cdot \cos(\omega_2 t) = \frac{1}{2} \left\{ \cos([\omega_1+\omega_2]t) + \cos([\omega_1 - \omega_2]t) \right \}$$

Response to comment

A frequency multiplier is a completely different kettle of fish. These are often highly nonlinear systems that take a signal $\cos(\omega_1 t)$ and produce $$\sum_{k=1}^{N} \cos(k \omega_1 t + \phi_k)$$

For example, a frequency tripler can be had by driving a transformer into saturation (so that it outputs something like a square wave) and low-pass filtering the result.


Is magnitude information enough to reconstruct an audio signal

Here is a MATLAB experiment for phase:

[x Fs Nbit]= wavread('test.wav',1024*100); % Get a piece of song ;)
x = x(:,1); % Get a single channel
Y = abs(fft(x)); % Take FFT of the entire piece at once!
y = real(ifft(Y)); % Throw away all-phase information and reconstruct y

y2 = zeros(1,length(x)); % NOW: we will process in BLOCKS!
for i=1:1024:100*1024 % get a block, take its fft, reconstruct a block from fft magnitude
 y2(i:i+1023) = real(ifft(abs(fft(x(i:i+1023)))));
end

figure,plot(x); % piece of a song
figure,plot(y); % obtained from IFFT of MAGNITUDE only of complete piece
figure,plot(y2);% obtained from IFFT of MAGNITUDE only of 1024 sample BLOCKS

sound(x,Fs,Nbit); % the original
sound(y,Fs,Nbit); % from single piece, this is garbage!
sound(y2,Fs,Nbit); % block based, robotic but inteligible!


Fequency Domain Oversampling - decimating when going back to time domain

Your time domain waveform after 6 will be zero-padded as in step 3 with the additional "smearing" due to your equalization to remove the channel effects that will go into the zero padded area. The entire waveform in time will be delayed by the group delay of your frequency domain equalization, which may not necessarily match the group delay of the channel distortion.

If you applied no equalization (by multiplying in frequency by your frequency domain equalization); the IFFT would recover your time domain waveform in its original form with the noise and channel distortion with the additional zeros at the end that were introduced by the zero-padding when taking the FFT. (Along with boundary conditions at the transition such as Gibbs phenomena due to the effective rectangular window applied). Multiplying by your equalizer in frequency is identical to convolving in time the original time domain waveform with the Fourier transform of your equalizer fucntion, and therefore will delay this waveform based on the group delay of this filter (where group delay is $\frac{d\phi}{d\omega}$). Equalizers are typically not linear phase (unless your channel distortion was which is unlikely) so the group delay is not constant versus frequency; but the point is the resultant delay of the recovered original waveform will be completely dependent on the net delay of the channel and equalizer. With knowledge of the channel and specific equalizer function, the delay can be determined.


Does a RHP zero imposes some limitation in the time domain?

If a RHP zero imposes limitation on the bandwidth in the frequency domain, it also imposes limitation in the time domain. This is due to time-frequency duality.

For example, if you want to have a really fast step response, a slow RHP zero will prevent you from reaching that goal.

I really recommand Murray's book for more information http://www.cds.caltech.edu/~murray/amwiki/index.php/Second_Edition


What is the magnitude of the impulse function at t=0?

The "impulse function" is not actually a function but a distribution, and its value at $t=0$ is not defined.

The impulse only really makes sense when multiplied by a signal inside an integral, where it is used to "select" the value of the signal at specific time instants.

Personally, I try to avoid impulses whenever possible. Look for Lapidoth's book on digital communications (available for free on his website) for a complete treatment of Fourier analysis (as relevant for communications) without a single delta. Not many people share my squimishness with deltas, though, and if you're taking a signals course they are pretty much unavoidable.

Note also that, once you go to the discrete time domain, everything becomes nicer: instead of the Dirac delta, you use the Kronecker delta, which is a regular function and has a value of 1 at $n=0$.


At what amplitude does spectral leakage start to be noticeable in the phase spectrum?

The issue is that you are not taking a large enough piece of your original x to get an accurate frequency domain representation.

Try setting

x = exp(-linspace(-20,20,N).^2);

If I do that, I get the bottom two plots below (the middle plots are for -10,10).  The top plot is a replication of your values.

In the three pairs below, the left plot is the semilogy plot of x and the right plot is the semilogy plot of its FFT f.

Derivation of the FFT of a windowed Gaussian

To answer your question in the comments, we need to look at how a windowed version of:
$$
x(t) = e^{-t^2/2}
$$
can be analyzed in the frequency domain.

The windowed version is:
$$
x_w(t) = w(t)e^{-t^2/2}
$$
where $w(t) = \mbox{rect(t/T)}$ if we use a rectangular window of width T (span from $-T/2$ to $+T/2$).

The Fourier transform will be:
$$
X_w(\omega) = \int_{-T/2}^{T/2} e^{-t^2/2}  e^{j\omega t} dt\\
= 2 \int_{0}^{T/2} e^{-t^2/2} \cos(\omega t) dt\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + j \mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} )
 \right] \right |^{T/2}_0\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}} )
 \right] \right |^{T/2}_0
$$
where this integral is used and the identity
$$
\mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} ) = \mbox{erf}(j \frac{ \omega - j t}{\sqrt{2}} ) / j = -j \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}})
$$

Unfortunately, I'm not sure where that gets us. I'll leave it for now and see if I can doodle more later.

scilab code below (Sorry, I tend to use scilab rather than Octave).

// Q26236;
N = 101;
x1 = exp(-linspace(-5,5,N).^2);
f1 = fftshift(fft(x1));
x2 = exp(-linspace(-10,10,N).^2);
f2 = fftshift(fft(x2));
x3 = exp(-linspace(-20,20,N).^2);
f3 = fftshift(fft(x3));

clf
subplot(321)
plot2d("nl",abs(x1))
title('x1')
subplot(323)
plot2d("nl",abs(x2))
title('x2')
subplot(325)
plot2d("nl",abs(x3))
title('x3')

subplot(322)
plot2d("nl",abs(f1))
subplot(324)
plot2d("nl",abs(f2))
subplot(326)
plot2d("nl",abs(f3))


Filter order vs number of taps vs number of coefficients

OK, I'll try to answer your questions:

Q1: the number of taps is not equal the to the filter order. In your example the filter length is 5, i.e. the filter extends over 5 input samples [$x(n), x(n-1), x(n-2), x(n-3), x(n-4)$]. The number of taps is the same as the filter length. In your case you have one tap equal to zero (the coefficient for $x(n-1)$), so you happen to have 4 non-zero taps. Still, the filter length is 5. The order of an FIR filter is filter length minus 1, i.e. the filter order in your example is 4.

Q2: the $n$ in the Matlab function fir1() is the filter order, i.e. you get a vector with $n+1$ elements as a result (so $n+1$ is your filter length = number of taps).

Q3: the filter order is again 4. You can see it from the maximum delay needed to implement your filter. It is indeed a recursive IIR filter. If by number of taps you mean the number of filter coefficients, then for an $n^{th}$ order IIR filter you generally have $2(n+1)$ coefficients, even though in your example several of them are zero.

Q4: this is a slightly tricky one. Let's start with the simple case: a non-recursive filter always has a finite impulse response, i.e. it is a FIR filter. Usually a recursive filter has an infinite impulse response, i.e. it is an IIR filter, but there are degenerate cases where a finite impulse response is implemented using a recursive structure. But the latter case is the exception.


Sorry, your understanding is mostly wrong.

Sound is a pressure & velocity wave in a fluid medium, i.e. air.
A human speaker will generate a 3 dimensional sound field around them. The sound pressure and particle velocity will be different at each location.
A microphone measures the "sound" at one location. Depending on the type of microphone it either picks up the pressure, the particle velocity or a mixture of both. Let's stay with sound pressure to keep it simple.
The microphone converts the local sound pressure into voltage. That is the signal that you see on the scope in the video.
Human speech consists of different types of sound. There are three major types
Voiced sounds. These are sounds that have an audible pitch. This includes vowels and nasals, such as "nnnn", "mmm", etc. Everything you can "sing" or "hum" to make a melody. They can be as long as you want.
Fricatives: these don't have any pitch and are more like noise. Example "ssss", "ffff", "shshsh" etc. They can be as long as you want.
Plosives: these are sounds that are short and just a single spike. "p", "k", "t" etc.

Your whole notion of fundamental + harmonics only applies to "voiced sounds", i.e. things that have a pitch. These can indeed be modelled by a fundamental and harmonics and, for example, the location of the harmonics determines whether a vowel is an "aaaa" or an "eee".
What you see on the scope is the sum of the fundamentals and the harmonics (and everything else: noise, reverb, air condition, ...). Technically it's not a wave, it's a signal.
Physically speaking: a wave is something that propagates over both over space and time.
The microphone has already resolved space, since it samples the sound wave only at one location and the dependency on space is lost. While it's still in air, it's a wave, but once it got picked up by a microphone it's just a signal. This is partially just "semantics", many people would still call it a "waveform" since it's a signal derived from a wave.
EDIT: answering the questions from the comments

(1) You speak of particle velocity; my understanding is that sound travels at 343 m/s in air. Are you saying this is not a constant, and that different waves/signals produce sound traveling at different speeds (lets say air is the fluid medium)?

No. Particle velocity and wave propagation speed are different things. Particle velocity is simply the speed with which  air molecules wiggle back and forth. The higher the intensity, the more they wiggle and the higher the velocity. Wave propagation speed is constant and doesn't change with intensity. It's a function of the compressibility and density of the media.

(2) My understanding is that pitch is the frequency of a wave/signal;

No. Pitch is perceived. If a human can tell whether something is "higher" or "lower" it has pitch. For signals like sine waves and periodic signals, pitch correlates highly with the period of the signal. Other signals have no discernable pitch.

wouldn't fricatives and plosives (which are sounds) also have amplitudes and frequencies, hence, pitches?

No. They don't have pitches since they are not sufficiently periodic.
(3) What would be the correct definition for "fundamental" and "harmonics" here?
Any periodic signal can be expressed as a sum of the sine waves that are integer multiples of the fundamental frequency. That's the whole concept behind the Fourier series.

(4) when you say "What you see on the scope is the sum of the fundamentals and the harmonics..." what would you call this signal? Is there a special name for it?

Not really. It's just a signal.

(5) whatever the name of this signal is, I'm intrigued by your use of the word "sum" here; that its the "sum" of all these other things. Again, I would have expected each distinct "sound" to form its own signal on screen. Any chance you could elaborate as to how all of these things gets "summed" into a single signal?

I think you have it backwards. The sound signal just is what it is: it's air molecules wiggling back and forth and the microphone detects whether the pressure goes up and down. The concepts of fundamentals and harmonics (etc) can be used to analyze the signal or break it down into different components, but it's not the way it's generated or propagates. It's more of a math tool to make predictions and analysis easier.


Basic sampling excercise - unsure of rationale informing np.arange parameters for time array

The Hz values are used to get your parameters for defining your signal.

I don't think there is any reason that they went from -1 to 1 instead of, say, 0 to 2.  Notice that the end point is not included in the range.

The spacing in your sampling is 1/100 of a second.  Yes, the array is giving you values on a time scale.

To generated the signal your are going to need:

f = 2
signal = np.sin( f * 2*np.pi * t )

The $2\pi$ comes as a unit conversion of radians per cycle.  Therefore the argument, in units is:

$$ \frac{cycles}{second} \cdot \frac{radians}{cycle} \cdot seconds = radians $$

Hope this helps.

Ced

================================

Followup:

Here is a sample program.

import numpy as np
import matplotlib.pyplot as plt

f = 2.0
fs = 100.0

t = np.arange( -1.0, 1.0+1.0/fs, 1.0/fs )
signal = np.sin( f * 2*np.pi * t )

plt.plot(t, signal, c='g')
plt.show()

Notice that I added a bunch of ".0"s to force Python into using floats.  Without it, the program passed 1/fs as the integer 0 to arange.


adequate sampling frequency in measuring acoustic emission

How are you determining your maximum frequency?

If you sample at $f_s$ the frequency resolution (the frequency difference between frequency components) of your FFT is $f_s/N$, where $N$ is the number of samples collected. And the maximum resolvable frequency is the Nyquist frequency corresponding to the sample rate. For more information see here


The calculations the OP is using is correct to predict the maximum carrier and symbol timing offset due to reference clock offsets in the transmitter and receiver.

For sampling clock offsets (both in the DAC in the transmitter and the ADC in the receiver), any offset of the sampling clock also introduces a carrier frequency offset that will be in proportion to the carrier that the analog signal is at when sampled according to $f_a/f_s$, where $f_a$ is the analog carrier frequency (of offset from $f=0$) and $f_s$ is the sampling rate. Thus a ppm change in the sampling clock will result in the same ppm change of the carrier of the sampled signal, given the actual carrier the signal is at when sampled. For example if an analog signal was translated down to a 25 MHz IF (meaning a 25 MHz carrier), and was sampled with 100 MHz sampling clock, then $f_s=100e6$ and $f_a=25e6$ in this case, and the fractional frequency of the digital IF would be $f_a/f_s = 0.25$ samples/cycle.  If the sampling clock was actually 1% higher, or 101 MHz, the fractional frequency would 1% lower given $f_a/f_s \approx 0.2475$. Be very careful to reference the actual analog carrier and not the digital IF which can be different in undersampling applications. For example, if the signal was at a 125 MHz analog IF and sampled with a 100 MHz clock to produce a 25 MHz digital IF, the fractional frequency with no clock error would be $(f_a-f_s)/f_s = 0.25$ but if the sampling clock was 1% higher, or 101 MHz, the fractional frequency would be $(f_a-f_s)/f_s \approx 0.2376$ which is 5% lower! This works out to to be the clock error in ppm multiplied by the ratio of the analog IF frequency to it's equivalent frequency in the first Nyquist zone, in this case $125/25 = 5$.

Additional carrier frequency offset can be introduced through channel effects such as Doppler due to the motion of the transmitter and receiver that also must be taken into consideration for determining the maximum carrier and symbol timing offsets.

Once at baseband, a symbol phase offset would typically be referring to the carrier phase offset itself (rotation of the symbol) since an added phase at the carrier directly results in an added phase of the symbol. Even though it would be a referred to as a carrier phase offset, it may not have originated at the original carrier of the signal but can be introduced by other factors in the processing of the signal. If this carrier phase is changing with time, then this is a carrier frequency offset given frequency is the derivative of phase with time. Symbol timing offset at baseband in contrast would refer to the offset between the ideal sampling location for minimum error rate versus the actual sampling location. (literally a sampling time offset, or delay), and there is no rotation involved when the sampling location on the baseband signal is delayed in time. Similarly if a waveform was sampled directly at baseband, only timing offset would be introduced if there were phase or frequency offsets in the sampling clock since the IF frequency at baseband = 0 (and a frequency offset in the sampling clock would generate a varying time offset on the baseband constellation but not a rotation). Any small residual carrier offsets would be affected just like the IF frequencies above if they were then sampled with a clock that has offset errors, but this will typically be negligible given the effective IF frequency is so small.


discrete wavelet transform matrix for vectorized image

For only one level on both direction, the classical DWT schemes are applied separately (rows then columns, or colums then rows).

For more levels, there exist two schemes:

one interleaving decompositions: row, col, row, col, etc. (also called non-separable/square/non-standard/isotropic/Mallat),
one on all rows, then all colums (also called separable/rectangle/standard/anisotropic/tensor).

The second one only deserves the "fully separable" mention.

References to related questions:

2D DWT Image Issue
What is the correct order of operations for a 2D Haar wavelet decomposition?
2D DWT computation order


Prevent negative image with BPSK modulation in a complex signal

An actual BPSK signal can be expressed as
$$\begin{align}
s(t) &= \operatorname{Re}\left[\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)e^{j(2\pi f_c t + \theta)}\right]\\
&= \left[\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)\right]\cos(2\pi f_c t + \theta)
\end{align}$$
where

$b_n \in {0,1}$ is the $n$-th data bit (note that $(-1)^{b_n} = \pm 1$),
$p(t)$ is the baseband pulse shape. Typically $p(t)$ is nonzero only for
$t \in \left[-\frac{T}{2},\frac{T}{2}\right)$, e.g.
$\operatorname{rect}\left(\frac{t}{T}\right)$.
$f_c$ is the carrier frequency which is sometimes chosen to be a harmonic of
the bit rate $T^{-1}$
$\theta$ is the carrier phase which might be arbitrary (depending
on initial conditions in the local oscillator when it is turned on) or carefully controlled to have
value $0$ or $\pm \frac{\pi}{2}$ etc. by MATLABi people.

Notice that the term in square brackets in the second displayed equation
above is a series of nonoverlapping pulses $\pm p(t-nT)$ and that
the value of $b_n$ affects the BPSK signal phase only during
$\left[nT-\frac{T}{2},nT+\frac{T}{2}\right)$.

With that as a preliminary setting of the stage, note that the negative
frequency image comes about from taking the real part in the first
displayed equation above.  We have
$$\operatorname{Re}\left[e^{j(2\pi f_c t + \theta)}\right]
= \frac{e^{j(2\pi f_c t + \theta)} + e^{-j(2\pi f_c t + \theta)}}{2}
= \frac{e^{j(2\pi f_c t + \theta)} + e^{j(2\pi (-f_c) t - \theta)}}{2}.$$
Furthermore,

  no filtering can remove the negative frequency image.

You might say "Well, let's not take the real part at all
since that's the source of the trouble!" and say
that
$$\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)e^{j(2\pi f_c t + \theta)}$$ is a BPSK signal
containing positive frequencies only. Unfortunately, that is
a complex-valued signal that cannot be transmitted on a single
wire or radiated from a single antenna.  So you need to live with
the negative frequency image: it is there to make your complex
signal more real.


Do FFT-based filtering methods add intrinsic latency to a real time algorithm?

You are correct. FFT based processing adds inherent latency to your system. However there are ways to tweak this.

Let's assume you have an FIR filter of length "N". This can be implement FFT-based using the standard overlap add or overlap save method, where the FFT length would be 2*N. Overall system latency will also be roughly 2*N: you need to accumulate a frame of N samples and then while you are accumulating the next input frame, do the math on the current frame. By the time the second frame has been accumulated (time offset 2*N), the first frame is ready to go out. If you have a fast CPU you can speed this up somewhat by using different alignment for input and output frames but that's typically more bother than it's worth.

You can also break down the filter into K smaller blocks of length M, i.e. N = K*M. The FFT only needs to be done over M input samples and the delay and accumulation over the multiple filter sections is done in the frequency domain. Bill Gardner described a few flavors of that here: http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF
It's often referred to as a "Block Convolver".

This allows to basically trade-off latency against efficiency and it provides a continuum between the direct FIR and the full-sized overlap add method. A nice side-effect of the Block Convolver is that the total filter length doesn't have to be a power of 2. For example you can implement a filter of 768 tabs as 6 blocks of 128 without any significant loss of efficiency.


Intuition behind image derivative using Fourier Transform for edges detection

In general, the time derivative property of the Fourier Transform is given as

$$\mathscr{F}[\frac{d}{dt}x(t)] = j\omega X(j\omega) $$

Notice that we can simply multiply by the frequency index in the Fourier Transform result.

For the 2D FT result:

$$\mathscr{F}[f(x,y)]= F(u,v)$$

Using the same property results in:

$$\mathscr{F}[\frac{d}{dx}f(x,y)]= uF(u,v)$$

$$\mathscr{F}[\frac{d}{dy}f(x,y)]= vF(u,v)$$

The reason this detects edges is clear if you observe how the differentiation of a 1D step function results in an impulse. Your edges approximate steps. The derivative is the slope, so where the slope is high (edges) the derivative is large.


what does frequencies in non periodic signals mean?

The physical siginificance of the Fourier transform for nonperiodic (finite energy) signals is not something clear from any point of view.
The limiting case of a Fourier transform of a periodic signal, as the period goes to infinity, consists of line components (dirac impulses at a countable set of discrete frequencies) which becomes more and more crowded, eventually approaching a continuum of frequencies, at which are no more dirac impulses.
This infinitesimal frequency component at any $\omega$, has no physical significance, unlike the dirac impulse at $\omega$, which has a direct implication of a sinusoidal signal component, extending from minus infinity to infinity (in time for example).
Non-periodic signals are also known as transients, or non-stationary signals, and for their analysis, Wavelet transform family can provide a better insight than the Fourier family.
Nevertheless, the Fourier family of transforms can still be used to analyse periodic as well as nonperiodic signals with great success.