Note that the condition

$$\int_{-\infty}^{\infty}|f(t)|^2dt<\infty\tag{1}$$

(i.e., that the signal $f(t)$ has finite energy) is very restrictive when we try to model signals, even though obviously any actually occurring signal must have finite energy. Modeling signals as random processes means that we ignore condition $(1)$. Models are always unrealistic to a certain degree, but many signals can be described very well by random processes even though the signals have finite energy and their models do not. This aspect of the model is often irrelevant.

One example which may serve to understand this fact a bit better is the frequently used model of a (wide-sense) stationary process. Certain statistical properties of such a process do not change over time, and, consequently, realizations of such a process will generally not decay as $t\rightarrow\pm\infty$, and $(1)$ will generally not be satisfied, even though we are only interested in the properties of that process during a certain finite time window. However, power and the power spectrum can be defined for such processes, and most practically useful processes have finite power (or can easily be made to have finite power).


how could window function really suppress those frequencies that we don't want?

You are right that the effect of windowing is the same for the design of FIR filters and for spectral analysis. You will always have a compromise between main lobe width and relative side lobe level (relative to the main lobe). The main lobe width determines the frequency resolution, and the relative side lobe level determines the amount of spectral leakage. A rectangular window has a narrow main lobe (desirable), and relatively large side lobes (undesirable). Other windows achieve different compromises between those two parameters. In FIR filter design, the equivalents to frequency resolution (main lobe width) and leakage (side lobe level) are transition band width and stop band rejection.

You should read the classic paper by Fred Harris: On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.


What is the filter with the less phase shift?

Technically, the most reactive filter is the all-pass filter with a gain of 1, this filter has no phase shift at all. But it is not a really useful filter.

Here's what you need to take in account :

1 - IIR versur FIR

An IIR filter will have less phase shift than a FIR filter for the same cut-off frequency

2 - Cut-off frequency

For an IIR low-pass filter,  a filter with 0.1 Fs cut-off frequency will have a higher phase-shift than a  with a 0.2 Fs cut-off frequency. An all-pass filter i.e. a gain of 1 has no phase shift.

3 - IIR filter order

For an IIR, the higher the order, the higher the phase shift.

4  - FIR number of taps

Assuming you a linear-phase low-pass FIR filter, the more taps you have the higher the phase shift. If you want a high cut-off frequency you can get by with less taps. If you want a really low cut-off frequency you will need more taps thus more phase shift.

5 - Minimum phase-filter

There is a class of filter called minimum-phase filter. They are the most reactive (fastest decay) for a given cut-off frequency according to this reference

https://www.dsprelated.com/freebooks/filters/Minimum_Phase_Filters.html


Assume noise is not a serious issue in your problem. I guess you can get pretty clear speech signal. If you have speech recognition part implemented in your system, I think you should be able to take advantage over the language model in your recognition system. According to the transition probabilities, you shall get some confidence to say at what moment someone stopped talking.

In case you donot know the language model. Here is a simple example. Suppose I want to introduce myself. I will say "Hello[verb] everyone [object]. My name [subject] is [verb] Bob [object]. I [subject] came [verb] from ... [object]". The probability for someone to stop talking right after a subject or verb is much lower than after an object. If you can recognize speech, then you shall be able to take advantage over these language models.

The other thing that is worthy to try is to see whether you can find some stop patterns when someone is talking. What I mean is that from word to word, sentence to sentence, everyone might have his/her stop patterns. If you could successfully retrieve some patterns for example duration, then they shall help you differentiate whether the current stop is alike previous stops or not.


Difference between delayed down-converted signal and a down-converted delayed signal

Let your original signal be (leaving out any constant multiplicative factor):

$$s(t)=e^{i\omega_ct}m(t)$$

Then the down-converted signal is

$$\tilde{s}(t)=e^{i(\omega_c-\omega_{lo})t}m(t)$$

and its delayed version is

$$\tilde{s}(t-\tau)=e^{i(\omega_c-\omega_{lo})(t-\tau)}m(t-\tau)$$

whereas the delayed and down-converted signal is

$$s(t-\tau)e^{-i\omega_{lo}t}=e^{i\omega_c(t-\tau)}e^{-i\omega_{lo}t}m(t-\tau)=
\tilde{s}(t-\tau)e^{-i\omega_{lo}\tau}$$

So the correction factor $e^{-i\omega_{lo}\tau}$ indeed depends on both $\omega_{lo}$ and the delay $\tau$.


Are convolution and multiplexing synonymous?

Multiplexing is the name given to the process of stuffing/switching/packing multiple signals, channel, bits, etc into a single stream/medium/band, whereas the opposite act of separating each individual packed signal, bit etc, back is named as de-multiplexing. In multiplexing each signal (or bits) remains completely intact and have no kind of intersignal interaction other than errors and ditortions whatsoever.

Convolution, on the other hand, operates on two signals and produces a new third signal which is somewhat like the child of the two parents. It's properties are related to its parents by the mathematics of the convolution process. Usually a convolution will be used to produce the output signal of LTI systems where the input signal is convolved with the system impulse response. Rarely a convolution between two arbitrary signals also considered (especially in frequency domain, which is more of a theoretical guide rather than a computational one, circular DFT convolution being an exception). The opposite operation of convolution is named as deconvolution and seeks for finding the driving input from a given output under different constraints.


What Are Intuitive Explanations for Shrinkage in the Context of Image Denoising?

In the image processing field shrinkage means we reduce the power of some components of the signal.
The idea is to have a basis which the signal can represented in in a sparse manner.
For instance, if the signal is harmonic, efficient basis would be the Fourier.
So if we look at harmonic signals and noise in Fourier we expect to see dominant component and other components are the noise, which will be shrinked in order to apply Denoising.
This is the idea behind what's known as Sparse Representations.
In the context of Image Denoising you may have a look at Image Denoising with Shrinkage and Redundant Representations.
In the general context you may have a look at the course Sparse Representations in Signal and Image Processing: Fundamentals.


I can't tell about SciLab but if it is not implemented, you can always get the circular convolution by

IDFT(DFT(x1) * DFT(x2))

where x1 and x2 are your signal vectors and the multiplication is element-wise. If the signals have a different length, zeropad the shorter one to the length of the longer one. This assumes that you want to do circular convolution with a modulus equal to the length of the longer vector.


Beginner level confusion regarding symbols used in modulation

An M-QAM modulation contains M different constellation symbols. Hence, in one symbol you can encode $\mu=\log_2(M)$ bits. (i.e. 6 bits for 64-QAM for example). Now, the bits are equally distributed along the real and the imaginary part of the constellation, i.e. you have $\frac{\mu}{2}$ bits for the real and imaginary part. To encode $\frac{\mu}{2}$ bits, you need $2^{\frac{\mu}{2}}$ different amplitude levels (i.e. 2^3=8 for 64QAM).

Another explanation for the relation of amplitude levels and M: Your M constellation points are arranged in a square. Hence, each side of the square contains $\sqrt{M}$ elements. These are your number of different amplitude levels.

Regarding why the amplitudes have the values of (-7, -5, -3, -1, 1, 3, 5, 7) for 64-QAM. These numbers occur due to two reasons:

the distance between adjacent symbols should be constant (with these numbers the distance is 2). Otherwise, different symbols have different "quality" or error rates. This would be unfair and also does not achieve the optimal distribution in terms of channel capacity.
the constellation should have zero mean (otherwise you would use power for transmitting the mean which does not contain information).

From these two criteria and the addition that you always have an even number of amplitude levels, you actually find that the amplitudes need to be symmetrically distributed around zero. Then, the simplest possibility is to use the values (-7, -5, -3, -1, 1, 3, 5, 7). However, also (-7/2, -5/2, -3/2, -1/2, 1/2, 3/2, 5/2, 7/2) yields a correct 64QAM constellation.


How do we calculate bandwidth from signal formula?

Your intuition is basically correct.  I'm sure they want the answer of 80Hz because that is the bandwidth that you would need to carry these 5 waves, but pure sinusoids that aren't modulated at all can not carry any information.  If they don't carry any information then they are infinitely narrow.  (0Hz rather than 5Hz).  You could modulate them individually to higher frequencies and put them very close together with any loss of information.


Minimum phase systems with pole at infinity

If I see it correctly, your transfer function is

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

This system is not a minimum-phase system because its pole is at infinity (as you've already pointed out). You could use a trivial all-pass filter with impulse response $h[n]=\delta[n-1]$, i.e. a simple delay by one sample to change your original (non-causal) transfer function to a stable and causal transfer function

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

which is also minimum-phase, because it has a pole at $z=0$ and a zero at $z=1/2$. This means that all the system's poles and zeros are inside the unit circle and the minimum-phase requirement is met.


Visualising a Z-transformed Transfer Function?

in continuous-time functions, like $x(t)$, the operation of the derivative makes sense, it is well defined.  the three components to a continuous-time LTI filter are adders (two or more signals being added), scalers (a signal is simply multiplied by a constant), and integrators (the $s^{-1}$ operators).  the first two components do not discriminate with regard to frequency, so with just those two components, you cannot make a "filter" that filters out some frequencies more so than others.  but the integrator does act differently on sinusoids of different frequencies.  higher frequencies come out of the integrator reduced in amplitude more than lower frequencies.

in discrete-time functions (or "sequences"), like $x[n]$, the derivative operator does not makes sense, it is not defined.  the three components to a discrete-time LTI filter are adders (two or more signals being added), scalers (a signal is simply multiplied by a constant), and delay elements (the $z^{-1}$ operators).  the first two components do not discriminate with regard to frequency, so with just those two components, you cannot make a "filter" that filters out some frequencies more so than others.  but the delay element does act differently on sinusoids of different frequencies.  higher frequencies come out of the delay element shifted more in phase than lower frequencies.

analog filters act on physical quantities in time.  digital filters act on numbers.  these numbers are samples of a continuous-time function ($x[n]=x(nT)$) and make up a sequence.  there is no way to do derivatives or integrals directly.  but we can delay any of these sequences by use of computer memory.  that's why we talk of $H(z)$ more so than $H(s)$ when designing and implementing digital filters.


You need to specify your filter design specifications parameters consistently for either an analog or a digital filter. With your posted code, the butterord computes the required order for a digital filter with cutoff frequencies near 1 (which would make sense as Nyquist-normalized cutoff frequencies, but not so much in Hz), then uses those directly to obtain analog filter coefficients.

Depending on your application, you can either design an analog filter with

N,Wn = scipy.signal.buttord(wp=[250.,1000.], ws=[200.,1200.],
   gpass=0.1, gstop=30.0, analog=True)
b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(0, 4, 500))

which should then produce a graph which looks like:

Alternatively you could design a digital filter with

N,Wn = scipy.signal.buttord(wp=[250./fso2,1000./fso2], ws=[200./fso2,1200./fso2],
   gpass=0.1, gstop=10.0, analog=False)
b, a = signal.butter(N, Wn, 'band', False)
w, h = signal.freqz(b, a)

Note however that the original specifications for a digital filter result in an 18th-order filter which may start to show some artifacts due to numerical accuracy. As such I've reduced slightly the gstop parameter to produce the following graph:


Algorithms for deformable image registration

Like Paul suggested, you may want to ask this question at the signal processing stackexchange as well. They may be able to provide you with more information. I am at best a talented amateur in the field.

Let's start with the biggest problem in image registration: there is no catch-all algorithm that just works for everything. Results are heavily dependent on image modality, site and also what you wish to use it for.

Luckily for you, brain mri is a fairly standard application. I would have a look at elastix, which is ITK based and used quite a bit at my department. It also has a nice database of parameter files for different applications, including several for brain mri. A quick glance reveals that most use B-spline registration for the brain.

I have mainly used registration for lung, where my personal experience was that Demons performed best. Note that Demons is NOT one single algorithm, although it is often presented as such. Rather it's a plethora of variations of the same basic idea.

Finally, some commercial software has registration algorithms that perform very well for specific applications. If you're associated with a hospital, ask researchers at your radiology or oncology department if they have access to something that works well.


Obtaining the Impulse response from the frequency of a discrete system

As the frequency response of the system is asked, the RoC woulbe be $|z|>\frac{1}{2}$. Thus, the impulse response is right sided and is given by $$h[n] = \left(\frac{1}{2}\right)^n u[n]$$ We obtain the frequency response by replacing $z = e^{j\omega}$ in $H(z)$. Thus, $$H(e^{j\omega}) = \dfrac{1}{1-\frac{1}{2}e^{-j\omega}}$$


Is the mean of visible color in CIELab color space can be non visible color?

You are implicitly using sRGB, a small color space. You might like to try using scRGB, a similar color space that allows such values.

The average of visible colors should also be visible, however the calculation should be performed in a perceptual color space for the result to be more meaningful; L*a*b* is better than sRGB.

Another question to ponder is whether the centers make sense in L*a*b*. Perhaps you have a bug?


Do waterfalls generate pink (1/f) or brownian (1/f^2) noise?

Interesting question. I juts happen to live right outside a waterfall so I can relate some hands on experience.

The sound depends A LOT and water flow and weather conditions. Water flow is the main engine but wind can generate significant modulations
Depending on conditions you can get strong harmonics in the signal. Sometimes you can actually see standing waves running up and down the fall. These tend to be very low in frequency: maybe around 1 Hz or so.
The fall is pretty wide and has strong directivity. So the spectrum heavily depends where you are.
Overall it sounds more brown than pink
Fun fact: if the fall hits the right resonance, it noticeably shakes the house and the dishes in the kitchen cabinet start to rattle.

This is just for one specific fall and I'm sure others behave quite differently. So all in all any simple model would be a gross over-simplification.


Please see Different way to separate a particular frequency from a signal where I explain the difference and similarity between homodyning the signal and homodyning the coefficients. The DFT is an example of homodyning the coefficients, such that multiplying the coefficients of what would otherwise be a moving average low pass filter by $e^{j2\pi f_c t}$ creates a bandpass response around frequency $f_c$. Twiddle factors are simply the discrete samples of $e^{j2\pi f_c t}$!

Observe specifically in the DFT matrix that the second row steps through all the samples in the roots of unity, so goes around the unit circle once (representing the normalized frequency of one sample per cycle).  The third row steps through every other sample in the roots of unity and in the process goes around the unit circle twice, the fourth row goes around the unit circle
three times, etc. (The first row is all ones representing no rotation, DC, or the base moving average filter). Note that going past $f_s/2$, where $f_s$ is the sampling rate, that the higher positive rotations alias to negative rotations. For example in the 4 point DFT as illustrated below, the normalized frequencies of 0, 1, 2 and 3 are also 0, 1, -2, -1; rotating twice around the unit circle counter-clockwise is the same as rotating twice around the unit circle clockwise when the twiddle factors are just 0, j, -1 and -j: (0, -1, 0, -1). Rotating three times counter-clockwise is the same as rotating once clockwise: (0, -j, -1, j).

The DFT as a filter bank is much clearer if you consider the case of a "streaming DFT" where we shift through a much longer sequence and compute the DFT once after each shift. In this case every point of the DFT would be identical to the output of a FIR filter with the coefficients as given by the DFT row associates with that output.


How to reduce time consuming for calculating pseudo inverse of large matrix in matlab?

Long story very short: No.

If that was possible, Matlab would likely be doing it already; in any case, the Matlab docs for pinv say that the singular value decomposition method is used; that might not be the fastest method, but it's relatively stable for many cases of P, so it's pretty much desirable unless you know exceptionally well what you're doing or your P has special structure (in which case you'd probably not want to use the Moore-Penrose Pseudoinverse in the first place, maybe?).

Longer story: typically, you multiply your Pi onto something to get the inverse operation to P executed.

To cite Matlab documentation itself:

  You can replace most uses of pinv applied to a vector b, as in pinv(A)*b, with lsqminnorm(A,b) to get the minimum-norm least-squares solution of a system of linear equations. lsqminnorm is generally more efficient than pinv, and it also supports sparse matrices.

But whether you can do that depends on the actual reason you're using the Moore-Penrose Pseudoinverse; not all uses are just replaceable with something else.


Several approaches exist to create a vinyl crackle effect [1]. Usually, this type of disturbance is assumed to be an additive one. Hence, the crackling effect is achieved by somehow generating/obtaining a disturbance signal and adding it to the clean signal, not by filtering the original signal.
If you just want to get hold of such crackle signal, you can find at least one generator online, for example here. If you don't mind to pay money for software that can simulate not only vinyl disturbances but a variety of effects, have a look at iZotope Trash 2. (This company used to offer a vinyl sound generator called "Vinyl" for free, but unfortunately this is discontinued.)
Furthermore, there are a number of free samples online (e.g. at freesound.org) that can be used.
If, on the other hand, you're interested in generating a plausible sounding but artificial vinyl click/crackle signal on your own, I would recommend having a look at an article by Välimäki et al., describing, among other things, how to do this.
[1] S. V. Vaseghi, Advanced Digital Signal Processing and Noise Reduction. Chichester, West Sussex: Wiley, 2008.


Sampling rate is only defined within discrete-time systems. So, there's really no sampling-rate equivalent that your integrator has.

What your integrator does have is finite bandwidth. It, for example, will not be sensitive for voltage variations at 2.4GHz.

The bandwidth off that circuit involves the characteristics induced by modeling the circuit based on ideal components, but also on limitations such as unintended (parasitic) inductance and capacitance of components and leads, and, most importantly: by the bandwidth of the opamp itself.

You can buy opamps with only a few Hz of bandwidth, you can buy them with multiple GHz of bandwidth, and very many values in between. Each of these have different strengths and disadvantages, so for a measurement system, choose wisely, based on an understanding of the requirements of your measurement system and signal, and at least a solid understanding of electronics and linear electrical networks; I'm not convinced designing a measurement integrator is trivial.

Generally, you might want to take a step back and actually ask over at electronics.SE how to measure a specific signal, describing exactly and in detail what you are observing, rather than assuming an integrator is the solution here. It might be - but it's equally likely you're not on the right track.


CRLB is usually used to assess the performance for an unbiased estimator because it is seen as the best variance you can achieve. If I understand your question correctly, you are asking if it is fair to use CRLB to assess the performance for a biased estimator? Well, this may not make much sense because you may come across situations where the biased estimator has better variance than the CRLB!

For example, Suppose the true value = 0, you can see one estimate is unbiased with large variance but the other has a bias with smaller variance.


How to estimate the poles of a linear combination of exponentials?

My transform, EST, finds exactly what you are looking for.

The steps of EST are as follows.

Let N be the number of exponentials.

Find the LPC (Linear Prediction Coefficients) for the signal.

For the sake of simplicity I will assume that all the exponential bases (the z values) are complex. The size of the LPC has to be 2N.

You must have at least 4N samples.
Case 1 is when you have exactly 4N samples. Create a system of linear equations with the Hankel matrix and the last 2N samples as the results vector. The LPC vector is the solution to this system.

Case 2 is when you have more than 4N samples. In this case you find the LPC using the covariance method.

Convert the LPC vector to a prediction polynomial, exactly as LSP (Line Spectral Pairs) is doing.
Find the roots of the polynomial.

The roots of the polynomial are the exponential bases you are looking for.

If the signal is noise-free, then the exponential bases are precise. If not, it is an estimate.

Once the exponential bases are known, finding the coefficients (the c values), is simply a matter of linear equations.

You also asked if you can find the number of exponentials.
In the case of a noise-free signal: yes.
Otherwise, you should find a large enough number of exponentials and filter out some of them.
Please note that each exponential base holds a frequency and the growth/decay rate of the signal at that frequency and each coefficient holds an amplitude and a phase. These parameters can be used for filtering criteria.


Here is a potential starting point: it is considered  that when a function behaves as a power law  ($x(t) \sim t^\alpha, t> 0$), then its magnitude spectrum as well, with $|X(f)| \sim f^{-(1+\alpha)}$. One can see for instance:

Computing Fourier transform of power law
An Interesting Fourier Transform - 1/f Noise
Sketches of a "proof" and references (since there are  convergence issues) are provided in an answer to I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$.

If I plug that in your case, since $|X(f)|^2 \sim f^{-5/3}$, so we seek $-5/6=-(1+\alpha)$. I obtain $\alpha = -\frac{1}{6}$. So you would have $x^2(t) \sim t^{2\alpha}$, with a magnitude spectrum $|\chi(f)|$ in $f^{-(1-\frac{1}{3})}$, or $f^{-\frac{2}{3}}$.
Now, is the discrete signal sufficiently well-sampled and clean to verify it in the discrete domain is another question.


Your specification of the particulars for your situation aren't clear to me, but I'll explain how you can arrive at the answer. The DFT is defined as:

$$
X[k] = \sum_{k=0}^{N-1}x[n] e^{-j2 \pi nk/N}
$$

The (normalized) frequency associated with the $k$-th bin is $\frac{2\pi k}{N}$. Therefore, the spacing between the $k$-th and $(k+1)$-th bins is $\frac{2\pi}{N}$. Recall that the conversion between unnormalized frequencies (i.e. before sampling) and the normalized angular frequency $\omega$ used in the DTFT is:

$$
\omega = \frac{2\pi f}{f_s}
$$

That is, a time-referenced frequency $f$ before sampling maps to $\omega$ in the normalized frequency domain, where $f_s$ is the sample rate of the signal (i.e. how many samples are taken of the signal per second). You can then map the bin spacing above into unnormalized frequency:

$$
\Delta f = \frac{2\pi}{N} \frac{f_s}{2\pi} = \frac{f_s}{N}
$$

This simple relationship just says that the spacing between bins is equal to your signal's sample rate divided by the DFT length. It's unclear from your question whether the interval between each sample is 3.5 microseconds or whether your block of 256 samples covers an overall duration of 3.5 microseconds. Based on whichever is true, you can calculate the appropriate value.


Is there a method to compute an average amplitude for 2D sinusoids using the F.C.'s?

Easy thing first: You might want to look into spatial statistics

A fast library to compute Spatial Statistics in Matlab using FFTs.
A good description of Spatial Statistics from a past student

Long Stuff:

The trouble with comparing materials science information (eg AFM) is that the piece of material was measured at a position in space and time.  It is difficult to quantitatively compare images in one place and an image in another place.

When you take the Fourier transform, you are measuring the frequency of that material at a position in space and time.

You can remove the spatial dependency by taking a convolution of the Fourier transform with itself.  This process is akin to Point Distribution Functions, N-Point Statistics, or Joint Probability Density Functions.  The Spatial Statistics Matlab library can be applied to most image data.


Found the answer in numpy documents for fft:

# python to perform dft
# from import numpy.fft import *

A = fft(a, n)

  A[0] contains the zero-frequency term (the sum of the signal), which is always purely real for real inputs.

  A[1:n/2] contains the positive-frequency terms

  A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

  For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency, and is also purely real for real input.

  For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency.

There are many ways to define the DFT, varying in the sign of the exponent, normalization, etc. In this implementation, the DFT is defined as

$$
A_k = \sum_{m=0}^{n-1} a_m exp \left\{ -2 \pi i \frac{mk}{n} \right\},\ \ \ \ k=0,...,n-1
$$

The DFT is in general defined for complex inputs and outputs, and a single-frequency component at linear frequency  is represented by a complex exponential , where  is the sampling interval.

The values in the result follow so-called “standard” order:

Correction to code above:

# python DFT Sin Test
# Bill Moore
# 4/11/2019

import sys
import numpy as np
import matplotlib.pyplot as plt

def dsp_mag_and_phase(DIN, N, filename):
    """
    Cleanup and plots DFT and FFT results

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FREQPLOT")

    XK = DIN
    K  = range(0, N)

    # Swap Positive and Negative Frequencies of DFT/FFT
    #  (to get standard spectrum of -pi to pi)

    # Assign a positive or negative frequency to each point
    KK = []
    for k in K:
        if k>=0 and k <= int(N/2):
            KK.append(k)
        else:
            KK.append(k - int(N))

    #print("KK:", KK)
    #print("len(KK)=", len(KK))
    #sys.exit()

    AM = np.zeros(N)
    AP = np.zeros(N)
    AM = np.zeros(N)
    AP = np.zeros(N)

    #Positive Frequences
    AX_pos = XK[0:int(N/2)+1]
    AK_pos = KK[0:int(N/2)+1]

    #Negative Frequences
    AK_neg = KK[int(N/2)+1:]
    AX_neg = XK[int(N/2)+1:]

    # Concatenate negative then positive
    AX = np.concatenate((AX_neg, AX_pos))
    AK = np.concatenate((AK_neg, AK_pos))

    # Convert Complex Frequency into Magnitude and Phase
    for k in range(0, len(AK)):
        AM[k] = np.sqrt((AX[k].real)**2 + (AX[k].imag)**2)
        AP[k] = np.arctan2(AX[k].imag, AX[k].real)

    return (AK, AM, AP)

def dsp_dft(DIN, N, filename):
    """
    plots N-point DFT

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_DFT_MAGNITUDE")

    DFT = np.zeros(N, dtype=np.complex)
    K   = range(0, N)
    w0  = 2 * np.pi / N #Fundamental Frequency (Radians)

    # DFT
    for k in K:
        DFT[k]=0
        for n in range(0,len(DIN)):
            DFT[k] = DFT[k] + DIN[n]*np.exp(-1j*k*w0*n)
        DFT[k] = (1/N)*DFT[k]

    # Covert Complex to Magnitude and Phase
    (AK, AM, AP) = dsp_mag_and_phase(DFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "DFT", filename + ".png")

    return (AK, AM, AP)

def dsp_fft(DIN, N, filename):
    """
    plots N-point FFT
    (FFT is just a computationally faster DFT)

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FFT_MAGNITUDE")

    # FFT Points equal to length of DIN (default for fft function)
    if N == 0:
        N = len(DIN)

    FFT = np.fft.fft(DIN, N)

    (AK, AM, AP) = dsp_mag_and_phase(FFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "FFT", filename + ".png")

    return (AK, AM, AP)

def plot_stem(x, y, title, filename):
    print("=> PLOT_STEM")
    fig = plt.figure()
    plt.title(title)
    plt.stem(x, y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def plot_mag_and_phase(k, mag, phase, title, filename):
    print("=> plot_mag_and_phase")

    fig = plt.figure()
    plt.clf()

    plt.subplot("211")
    plt.title(title)
    plt.ylabel("Magnitude")
    plt.plot(k,mag)

    plt.subplot("212")
    plt.ylabel("Phase")
    plt.plot(k,phase)

    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def dsp_range(sampleN, sampleF):
    """
    return  : time points for sampling at a given frequency

    sampleN : Number of Samples
    sampleF : Sample Frequency in Hz
    """
    sampleT       = 1 / sampleF
    time_range    = np.linspace(0.0, sampleN*sampleT, sampleN)
    return time_range;

def run():
    print("=> RUN")

    sampleF   = 1E3
    sampleN   = 100
    Fsin      = sampleF/4

    # sample function
    tvec       = dsp_range(sampleN, sampleF)
    y1         = np.sin(2*np.pi*Fsin*tvec)
    plot_stem(tvec, y1, "samples", "plot_test0.png")

    dsp_dft(y1, 100, "plot_test1")
    dsp_fft(y1, 0,   "plot_test2")

    sys.exit()

run()


resonance demodulation for bearing faults

This presentation has a pretty good overview of characterizing bearing faults. In essence, once the measured frequency content does not match what is considered a good response from a bearing, it can be flagged as faulty or nearing fault.


What do the Contour Levels Produced by Matlab Function Contourc Represent

So you:

  I constructed a picture with 10 rectangles, each of which has a distinct depth starting from 0.1 ending up at 1.0

and then:

  I ran the code on 0.1 contour level

The matlab function is doing the right thing. A rectangle of height 0.1 has a contour at 0.1.  A rectangle of height 1 also has a contour at 0.1 (and 0.05 and 0.5, and 0.9999999, but not at 1.000001).  The rectangle exists from 0 to 1 (in this case).


How is the Tx/Rx Gain calculated in USRP?

The term "gain" is pretty appropriate. It's a number describing the output power relative to some arbitrary base power.

So I would like to know how this transmitter channel gain process is implemented,

Depends on the USRP (or if you're using a modular model, on the daughterboard), but speaking generally: A mixture of selectable/bypassable amplifiers, variable gain amplifiers and adjustable attenuators.

how it is calculated,

Not calculated. Set by configuring the things mentioned above.

how it works in practice (since a particular gain value doesn't correspond to a particular power level).

It sets how much stronger or weaker your amplification is. Simple! You add 3 dB of gain, you get twice as much output power. Of course, that's approximate, and only within the linear range of all the signal chain.


Impulse response of a system in z domain

You are right.

In the solution, the second line which performs the long division

$$ H(z) = -4 + \frac{ 5 + \frac{7}{2} z^{-1} }{1 - \frac{3}{4}z^{-1} + \frac{1}{8} z^{-1} } $$

is wrong and should be corrected as:

$$ H(z) = -4 + \frac{ 5 - 3 z^{-1} }{1 - \frac{3}{4}z^{-1} + \frac{1}{8} z^{-1} } .$$

However, the partial fraction expansion at the following third line is correct; indicating there was a typo in the second line. Finally, the impulse response (which is based on line three) is also correct.


Typically a known sequence is transmitted (often as a preamble) that has properties to aid in determining the channel characteristics (by spreading its energy across all frequencies within your signal bandwidth) and in determining the symbol clock (many transitions at the symbol clock boundary). Both can be accomplished by using pseudo-random sequences which also have the desirable property of an autocorrelation that approximates an impulse function (useful for establishing a timing boundaries under low SNR conditions). When shorter sequences are sufficient for  acquisition and channel estimation, Barker codes also have the desirable autocorrelation and signal energy spreading properties mentioned (these go hand in hand: the autocorrelation of white noise is an impulse).

The details on estimating the channel (or equalizing your signal to remove the channel distortions) for single-carrrier modulation schemes are identical to the operations detailed in this post here:  Compensating Loudspeaker frequency response in an audio signal


How to turn off calculating pyramids during image processing?

Pyramids can be toggled on and off using the Erdas preferences: File > Preferences > File Formats > Image Files (General).


How to remove wind noise from video file

Extracting Audio from Video

There is a simple cross-platform tool for manipulating / accessing video & audio files called ffmpeg. I've linked you the documentation available for the whole toolkit, just so you can glance around and get a feel for what sort of wonders can be performed with it.

But, all you need to do is extract audio from your video file. You can do this with ffmpeg and a pointed explanation of how to do exactly this, following this link. It's more of a tutorial / user's guide and less overwhelming than the documentation I linked above. You have two choices...

Extract the audio straight from the terminal:

$ ffmpeg -i C:/path/to/file/your-video-file.mp4 -vn -acodec copy your-audio-file.aac

Remember to change the file extensions depending on what video file extension you're working with.

Extract the audio from the terminal through Python:

import subprocess as subproc
command = "ffmpeg -i C:/path/to/file/your-video-file.mp4 -vn -acodec copy C:/path/to/your-audio-file.aac"
subproc.call(command, shell=True)

Removing Wind Noise By Analyzing FFT Data

Built-in tools for audio noise reduction / removal, like the Noise Removal Tool in Audacity, are all various implementation of the same concept: spectral subtraction. Spectral subtraction is so widely used because... it's relatively simple. All it takes is to compute the FFT of the audio signal in question, split the FFT into smaller "frames", identify the power level of the steady background noise within each frame, and then filter out that energy from the entire frame. Coming up with what filter, what frame size, and what energy level works best for your case is part of what makes each spectral subtraction technique different. So, there's no one, static solution for all audio files that I could show you - just analyze the magnitudes of your audio's FFT frame-by-frame and figure out what works best for your audio. In the end, you're left with whatever additional energy is present in the signal, which is usually the part of the interesting part of the audio.

If you don't know how to compute & view the FFT of a signal in Python, consult the Examples section of the numpy documentation on the numpy.fft.fft function.

Patching Together Reworked Audio with Original Video

Again, ``ffmpeg` is your friend here. This is as simple as a one line command:

ffmpeg -i your-video-file.mp4 -i your-audio-file.aac -c:v copy -c:a copy output.mp4

I hope that this approach helps give you the DSP experience you're looking for, while also solving the noise reduction problem you have with your current video. Good luck!


Is there any relationship between the FFT and linear prediction?

One thing that comes to mind is that a linear predictor can be implemented using the LMS algorithm, and the LMS algorithm can be implemented in block fashion using FFTs (see Chapter 7 of Haykin's Adaptive Filter Theory). I wouldn't say that implies any kind of relationship, but FFTs can be used to accelerate the implementation of some linear predictors.


How to construct Neyman Pearson statistic to detect discrete signal in AWGN with unknown location?

I'd suggest trying to get rid of the $p$ index. You can write $\mathcal H_1$ as $Y_i = W_i + A$, then you can do the detection on each sample. You put the self-study tag, so I'll let you do the math but it'll come down to detecting the pulse if $y_i > \frac{A}{2}$. This will only work, if you know $A$ of course. You'll need a fairly high SNR for this to work well though. If the pulse is more obscured by the noise, you can try the next approach.
Another way to approach this problem, is to use filtering. Since you know you're looking for a rectangular pulse of length $100$, you can run the noisy data through a matched filter which will at least give you some processing gain (https://ccrma.stanford.edu/~jos/sasp/Processing_Gain.html). If you plot the samples after the filtering, you'll see a peak somewhere. Since the matched filter will have some delay, the start of the pulse is $\text{Peak Index}-\text{Pulse length}+1$. This is another way to do the processing, especially if you are able to gather a batch of data then do the processing.
Edit
The first approach requires a single scan through the input data stream leading to $10^5$ comparisons, $\mathcal{O}(N)$ operations where $N$ is the length of the input data.
The second approach involves using a FIR filter which can be done in $\mathcal{O}(NM)$, where $M$ is the length of the filter.
The second approach can also be implemented using FFT filtering which can be done in $\mathcal{O}(N\text{ log}N)$ operations.


Phase shift is often confused with time delay, although related these are two different quantities. Phase shift specifically refers to the rotation of a complex number. For example given the complex number $x = Ae^{j\theta}$ which through Euler's formula has real and imaginary components given as $A\cos(\theta)+jA\sin(\theta) = I + jQ$, a phase shift of $\phi$ would rotate $x$ according to:
$$ x e^{j\phi}= Ae^{j(\theta+\phi)}$$
Resulting in real and imaginary components as:
$$xe^{j\phi}= A\cos(\theta+\phi)+jA\sin(\theta+\phi)$$

Time domain and frequency domain waveforms can be represented as complex numbers where each sample has a magnitude and a phase. A phase shift would be a rotation of those samples, in either domain (so we can have a phase shift in time, or a phase shift in frequency). A constant delay in time has the property of a linear phase in frequency: the phase shift versus frequency will be linearly increasing in the negative direction for a fixed delay in time. For instance, given the observation of a sinusoidal tone at the input and output of a cable (fixed delay), for an extremely low frequency we would see the result of a very small phase difference between the sinusoid measured at the output compared to the input, but as the frequency is increased, this measured phase would increasingly larger in direct proportion to the increase in frequency.
For further related details, see:
What is an Intuitive Explanation of the Phase of a Signal
How does a shift in time domain result in phase shift in frequency spectrum?
Meaning of Phase response of a filter? In simple words?


Why does the Hilbert-Huang transform not suffer from the uncertainty principle, and what are the implications of this?

I'm not that familiar with the Hilbert-Huang transform, so I won't comment on that.

You seem to be under the impression that the DFT "suffers from the uncertainty principle".  This is not true.  Don't feel bad about thinking that as it is a widespread misunderstanding.  This is a quote from an email to me from a well known expert in the field after I tried to get them to look at my exact frequency formula for a single (non-integer) pure real tone in a DFT:

  An "exact" computation of frequency based on FTs is not possible for
  some very simple reasons:  The Fourier Transform is itself an
  estimator, and the time-frequency uncertainty principle means
  frequency can't be determined exactly without an infinitely long
  observation window.  So essentially all frequency computations of time
  series or signals are estimates, it's just that some are more accurate
  or more efficient than others.

  So claiming that you have an exact solution just isn't going to get
  the attention of everyone interested in the field.

The first paragraph is entirely hogwash.  The last sentence was prophetic, hence I started my blog so no one could dispute the math.  You will find my original exact solution here:

Exact Frequency Formula for a Pure Real Tone in a DFT

I have discovered, and written about, several more.  You will find them among my articles.  There are also a few in the pipeline, including a correction to Macleod's three bin formula which makes it exact also.

The FT does apply to the uncertainty principle in physics, but in a different manner.  There they are talking about the transform of a Gaussian curve which happens to be an eigenfunction of the FT (but not the DFT, close though) and how it is the "narrowest function" in the transform.

I'm not an expert in Physics either, but I dabble.  I consider this to be my greatest mathematical discovery ever (exclipses the simple math of the frequency formulas by a mile):

Off Topic: Refraction in a Varying Medium

This is the vector form of Snell's law in a continuously varying index of refraction medium.  It is a theoretical particle model so its applicability to light is still a subject of study for me.  I believe it has the potential to illuminate subatomic behavior (i.e. the quantum world) but I haven't gotten there yet.


This is known as deconvolution. A typical approach would be to apply a Fourier transform to both sides
$$
y = x \circledast h \overset{\mathscr{F}}{\Leftrightarrow} Y=XH
$$
and divide both sides by the known $X$ to obtain $H=Y/X$. Then apply an inverse Fourier transform to both sides, yielding $h$.

The linked wikipedia article goes into more detail on variations of this approach.


So is there such a thing as a DSP engineer nowadays?

Yes there is. There definetely is there such a thing called DSP  Engineer, but it is not exactly what the article described, i think it even goes more high level.

Nowadays, most companies relies on already established Knowledge such as standards, or papers from IEEE or AES, or such places. So a DSP engineer requires now more than ever, a lot of knowledge about Programming, in a low level fashion, programming in Assemly, C/C++ to producte robust, reliable, mantainable, understandable and expandable code. It is true what the artcile says that the are cookbooks, so if you don't remember exactly how LMS algorithm goes, you can pick up the book and refresh it in 5 minutes. But the concept needs to be there. If you are completely new to DSP, but you know a lot about programming and Kernel and operating systems, you will have a hard time implementing algorithms if you don't know the basics.

There is no real answer to your question but i state again my first comment. Yes there is, but it has incorporated a lot more about embedded Systems that what it had bewteen 1960 and 1980.


Why noise canceling headphones do not attenuate the signal itself?

The microphone will record the noise on the outside and the speaker sits on the inside with a bit of passive sound shielding in between. Outside noise will reach the microphone unattenuated, but the actual signal will be attenuated by the passive shielding. The noise that reaches the speakers attenuated by the passive sound shielding and the circuit of the active noise canceling will simulate the effect of the passive shielding to generate the right amount of anti-noise to add to the actual signal for playback on the speaker. The circuit can also simulate how much of the actual signal reaches the microphone and compensate for that or it can just ignore that issue because this effect would only cause a small attenuation of actual signal.

Example:
Lets say our passive shielding reduces sound levels of signals traveling through it by 75%, so 25% reaches the other side. Actual signal on the inside without noise canceling would be 1.0, noise on the outside would be 1.0 as well. Without noise canceling noise level on the inside would 0.25 and actual signal on the outside would 0.25. So the microphone would record 1.0 noise + 0.25 noise, now when we turn on active noise canceling, as a first order model, the circuit would simulate the attenuation of the shielding and subtract 0.25 noise and 0.0625 actual signal. If the attenuation is simulated perfectly, on the inside we would get 0.9375 actual signal and no noise. SNR is significantly improved. Of course this is just a first order model, the attenuation provided by the passive shielding does not have a flat frequency response and the microphone has also some electrical noise, that will be added to the speaker signal. If the simulation of the passive attenuate is not perfect, we will also add extra noise. We also need to adjust for the frequency response and gain of microphone and speaker. Small errors can easily cause us to add more noise to the signal or subtract to much from it.


Resampling time series to regular array, then downsampling (Butterworth)

I am not sure why you concluded "A Butterworth Filter" is the way to go, and looks like you have come to the same question- what originally led you to that?

A Butterworth Filter implementation stems from copying analog filter techniques and from what I have learned with DSP is that you have the wonderful opportunity to NOT copy the analog. This is because analog design techniques have been derived from the mathematical limits of what standard R L and C components can give you, while in DSP you can more likely address the problem from an optimum mathematical perspective within the greater space of delays, adds and multiplies (to note and hopefully not confuse; delays, adds and multiplies are NOT strictly limited to digital implementations; for example SAW filters). This may be a good case in point if we look closer at why you would choose a Butterworth filter versus what optimized algorithms such as least squares (FIRLS in Matlab) can give you. The least squares algorithm will result in a linear phase FIR solution, with a stop band roll-off which is ideal for multi-rate resampling applications. fred harris has made a good case (I cannot find the specific reference) that for resampling applications, and FIR filter will always out-perform an IIR filter.

Further what is really great about FIRLS for multirate applications is that you can easily implement a multi-band filter, concentrating the rejection bands to be just at the alias locations, as in the example shown in the figure below showing the same filter design for either a decimation (DFIR) or interpolation (IFIR) design, in this case using a factor of 4:

To answer your question about upsampling: If the signal has been upsampled, this is typically done by inserting zeros (which creates all the aliases in frequency) followed by a low pass filter which will eliminate the aliases and in the process grow the zeros to their optimum interpolated value, so yes in this case, as long as the intepolation filter is present, you will not need an additional low pass filter.

I show this in the following figure showing that if you do the interpolation first in a rational rate conversion (both interpolate and decimate) you can combine the IFIR and DFIR designs as one filter.

Also for resampling given its simplicity, consider using a CIC structure. If you are concerned about passband roll-off- the inverse Sinc compensator is quite simple (3 taps). Both are detailed in the following links:

Decimation concepts: Fast Integer 8 Hz 2nd Order LP for Microcontroller

Upsampling Spectrum: Upsampling and comb function

Inverse Sinc Filter: how to make CIC compensation filter

More about CIC: CIC Cascaded Integrator-Comb spectrum


How to practically correct for group delay imposed by a digital filter (Python)

A causal filter such as what must be used in real time filters will always have group delay. Because it is unavoidable in such a circumstance it’s beneficial to know what it is.
Group delay cannot be corrected for in real time, as it would require knowledge about signal which hasn’t happened yet.  Variable group delay can be made constant by applying a filter with a phase response such that the sum of the applied filter and the incident filter is linear, but this will result in a constant, almost certainly non-zero and positive group delay.
In rendered applications, a method such as filtfilt will apply the filter, then flip the filter about zero time and apply it again, which will square the magnitude responds and result in no phase shift. There should not be anything about this process that results in stability problems.
Hypothetically, one could construct a zero phase version of a filter in a rendered application by way of some maths. One way would be to take the frequency transform of the resulting signal and subtract the phase response of the filter, and transform it back to the time domain.  This might have some gnarly edge artifacts.  You might try creating an all pass filter with the negative of the filter’s phase response, use that to filter the signal.  That would result in a signal that’s delayed by half the filter length, which you could remove manually. However, it sounds like you are using an IIR filter, so picking an appropriate length for prototyping the all pass might be tricky.


Deriving the Langrangian interpolation polynomials in Cook-Toom convolutions

If you represent a second-order polynomial $s(x)$ with Lagrange polynomials $L_i(x)$ and interpolation points $\beta_i$, $i=0,1,2$, such that

$$s(x)=s(\beta_0)L_0(x)+s(\beta_1)L_1(x)+s(\beta_2)L_2(x)\tag{1}$$

then for the equality in $(1)$ to be satisfied, the polynomials $L_i(x)$ must have zeros at $x=\beta_j$, $j\neq i$, and they must equal $1$ at $x=\beta_i$. This completely determines the polynomials $L_i(x)$.

E.g., the polynomial $L_0(x)$ must have two zeros at $\beta_1=1$ and $\beta_2=-1$, respectively, and it must satisfy $L_0(\beta_0)=L_0(0)=1$. The latter condition just determines the scaling. This results in

$$\begin{align}L_0(x)&=\frac{(x-\beta_1)(x-\beta_2)}{(\beta_0-\beta_1)(\beta_0-\beta_2)}\\&=-(x-1)(x+1)\\&=1-x^2\end{align}$$

The other polynomials can be derived in the same way.

Note that the polynomials $L_1(x)$ and $L_2(x)$ have a scaling factor of $\frac12$, which is later in the derivation absorbed in other constants for reasons of computational efficiency.


Even & odd refer to the symmetry around $n = 0$.

Even means $x[n] = x[-n]$; you can get the part for $n < 0$ by simply mirroring the part for $n > 0$ at the $n=0$ line.

Odd means $x[n] = -x[-n]$; you can get the part for $n < 0$ by simply mirroring the part for $n > 0$ at the $n=0$ line and multiplying it by $-1$.

A cosine wave is even, sine wave is odd.

These are all just special cases of the general symmetry that says

  if it's real in one domain, it's conjugate symmetric in the other.

Conjugate symmetric means that the real part is even and the imaginary part is odd. Most people know that a real time domain signal as a conjugate symmetric spectrum, but it also goers the other way around: a conjugate symmetric time domain signal has a real valued spectrum.


Normalizing audio waveforms code implementation (Peak, RMS)

What is the best method to solve this?
  Intuitively, it seems I might need to calculate a local max or average based on a moving window (rather than the entire set) but I'm not entirely sure. Help?

What you are describing with your second of these two quoted bullet points is Automatic Gain Control, which eventually will take us to compression. What you might find useful in your case, i.e for visual purposes only, is a simple logarithmic transformation which is very close to what compression is about but without the complex controls it demands.

The most common way to solve such problems in the audio world is to use a compressor.

The compressor is monitoring the amplitude of the incoming signal. If that amplitude crosses a user-defined threshold it decreases the amplification in an attempt to keep the signal within certain limits. As you can understand, imposing such limitations to a signal without special care could result in annoying audible artifacts. For this purpose, the compressor has an additional set of controls called "Attack" and "Release". These are two parameters that control how fast does the compressor impose its limitation once the threshold has been crossed (Attack) and how fast does the compressor return to the initial amplification of 1:1 after the amplitude of the signal has returned below the predefined threshold. These three parameters, the threshold and attack and release times are what make the compressor "complex" in setting it up properly. Extreme values in these settings can result in more artifacts.

You may have "heard" a compressor working in radio shows. It is the effect that you hear when a speaker suddenly stops talking and after a very brief period of time the background sounds seem to rise in volume. When the speaker starts talking again, the background sounds seem to be diminishing abruptly again and the speaker's voice now dominates the dynamic range of the signal. That's the effect of a compressor.

The main idea here is that the compressor is making soft sounds louder and loud sounds softer.

  Note: The waveform is purely for visual purposes (the audio will not be played back to the user).

If you are not looking to actually apply dynamic range compression and you simply want to make sure that your viewer can see all details within a signal, you can apply a logarithmic transformation, which is very similar to how a compressor / limiter would operate but without the complexity of setting it up.

Given a signal $x(n)$, you can generate a signal $x_{transformed}$ with:

$$ x_{transformed}(n) = c \times log(1+x(n))$$

Where $c$ is a normalisation parameter. Due to its shape the $log(\cdot)$ function will boost small amplitude values and limit high amplitude values. To get a feel about how the waveform might look like, you might want to have a quick look at the way Audacity does its 'dB View' of waveforms.

You might find both the compression and logarithmic transformation concepts useful to what you are trying to achieve.

Hope this helps.


Generating drum sample in programming language

There are different approaches to generate artificial (I guess that's what you're aiming at here) drum sounds. The approach that you describe is certainly a good way to start - i.e. look at an existing sample that you're trying to approximate and find the properties that make up it's main characteristics. For example

is it tonal (e.g. like a tom or bass drum) or noise-like (e.g. like a high-hat) or a mixture (e.g. like a snare)?
what's the frequency/spectrum of the tonal part
what's the bandwidth of the noise part
what's the attack / decay behaviour
etc.

Knowing these important aspects of many percussive sounds, you can try to modify - e.g. - a sine wave by changing the frequency and amplitude over time, adding (filtered) noise and so on. Some experimenting will yield sounds similar to that of the classic TR-808. This machine may also be a nice source of inspiration on how to approach drum modelling using basic synthesis elements.

Another - somewhat more involved - approach is called physical modelling. This method starts from trying to emulate and - as the name says - physically model the elements that are involved in actually generating a sound. In the case of drum synthesis this comprises, for example:

the drum head
the stick that is used to hit the head
the drum shell.

A thesis on exactly this topic can be found here.

Resources on how to actually program these models can for example be found in the communities of Pure Data, Csound or SuperCollider as these are programming environments that have been specifically designed for audio/music signal processing (for example, see here for a short tutorial on drum synthesis in Pure Data).

For some theoretical background, I have been recommended the musimathics books, but have not read them personally.


Implementation of gain scheduling for a PI controller

if you want the transition between two sets of controller parameters to be smooth, you need to apply the integral gain before the integrator. The controller output will perform a step otherwise, which is the exact opposite of seamless.

Your speaking of a PI controller and you have to apply a suitable strategy for the proportional gain too. Possible solutions are, but not limited to:

apply a 2nd order filter to the controller parameters to ensure smooth transistion, this might affect controller stability during transistion but this shouldn't be a problem with your simple control loop.
Implement a controller per parameter set, one of the controllers is active while the others are in tracking mode, switch over when another set of parameters is required. This will ensure, that the controller runs instantly with the correct parameters which typically doesn't destabilise the control loop. (Simulink PID Block, search for tracking)
Try to find a function that expresses the perfect controller parameters based on the current load situation and adjust the parameters accordingly. This might be tedious though.


Conversion from Floating point representation to Fixed point representation

I think you need to use the floor() or round() functions in MATLAB, to emulate fixed-point variables and operations.  So you have to know the range of the fixed-point value and the precision of it.  The ratio of the range to the precision is the dynamic range and you get 6.02 dB and one bit of word width every doubling of that ratio.

If you convert:

x_fixed = precision*round(x_float/precision);
if x_fixed > fixed_max
    x_fixed = fixed_max;
elseif x_fixed < fixed_min
    x_fixed = fixed_min;
end;

then fixed_max - fixed_min = (2^N - 1)*precision  for N being the number of bits in your fixed-point word.  If you do this after every mathematical operation, you will be emulating fixed-point arithmetic.

there might be better ways to round x_float than i shown.  sometimes the floor() is better than round().  and "fraction saving" is a better way to round an intermediate result with more bits to your N-bit fixed-point word, but requires an error state for each variable.

for every variable and do your arithmetic with their fixed-point values and round the result to a fixed-point value with code like above.

that's how to convert floating-point MATLAB to emulating fixed point math.


Filtering - multiplication in frequency domain

To expand on @Maximilian Matthé's answer, you can visualize the spectral leakage effects (convolution of by sinc function in the frequency-domain) by zero-padding the inputs used in calculateFFT. For example, the following function zero-pad the inputs to a length $k$ times longer than the original input (where in this instance $k=4$):

def calculateFFT(time,signal):
    k=4
    N=k*len(signal)
    df=1/(k*(time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [k*2.0/N*abs(i) for i in fft(signal,N)[0:int(N/2.0)] ]
    return frequencies,fftValues

Plotting the result you should see that the low-frequency components actually overlap (which indicates that the filtering does not actually increase the signal there):

So, where did those ripples around the spike came from? It turns out they were always there, but when computing the FFT you were obtaining the value of the spectrum at a discrete set of frequency values and those ripples happened to pass through exactly zero at those frequencies. If you had chosen a slightly different signal frequency that is not an exact multiple of the FFT frequency resolution (0.2Hz in your case), for example 1.25Hz instead of 1.2Hz, the sampling of the frequency spectrum would have looked quite different (due to the frequency sampling at different points in the oscillation of the ripples):


Derive chernoff bound for $\mathrm{erfc}(x)$

I don't know why you are approaching the problem the way you are by going to a two-dimensional approach.

The Chernoff bound starts with the observation that for any random variable
$X$ and any $\lambda > 0$,
$$P\{X \geq a\} \leq E[\exp(\lambda(X-a))]\tag{1}$$
because $\exp(\lambda(x-a))\geq  \mathbf 1_{\{x\colon x \geq a\}}$.
Consequently,
$$P\{X \geq a\} \leq \min_{\lambda > 0}E[\exp(\lambda(X-a))].\tag{2}$$
For a standard Gaussian random variable, $E[\exp(\lambda(X-a))]$
can be found easily by using a method called "completing the square"
(in the exponent which as terms like $-\frac{1}{2}x^2+\lambda(x-a)$ in
it) and noting that the resulting integral is the integral
of a Gaussian density times a constant. Alternatively, you can
use the moment-generating function and read off the answer.
After minimization with respect to $\lambda$, the result is
$$P\{X \geq a\} = Q(x) \leq e^{-\frac{x^2}{2}}.\tag{3}$$
The two-dimensional approach leads to the slightly better result
$$Q(x) \leq \frac 12e^{-\frac{x^2}{2}}.\tag{4}$$
Note that you have proved already that for $x > 0$,
$$P\{X_1>x, X_2>x\} = Q^2(x) \leq \frac 14 e^{-{x^2}}$$
from which $(4)$ follows immediately, and $(4)$ is not the
Chernoff bound on $Q(x)$.


What does it mean for an image to be "Markovian"?

It's a probabilistic model that is assumed to generate the image. To make theoretical computations on the performance of certain image processing tasks such as DCT transforming the image for data compression, the image is modelled with a simple mathematical equation.

The Markov image generation model, is used to produce images where there will be high enough correlation between neighbouring samples. This high correlation is therefore efficiently utilized by a transform such as DCT to achieve compression ratios that is close to the maximum available with a nonpractical KLT (Karhunen Loeve Transform)

As far as I know, it is not much related with markovian fields.
The simplest 1-D example of such a markov source model, assuming linear dependence between samples, is like: $$x[n] = \rho x[n-1] + {\omega}[n]$$

This simple dependency clearly shows the fact that current sample x[n] is produced from past sample x[n-1] plus a white noise sample w[n], where $\rho$ is the correlating factor and as it is close to 1, samples are highly correlated, when it is close to zero, samples become more noisy and uncorrelated.

For a real image, such a model can only predict samples in a very short duration, because of the non-stationarity of the image. This is one of the reasons why DCT based transforms are applied for short lengths of blocks such as 8x8 pixels.


How can I determine the complex envelope of a signal when knowing the analytic signal?

With constant amplitudes and phases, this becomes a simple exercise:

First, go to the spectral domain
$$
X(f) = \mathcal{F}\{x(t)\}=\sum A_i \frac{1}{2j} (\delta(f-(f_c+f_i)-\delta(f+(f_c+f_i))\exp(j2\pi\tau_if)
$$
with
$$
\tau_i=\frac{\theta_i}{2\pi (f_c+f_i)}
$$
being phase offset translated to the equivalent time shift.

Now, you can shift this such that the carrier becomes DC:

$$X_d(f)=X(f-f_c) = \sum A_i \frac{1}{2j} (\delta(f-2f_c-f_i)-\delta(f+f_i))\exp(j2\pi\tau_i(f-f_c)$$

Applying a low-pass filter kicks out the high-frequency Dirac at $2f_c$. Then, transforming back to time, you get

$$
x_d(t) = -\sum A_i\frac{1}{2j}\exp(j2\pi f_i(t-\tau_i))\exp(-j2\pi\tau_if_c).
$$

This should be your complex envelope.

In case yo want to do this in Matlab, I would go for the following approach (which is completely different from the math above. It is the straight-forward implementation of shifting in frequency by multiplication in time):

T=1e-6;
fs=1e11;
t=[0:1/fs:T]';
fc=11.5e9;

B = 0.5e9;  % Signal bandwidth

carrier = exp(2j*pi*fc*t);

xpassband=1*sin(2*pi*(fc+0)*t)+2*sin(2*pi*(fc+0.1e9)*t);

xbaseband = carrier .* xpassband;

[b,a] = butter(6,B/(fs/2));
xbaseband_lowpass = filter(b, a, xbaseband);

subplot(2,1,1);
hold off;
plot(t, real(xpassband));
hold on;
plot(t, real(xbaseband_lowpass), 'r', 'LineWidth', 2);
plot(t, imag(xbaseband_lowpass), 'k', 'LineWidth', 2);
xlim([1e-7, 2e-7]);


MMSE- Channel estimated based - DFE Equalization

It seems that this PhD dissertation covers part of your question:
« Improved receivers for digital High Frequency communications: Iterative channel estimation, equalization, and decoding (adaptive turbo equalization)»
https://ntnuopen.ntnu.no/ntnu-xmlui/bitstream/handle/11250/249838/126139_FULLTEXT01.pdf?sequence=1&isAllowed=y

Edit:
4.2 tapped delay line model, 4.11-4.13 relates the (complex baseband) tapped delay line to physical channel multipath. (p.p 57)


What is the value range of log magnitude spectrograms?

Let us consider a standard discrete "spectrogram" or short-term Fourier transform with a window $w$ and  frames of length $N$. Let $\|x\|_r = (\sum_{n=0}^{N-1}|x_n|^r)^{1/r}$ denote the $r$-norm of $x$ ($r\ge1)$. Then for each frame of $x$, the first windowed Fourier transform is:
$$F(\omega)=\sum_{n=0}^{N-1}x_n w_n e^{-j\omega n}\,.$$
Hence, by the Rogers-Hölder's inequalities, for all $\omega$, and any couple  $(p,q)$ such that $\frac{1}{p}+\frac{1}{q}=1$:
$$|F(\omega) |\le \|w\|_p\|x\|_q\,.$$
In particular, for classical positive unit windows, $\sum_{n=0}^{N-1} w_n =1 = \|w\|_1$. Choose $p=1$ and $q=\infty$, and you get (over all frames):
$$|F(\omega) | \le \max_n | x_n|\,.$$
Not the tighest upper bound, but simple enough. And  honestly, using Rogers-Hölder was seen a bit of an overkill here, since simply:
$$\left|\sum_{n=0}^{N-1}x_n w_n e^{-j\omega n} \right| \le \sum_{n=0}^{N-1}\left |x_n w_n e^{-j\omega n} \right| \le \max_n | x_n|\sum_{n=0}^{N-1}\left | w_n  \right|\,, $$
but for "unit-energy" windows, you will get other estimates with $p=q=2$.
You can now plug it inside your "log of one plus" formula.
Similarly, you could find a tighter lower bound.


interpolate 1D signal depending on 3D parameter space

What you are trying to do is "guess" the behaviour of $d$ given $p_1,p_2,p_3$ without going through a, potentially costly, full evaluation of the expression that generates it.

Is this interpolation? No.

Inteprolation will be handy, to a certain extent, but the problem sounds more like one of classification. If you treat it as classification then the Discrete Fourier Transform will also be useful.

So, given the values of the set of $p_i$'s, is $d$ spiky, bumpy, smooth, other? "Spiky" now becomes a region of $p_i$. It doesn't have to be continuous, it doesn't have to be smooth. All you need to know is that as $p_i$ tends to certain places in this $\mathbb{R}^3$ space, $d$ is more likely to be "Spiky".

How do you decide if something is Spiky, Smooth or Bumpy? You look at its spectrum using the Discrete Fourier Transform. You use the spectrum as a feature that tells you what $d$ looks like. Specifically, sharp waveforms in the time domain give you wide spectra and vice versa. So, a spiky $d$, has a wide spectrum. A smooth slowly varying $d$ might show a few bumps at the low end of the spectrum but then die down quickly at the top end and of course flat $d$s will only have a DC component.

You will still not be able to do a full blown evaluation of $d$ across all combinations of $p_i$. First of all, some combinations of $p_i$ might be nonsensical, which is a good thing because it reduces the size of your search space. This is where interpolation comes in, by assuming that nearby values of $p_i$ lead to qualitatively similar $d$s. This is a judgement call but of course if the expression that leads to $d$ has $p_i$ involved in subtractions in some denominator, then small changes in $p_i$s can produce wild swings in the values of $d$. So, in the end, you might want to produce a "landscape" of $d$ for a quantised grid of $p_i$s.

Here is such a landscape but from another field:

The figure is reproduced from this paper

The authors take Jansen's model which models the behaviour of a population of neurons and want to see what happens when more than one populations are coupled together.

The model has a number of parameters but there is no function that tells you what parameters to plug to the model so that its output oscillates at a particular frequency. The model is capable of a wild repertoire of outputs from smooth, to fast oscillating to spiking to even just intermittent, a few spikes and then a saturated value.

So, the landscape that you see here was necessary to examine the behaviour of the coupled neuronal population when either of them was supposed to be oscillating in different ways.

(And yes, evaluation of Jansen's model for realistic waveforms (e.g. 20-30 seconds) is costly, the model is defined as a set of coupled differential equations which you have to integrate (carefully) across all the different populations you are trying to simulate. This can be very slow, especially when you try to do it at a very high level language with lots of layers between you and the hardware).

Your application is very similar to this but you are not as lucky as only being interested in oscillatory activity which has a well defined spectrum.

So, the challenge here will be how do you quantify the behaviour of $d$ using the right features to characterise the "bumpiness" (or smoothness) of $d$.

Just out of curiosity, what does the expression that lead to $d$ looks like?

Hope this helps.


Estimating the amplitude of a particular frequency

The usual way to estimate the amplitude of a particular frequency is to use the Goertzel algorithm.  There is a good write-up by Rick Lyons here.

Even though Rick's writeup is about single tone detection, it can be applied when multiple tones are present, too.


How to take the FFT in a limited range of frequencies?

If your sampling frequency is $f_s$, and you calculate the FFT with $N$ samples, then the FFT calculates the signal's spectrum at $N$ frequencies from $-f_s/2$ to $f_s/2-f_s/N$, inclusive. In Matlab notation, this would be written as:

freq_vector = (-N/2:N/2-1)./(f_s/N);

(Note that you have to shift the spectrum returned by the FFT operation to correspond to this vector; in Matlab, you would use fftshift.)

I don't know of any FFT algorithm that lets you restrict the frequency range. I would suggest using Goertzel's Algorithm to find the spectrum at the particular frequencies you're interested in. You may also use a bandpass filter to eliminate the unwanted frequencies before performing the FFT.


How to prove that this impulse response is stable?

If you check graph of h(t) which is underdamped oscillations presrnt  in 1st and 4th quadrant only. And integration of this h(t) from limits (-infinity to + infinity )  Is finite value. therefore,  system is stable


Simplest way to generate AR(2) process on MATLAB

The simplest way to approximate an AR-2 process in Matlab / Octave is the following:

N = 1024;                    % number of process samples.
a = [1, -0.9, 0.2];          % denominator coefficients, p = 2.
b = [1.0];                   % numerator coefficient.
x = filter(b,a, randn(1,N)); % generate N sample of AR-2 x[n].

Note: an AR process requires a true-white noise sequence $v[n]$ at the input of the filter but here we input a single instance of a crude approximation of it. Hence the process is not truly an AR-2 but an approximation...


The FFT size can be longer than the window size.  For instance you can use audio sample windows of 5512 (about 8 per second) and zero-pad them to length 16384 for each FFT.  If noise level is low enough and the separation between frequency peaks is great enough, the interpolation provided by this zero-padding can provide better frequency estimation and plot resolution.

The "highest resolution" FFT reasonable depends on the signal-to-noise ratio as well as the window length, as after a certain point, interpolation and frequency estimation methods will be affected as much by nearby spectral noise as by the spectral peaks that one is trying to estimate at "high resolution".


Inconsistency Between Analysis and Simulation Results of Fourier Transform of a Sinc Function

I have recently answered a related question.
What happens is that MATLAB fft is calculating a sum in the sampled points, the integral would be approximated with a (Riemann) sum.
$$ \int_{0}^{T} f(x)dx \approx \sum_{k=0}^{N-1} x(k \Delta t) \Delta t$$
The $\Delta t=0.01$ factor matches the two results.


Since the blackman window is defined in relation to its length (at least, that is, with the parameter chosen canonically), there's too few parameters open to match the band edges. You could, however:

Determine $\omega_0 \over \pi$ (the $\pi$ being the scale built in to the $\operatorname{sinc}$ function). $\omega_0$ could be picked as the arithmetic mean of the band edges.

> w0=600/7200
Window a sinc scaled accordingly by a blackman window.

 > F=sinc((-71.5:71.5)*w0)).*blackman(144)'

F now contains a sinc windowed by the blackman function, centered about the middle of the filter (giving a constant phase filter).

The resulting filter is quite steep, but it has a rather leaky stopband. If the task at hand was a different one, please excuse and explain in a little more detail.


$x[n]$ after sampling of $cos(16\pi t+\phi)$ at 12kHz

First of all, as you said the sampling rate is probably $12$ Hz, rather than $12$ kHz, and perhaps they want to demonstrate an aliasing example.

Given a bandlimited continuous-time periodic signal
$$x(t) = \cos(16\pi t + \phi)$$
the samples taken at the rate $F_s = 12$ Hz will be denoted as $x[n]$ and will be obtained by via $x[n] = x(t_n)$ with $t_n = n T_s = n/Fs$ :

$$x[n] = \cos( 16 \pi \frac{1}{12} n + \phi) = \cos(\frac{4 \pi}{3} n + \phi)$$

$x[n]$ is a discrete-time periodic sequence with frequency $w_0 = \frac{4\pi}{3}$, however, in the discrete-time case we are interested in the range of frequencies that fall in $ [-\pi,\pi)$, and therefore this signal's frequency will be cast as:

$$ x[n] = \cos(\frac{4 \pi}{3} n + \phi) = \cos(\frac{2 \pi}{3} n + \phi)$$

which indicates that the original continuous-time frequency of $8$ Hz is lost, and a new aliased frequency of $4$ Hz is attained...


Getting stuck into an open source image processing/computer vision project

I believe it would if you intend to practice your knowledge, rather than go into academia.

I suggest you look into the code base of opencv, simplecv, and scikit-image, as suggested. Once you become comfortable with using them, you can submit patches for bugs and new algorithms.

Nevertheless, you will also need to understand what you are doing at a theoretical level, so get some books too. There are also courses online, such as this and that.

Good luck!


Obtaining normalized matrix for the Haar Wavelet Transform

After a few quick calculations, it seems to me that the trouble comes from poor notations for the root in your reference. If you read, in the final normalized matrix, $\sqrt{8/64}$ and $\sqrt{2/4}$ instead of  $\sqrt{8}/64$ and $\sqrt{2}/4$ (along with the $\pm$ signs), then  the final result is correct.

The matrices $V_i$ are orthogonal. To normalize them, you only need to multiply them by a diagonal matrix $D_i$ made from the inverse of the norm of each column.

For $V_1$, the norms are all equal to $\sqrt{1/2^2+1/2^2}=\sqrt{1/2}$. Hence, you can multiply $V_1$ by the diagonal matrix
$$D_1= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2}\}\,.$$
and get an orthonormal matrix. Similarly, you get that:
$$D_2= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,1,\,1,\,1,\,1\}\,,$$
and
$$D_3= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,1,\,1,\,1,\,1,\,1,\,1\}\,.$$

Finally, because of commutativity, you get a matrix $D=D_1 D_2 D_3$:

$$D= \operatorname{Diag} \{2\sqrt{2},\,2\sqrt{2},\,2,\,2,\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2}\}\,.$$

If you now multiply the matrix $W$ (the one with the $1/8$)

with that matrix $D$, you get a final matrix $W$:

provided the notation interpretation for the square root. For instance, you can recover $1/8\times 2\sqrt{2} = \frac{\sqrt{8}}{\sqrt{8^2}} =\sqrt{8/64} $.


Just use the definition of the $\mathcal{Z}$-transform:

$$B(z)=\sum_{n=-\infty}^{\infty}b[n]z^{-n}\tag{1}$$

from which it follows that

$$B(-z)=\sum_{n=-\infty}^{\infty}b[n](-z)^{-n}=\sum_{n=-\infty}^{\infty}b[n](-1)^nz^{-n}\tag{2}$$

From $(2)$, the sequence that corresponds to $B(-z)$ is $b[n](-1)^n$. Now taking the inverse $\mathcal{Z}$-transform of

$$B(z)+B(-z)=2c\tag{3}$$

gives

$$b[n]+b[n](-1)^n=2c\delta[n]\tag{4}$$

For odd $n$, the left side of $(4)$ equals zero. For even $n$ it equals $2b[n]$. Consequently, equation $(4)$ can only be satisfied if $b[n]=0$ for even $n$ except for $n=0$, where we require $b[0]=c$:

$$b[n]=\begin{cases}c,&n=0\\
0,&n\text{ even}\end{cases}$$


QPSK, from Square Root Raised Cosine to Rectangular Impulse Shaping

You can use rectangular pulses, but it is not necessary to do so.  Rectangular pulses have advantages and disadvantages.  The advantage is that they are very simple and easy to use, both for the modulator and the demodulator.  It makes life easier for the demodulator because it makes symbol timing very forgiving- you can be off on where you sample the symbol and still get a good measurement.

There are two disadvantages to using rectangular pulses.  The first is that it uses a lot of bandwidth.  The Fourier transform of a rectangular pulse is a sinc function, which drops off very slowly.  The second disadvantage is related to the first- very few real world solutions use rectangular pulses because they use a lot of bandwidth.  If you want something realistic, in terms of what most designers would do professionally, you should not use rectangular pulses.

You could, perhaps, get it working with rectangular pulses, and then when it is solid reintroduce the root raised cosine pulses.


Frequency response of a digital state variable (Chamberlin) filter

The main reason you're seeing these results is because your transfer function is in the $z$-domain. The unique frequency response is limited to $[-\pi, \pi]$ and is periodic. In your case, you are using invalid values of $f$ to define where the $z$-domain spectrum exists. This is why you get those spikes. Using your frequency vector, I've replicated your result below

What you're generating is somewhat actually correct, but you're looking at it in the wrong way.
If you choose the range to be between $[0.1, \pi]$, you'll get the following frequency response

It's a little cleaner and you can play around with how you define your function and frequency axis to make it look more ideal. If you look at the data marker, the 3-dB point is still at that pesky $2 \times10^{-5}$ Hz. This is still correct! It's just that this is a "discrete" frequency.
In order to yield the continuous-time frequency $f$, you must convert from the discrete-domain frequency $f'$ by using the formula
$$f = f'fs  \space \text{Hz}$$
Using the 3-dB cutoff frequency we get
$$f = (2 \times 10^{-5})(500 \times 10^3 )= 10 \space \text{Hz}$$
Which is exactly what you expect. There's a related answer I have regarding continuous vs discrete frequencies along with MATLAB-specific details here.
EDIT: Using the filter
To show that the filter works, we generate and filter the signal
$$x(t) = \cos(2\pi(5)t) + \cos(2\pi(10)t) + \cos(2\pi(50)t)$$
Which are sinusoids are $5 \text{ Hz}$, $10 \text{ Hz}$, and $50 \text{ Hz}$.
I used MATLAB's filter() function which takes the transfer function's coefficients and generates a difference equation to perform the filtering. Identifying the coefficients is trivial given that your function is already in rational form. Python should have an equivalent.
Below is the signal $x(t)$ before and after filtering

You can see that the $5 \text{ Hz}$ component survives, the $10 \text{ Hz}$ component is partially attenuated since it's at the 3-dB cutoff, and the $50 \text{ Hz}$ tone is suppressed. Below is the MATLAB code I used to generate these results.
%% Sampling and constants

fc = 10;
fs = 500e3;
fn = fs/2;

q = 1/sqrt(2);
Q = 1/q;
Fc = 2*sin(pi*fc/fs);

%% Manually define the transfer function. Uncomment to generate and manually plot the frequency response.
% f = logspace(-10, pi, 1e5);
% w = 2.*pi.*f;
% z = exp(-1i.*w);
%
% freqResponse = (Fc.^2.*z)./(1 - z.*(2 - Q.*Fc - Fc.^2) + z.^2.*(1 - Q.*Fc));
%
% figure;
% semilogx(f, 20*log10(abs(freqResponse)));
% xlabel("Normalized Frequency (Hz/sample)");
% ylabel("Magnitude (dB)");
% axis tight;
% ylim([-120 10]);

%% Using built-in function filter()

b = [0 Fc.^2];
a = [1 -(2 - Q.*Fc - Fc.^2) (1 - Q.*Fc)];
[h, w] = freqz(b, a, 1e5);

figure;
semilogx(w./(2*pi), 20*log10(abs(h)));
axis tight;
ylim([-120 10]);

%% Use the filter to process a signal

t = 0:1/fs:2;

x = cos(2*pi*(5).*t) + cos(2*pi*(10).*t) + cos(2*pi*(50).*t);

nfft = 10*numel(x);
f = fs.*(-nfft/2:nfft/2-1)./nfft;

figure;
subplot(2, 1, 1);
plot(f, abs(fftshift(fft(x, nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Original Signal");

subplot(2, 1, 2);
plot(f, abs(fftshift(fft(filter(b, a, x), nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Filtered Signal");


How different methods of getting spectra in Python actually work?

You made a vector of time points spaced at 0.01 seconds. So your sampling frequency is 100 Hz.

You multiplied this time vector by $0.2\pi$ for the $\sin$ term and $3\pi$ for the $\cos$ term. So the two frequency components are 0.1 Hz and 1.5 Hz.

In your first attempt with FFT, you can get the power spectrum by calculating the square of the spectrum:

power = np.abs( np.fft.fft(y)/N )**2.0

The first half of this vector corresponds to analog frequencies from 0 to Nyquist, which in your case is 50Hz. So:

freq = np.linspace(0, 100./2, N/2)

and then plot

plt.plot( freq, power[0:N/2] )

You should see two spikes, one at 0.1 and the other at 1.5.

When using the Welch method, remember that your spectrum will get smeared and attenuated due to windowing. The right way to use that function is to provide it your sampling frequency as an input:

f, Pxx = sp.welch( y, fs=100.0 )

Since your sine wave frequency is 0.1Hz and 10 seconds of data barely fits 1 full cycle of the sine wave, Welch method will not be able to pick out that low frequency very well. If you increase your signal length to, say, 100 seconds and try the same thing with a longer window parameter, say nperseg=1024, you'll see a much cleaner double peak using sp.welch.


Why can the :$\mathbf W^*_{MMSE} $ be written as $E[\mathbf s \mathbf y^*]E[\mathbf y \mathbf y^*]^{-1}$

The mean squared error equalizer is one that minimizes the mean squred error

$$J(\mathbf{W})=\mathbb{E}\left[\|\mathbf{s}-\mathbf{W}\mathbf{y}\|^2\right]$$

which is the equalizer that minimizes the Euclidean distance between the transmitted signal $\mathbf{s}$ and the equalized received signal $\mathbf{Wy}$. To find the value of $\mathbf{W}$ the minimizes $J(\mathbf{W})$, take the derivative with respect to $\mathbf{W}$ and put it equal to $\mathbf{0}$, which will result in

$$-2\mathbb{E}\left[\left(\mathbf{s}-\mathbf{Wy}\right)\mathbf{y}^H\right]=\mathbf{0}$$
Solve the above equation with respect to $\mathbf{W}$ you get the MMSE equalizer.

Note: $(.)^H$ is the conjugate transpose, which is equivalent to $(.)^*$ in the paper.


Under what conditions do the phase margin and Nyquist criteria give the same results?

As you have pointed out, the nyquist stability criterion is more general, moreover, is the only stability criterion. Nevertheless, the Bode plot give us the exact same information that the nyquist plot does, so the stability of a system can be determined with either the bode or the nyquist plot, some times its easier with the former and other times with the latter. for minimum phase system this is usually easier with the bode plot, since you only have to see that (as you said) $$\angle L(j\omega) > -\pi/2 $$
at all frequencies for which
$$|L(j\omega)| > 1 $$

Answering your question explicitly, this is only valid for minimum phase systems. For systems that are not minimun phase the gain and phase margins are not as straightfoward to find.


Modifying sizes and thickness of letters to "embed" an image in text

...give the text an appearance of an image when looked from far away.

The technique you are looking for is called "Halftoning".

It works with all sorts of tiles as long as you can get an indication of their "grayscale value" or, in the case of coloured images, the "grayscale value" for each colour (RGB).

Consider for example these two letters in a very elementary and rather retro "dot matrix" representation. Grey values represent background, White values represent foreground.

The amount of "dots" that each letter turns on is 10 and 5 respectively. Therefore, "A" would be brighter than "I" if they were to be used as pixels.

If you count how many "dots" each character (in a particular alphabet of characters you may have available) turns on and call that "brightness", you can then order the letters from darkest to brightest and create a correspondence between pixel value and letter that has to be typed.

Similarly, you can do this for any coloured image and create a colour palette (a set of RGB values, potentially ordered by their brightness).

This is the simple version of half-toning.

There is a tonne of software that do this (or slight variants) already, such as caca and hiptext  (or even img2txt).

The next step is to apply dithering which is a way of reducing the impact of the quantisation of the image from fine grain pixels to the tiles or maintain a list of alternative tiles for the same "grey value" to avoid repetition. Some of the tools above already cover these.

Beyond this simple method there are variants that given a set of tiles they will also try to match the "shape" of the tile to the shape of the image patch they are placed over. The output of those programs is really impressive. Unfortunately, I could not find the one I had in mind but there seems to be this one and one of the papers describing the technique is this one

Hope this helps.


Are there any computer vision algorithms specifically targeting depth images?

My favourite 2.5D/3D feature descriptor for registration and recognition is the spin image (original paper + more details in Ph.D. thesis and software available from CMU).

Other recent advances (all searchable on-line for suitable algorithims) include: 3D-Sift, Fast Point Feature Histogram,  Normal Aligned Radial Features (NARF), Depth Kernel Descriptors. Older methods simply used surface properties such as curvature and edges to identify region patches.

Which is best? Depends on what you want to find, viewpoint invariance, additional clutter, etc.


OK, this is a slightly constructed situation, but as far as I can see, the following is the only thing that makes sense, kind of ...

  If $A$ is your DC voltage, then a phase of $45$ degrees means multiplying it with $e^{j\pi/4}$, i.e., you get $Ae^{j\pi/4}=\frac{A}{\sqrt{2}}(1+j)$. So scale the voltage with $1/\sqrt{2}$, and - apart from the ground wire - use two wires connected to '+' coming out of that box. Stick a tag with $j$ on it on one of the two.


Are there any realtime voice activity detection (VAD) implementations available?

So, this is not an answer in the quality I'd like it to be, but the Mumble VAD is in AudioInput.c , as a quick browsing yielded:

speex_preprocess_ctl(sppPreprocess, SPEEX_PREPROCESS_GET_PROB, &prob);
fSpeechProb = static_cast<float>(prob) / 100.0f;

// clean microphone level: peak of filtered signal attenuated by AGC gain
dPeakCleanMic = qMax(dPeakSignal - gainValue, -96.0f);
float level = (g.s.vsVAD == Settings::SignalToNoise) ? fSpeechProb : (1.0f + dPeakCleanMic / 96.0f);

bool bIsSpeech = false;

if (level > g.s.fVADmax)
    bIsSpeech = true;

What's done here is that Speex, a FOSS voice codec and utility pack is used to get an estimate for the silence->voice transition. Whyever, Mumble converts that probability value (probabilities: $\in[0,1]$) to a percentage, and then compares it to the internal setting.

So the algorithm of interest here is definitely not within Mumble, but libspeexdsp!

Funnily, that's basically where it goes down in speexdsp's preprocess.dsp

/* FIXME: This VAD is a kludge */
st->speech_prob = Pframe;
if (st->vad_enabled)
{
    if (st->speech_prob > st->speech_prob_start || (st->was_speech && st->speech_prob > st->speech_prob_continue))
    {
        st->was_speech=1;
        return 1;
    } else
    {
        st->was_speech=0;
        return 0;
    }
} else {
    return 1;
}

So, this really just boils down to an SNR estimate of the current audio frame:

/* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */
Zframe = 0;
for (i=N;i<N+M;i++)
    Zframe = ADD32(Zframe, EXTEND32(st->zeta[i]));
Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands)));


Baseline Correction: What is the concept of a baseline shift and baseline correction?

Baselines are visually simple, most of the time, but can be a pain to filter.
With your graph, for instance, what should happen between 380 and 400 is not clear.

With fluorescence, it is not sure you can get a fine model of the baseline. So effectively, in many analytical instruments, there are more or less automatic baseline fits: constant, linear, parabolic.

Without hard models, it is common to use soft ones, but we need to know what you want to quantify in your data. If it looks like other analytical chemistry data, where you want to measure peaks above a baseline, you may try the technique given in Chromatography Baseline Placement, combining low-pass filtering and peak preservation.

Is there any chance you could share the dataset?


I'm not sure if there are clear definitions for each term but here are my suggestions:

Symbol rate offset is the difference between the symbol rate entering a receiver and the symbol rate used by the demodulator. It could be measured in Hertz or parts per million (ppm).

Symbol timing offset is the difference between the time that a symbol should have been sampled and when it was sampled. It could be measured in seconds or fractions of the nominal baud duration.

Sampling clock offset is symbol rate offset multiplied by samples per symbol.

John


Calculating transmitter power from spectral density

You're right, to convert from a power density to a power, you integrate over frequency!
However:
Well, if neither of your transmitter nor receiver are calibrated, then you simply can't find the "real" transmitted or received power. The plots you have are "amplitude relative to maximum ADC input". You don't know what power "maximum ADC input" corresponds to without calibration.
End of story - there's not a single absolute power in your system that you can use as reference!
The solution here is to calibrate; i.e. send a signal, measure it with a calibrated spectrum analyzer, measure the same signal with your USRP, and then calculate the factor between these two powers. Due to linearity, that factor (==corrective additive term in dB) stays the same for all powers.


I think you could say that in simple terms this is the essence of the DFT. Just note that the DFT correlates against complex exponentials:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

You've just computed the (negative) imaginary part of $(1)$ (for $k=2$):

$$-X_I[k]=\sum_{n=0}^{N-1}x[n]\sin(2\pi nk/N)\tag{2}$$

For your example that's also sufficient because the real part is zero since your signal $x[n]$ is one full period of a sine wave. However, imagine you added a phase to your signal, then you would also need the real part of $(1)$. In the worst case you could choose a cosine for $x[n]$, which would give you a correlation of zero, even for the same frequency, if you just correlated with sine functions as in $(2)$.


Poles and zeros in time domain of analogs filters

I think you calculated the poles incorrectly, namely for $\omega$ instead of the Laplace variable $s=j\,\omega$. You can find the appropriate poles in terms of $s$ by multiplying your poles by $j$. This will give you conjugate poles, with negative real parts.


Is this system an FIR filter or an IIR filter?

It is an IIR filter.  I think you misunderstand what is meant by "finite impulse response".  It does not mean finite amplitude, it means finite time- i.e. that the impulse response eventually becomes zeros.

The impulse response of the filter is a step function, which is infinite in time.

x[-1] = 0
x[ 0] = 1
x[ 1] = 0 ...

y[-1] = 0
y[ 0] = 1
y[ 1] = 1
y[ 2] = 1 ...


Why would one use a Hann or Bartlett window?

In reviewing Fred Harris Figures of Merit for various windows (Table 1 in this link) the Hamming is compared to the Hanning (Hann) at various values of $\alpha$ and from that it is clear that the Hann would provide greater stopband rejection (The classic Hann is with $\alpha =2$ and from the table the side-lobe fall-off is -18 dB per octave). I provided the link as you can see many more considerations involved when choosing a window for various applications.
The result of this is apparent when comparing the kernels for a 51 sample Hann and Hamming window using Matlab/Octave. Note the higher first sidelobe level with Hann but significantly greater rejection overall:

Personally, I would not use either window for filter design. If any window, I would use the Kaiser window, or preferably firls. See FIR Filter Design: Window vs Parks McClellan and Least Squares for the related discussion.
I convolved a 26 sample Hann with a 26 Hamming to come up with an alternate 51 sample "Hann-Hamming" with the following result:

UPDATE: This Hann-Hamming does not (generally) out-perform a Kaiser window of similar main-lobe width:

I then tried what I call a "SuperKaiser" where I convolved two shorter length Kaiser windows to come up with an alternate 51 tap window with the following result. This was done by convolving Kaiser(26,5.5) with Kaiser(26,5.5) such that SuperKaiser(51,5.5)= conv(kaiser(26,5.5),kaiser(26,5.5). At first glance it appears to generally outperform the kaiser(51,12), matching the main lobe width and providing superior stopband rejection over most of the stopband. An integration of total stopband noise under the assumption of AWGN is of interest to see if this new window is superior under that condition (does the relative area under the first two sidelobes where SuperKaiser is inferior completely offset all the remaining stopband improvement?). If I have time I will add that assessment. Interesting! As @A Concerned Citizen astutely pointed out, as we keep convolving we should approach a Gaussian window.


what is difference between Alias vs. down sampling?

"Alias" is to "downsampling" what "a broken leg" is to "jumping from a 40-feet height". An effect from a cause, with a very probable risk in general, but could be harmless if you are lucky.

I am not sure there is a clear shared difference between uses of downsampling and undersampling. To me, used as a verb, "to downsample" may suggest a rate reduction by an integer or a fraction. While "to undersample" could be more general (at some sample rate below its Nyquist rate). Some definitions equate undersampling and bandpass sampling.

Used as an adjective as in "the signal is downsampled" or "the signal is undersampled", I understand different meaning. In the first one, I suppose a voluntary action. In the second one, this can be interpreted as an observation: the signal has been sampled at a rate lower than expected.


Morphological Gradient Versus Linear Gradient in Image Processing

I've spent few years working with specialists on mathematical morphology, and I didn't get a mathematical answer why morphological gradient was better than a classical convolution. It was more a choice of ideology.

However, the results are a little bit different and according to the problem, you can choose one over the other.

But morphological gradient is super fast, because it is usually performed with a small structuring element (radius 1 or 2) of type square or hexagon, so you can use SIMD/SSE optimization.


Why does applying RRC to the signal increase bandwidth?

I wouldn't really call it "bandwidth expansion". If you do 2 samples per symbol, from a purely theoretical point of view (i.e. using perfect infinitely long sinc pulses), you could do with exactly half of the Nyquist bandwidth of the output signal.
However, RRC is not infinitely long sincs, so you get excess bandwidth.