What is the advantage of measuring IQ data non-digitally?

IQ demodulation is normally done with the "I" component being sampled at a rate that is at or near (including below) the carrier frequency of the signal, or the center of the band of interest.  Since this is well below the sample rate required to meet the Nyquist criterion (e.g. sampling at a rate greater than 2X the highest frequency of the signal spectrum), using the "I" component samples alone at that 1X  rate won't work, as you will get aliasing (likely scrambling together both positive and negative frequency spectrum above and below the sample rate into the one channel of strictly real I.F. or baseband data samples).  This will alias (make a mess of) your spectrum before you can even try to do any sort of Hilbert transform filtering.  (Unless you have a really steep filter before sampling that does the split into just USB or LSB for you.)  Sampling at exactly 2X the carrier has a similar problem.  Using the actual sampled Q component data adds information that properly keeps the positive and negative spectrum (the spectrum both above and below your "I" component sample clock) separated in the IQ I.F. or baseband samples.

IQ sampling can be done using a sample clock (oscillator) that isn't at a frequency much higher (e.g. over 2X higher) than the frequency band of interest.  Instead, only a 90 degree phase delay is required at that much lower sample clock rate (Possibly also simplifying the required filtering as well as the sample clock generation.)


python3: speech recognition - calculating features of a WAV file

you can extract most of the well-known features such as MFCC, filter bank, PLP, Prosodic, chroma features and so on ... using:

speechpy or librosa libraries in python environment. You can extract any features from *.wav file. the list of functions is here.you can choose the best of them depending on your need.
for example, you can try for MFCC by this code. also, it is very easy by speechpy functions.
Another great tool for feature extracting is OpenSMILE by AudEERING. You can use it via command line in Windows or terminal in Linux. For more info about it, you can use this link. This tool is most cited in speech research papers.


Simply put, you need a bank of passband filters.  You feed your signal through each of the filters, and sum up the outputs from the filters.

Designing the filters is where the fun comes in.

First off, assuming this is just audio (music or the like) then there's no need of special filters.  You can use the simplest and fastest and not worry too much about phase and ripple.
You will want to use IIR filters, as they require the least amount of calculation - since you have 28 bands and presumably are working in stereo, you will want to reduce the load on the CPU else the computer won't be able to do much else.  IIR uses less CPU time than FIR filters.

There are online tools that you can use to create the needed coefficients for the filters.  Matlab or Gnu Octave can also be used.
You will need to look up algorithms for implementing IIR Filters, and how to make use of the calculated coefficients.

Specifying the bands is the next fun step.
When a filter passes audio, it doesn't simply remove everything out side of its passband.  It just isn't possible to do so.  Frequencies outside of the passband of a filter sort of slope off down towards zero amplitude.
At some frequency outside of its passband, the audio amplitude will drop by one half.  This is the -3dB point, and is used in specifying filters.
What you want to do is to arrange your passbands such that the -3dB points of adjacant filters are at the same frequency.
Say you have two filters.  One should pass frequencies from 500 to 1000Hz and the next from 1000 to 2000Hz.  You would want to make the high cutoff of the lower filter to be 1000Hz, and the low cutoff of the high filter also be 1000Hz.  Design the filters such that 1000Hz is down by half in each filter (usually, just specifying it as the cutoff is enough.)  When you now sum the outputs of the two, they add up to just the audio from 500 to 200Hz, with everything in that range at (more or less) the same amplitude as before.

Now, with your audio split into banks, you can multiply the output of each filter with some constant to dampen or emphasize that band.
This image demonstrates what I mean about the way the bands have to be layed out.


What is the coherence bandwidth of the channel?

If you have a multipath channel with impulse response

$$h(t)=\sum_{k=0}^{m-1}a_k\delta(t-\tau_k)\tag{1}$$

the corresponding frequency response is

$$H(j\omega)=\sum_{k=0}^{m-1}a_ke^{-j\omega\tau_k}=e^{-j\omega\tau_0}\sum_{k=0}^{m-1}a_ke^{-j\omega(\tau_k-\tau_0)}\tag{2}$$

Assume that the delays are ordered such that $\tau_0<\tau_1<\ldots <\tau_{m-1}$. Note from $(2)$ that the periods of the complex exponentials are given by $1/(\tau_k-\tau_0)$, and the smallest period is $1/(\tau_{m-1}-\tau_0)$, which is simply the inverse of the largest delay difference. This largest delay difference is called the channel's delay spread. Its inverse, i.e., the smallest period of the complex exponentials in the frequency domain is referred to as the coherence bandwidth.

The importance of the coherence bandwidth is the following: in a frequency range which is considerably smaller than the coherence bandwidth, the channel's frequency response can be considered flat (i.e., it doesn't change much). So if your signal's bandwidth is small compared to the channel's coherence bandwidth, it will not be distorted very much. If the signal bandwidth is larger than the coherence bandwidth, the signal will experience frequency-selective fading.


What is importance of the conjugate reciprocal and none reciprocal zeros in FIR Type 2 Filters?

I would like to provide a short answer.

Consider

$$ h[n] = \sum_{k=0}^{M} b_k \delta[n-k] ~~~~\longleftrightarrow ~~~~H(z) =  \sum_{k=0}^{M} b_k z^{-k} $$

Left: LTI FIR filter impulse response coefficients $b_k$, and right: its Z-Transform.

If the FIR filter is real; i.e., coefficients $b_k$ are real, then the Z-transform is a polynomial of real coefficients, and the roots of such polynomials are either real complex-conjugate pairs. This explains why a real FIR filter's Z-domain poles & zeros are in conjugate pairs when they are complex.

To explain why, also, they are in reciprocals, we need to consider linear phase property (without losing generality) : $h[n] = h[-n]$; then since the Z-transform of $h[-n]$ is $H(1/z)$; we have the following :

$$ h[n] = h[-n]  ~~~\implies ~~~ H(z) = H(1/z) $$

And this means that if $z_0$ is a zero (or pole) of $H(z)$; i.e., $H(z_0) = 0$ then its reciprocal $1/z_0$ is also a zero of $H(z)$.

To sum up; for a real and linear phase FIR filter, if $z_0$ is a zero of $H(z)$ then $z_0^*$, $1/z_0$ and $1/z_0^*$ are also zeros; hence zeros and poles are in pairs of four when they are complex.


Is this a correct method for fixed pattern noise correction?

The black image is the sum of a fixed pattern and dark noise (which most likely follows a normal distribution since it usually arises from current fluctuations). You want to subtract the fixed pattern, but not the dark noise - subtracting random noise from a signal simply increases the overall noise, and thus decreases signal quality.

To get a good estimate for the fixed pattern, you should capture a sizeable number of frames (say 25, though 100 will of course leave you with only half the noise), and average them. Since the dark noise is (should be) uncorrelated in time, it averages out, so that you're left with a low-noise fixed pattern that you can subtract from your future images, and that will not increase the noise in your image.

Note that the fixed pattern usually depends on the exposure time (a CCD camera, for example, may accumulate electrons during shift operations), thus you'll have to do a calibration for each exposure time. If you vary exposure times often, and if it's feasible, you can set up your experiment to capture a series of dark frames after each experiment, which means that you get a calibration for each experiment.

If you subtract a low-noise (i.e. averaged) dark frame, you will get some negative values (because the dark noise occurring during image acquisition can have negative values), but the range of your image should not increase significantly. If it does, it is a sign that you have either not averaged enough dark frames, or that the fixed pattern has changed since you're using a different exposure time.


How to find SPL of a WAV file given a calibration tone in WAV format

The calibration file and the file you want to analyze should be captured on the same equipment and the same conditions. A reference will be obtained in a free space environment, no echoes, and far field. If you have a room impulse response, you can usually do a reasonable correction but the collection equipment should be nearly the same.
Most cheap SPL meters are usually A weighted.

There are also different ways conventions on how to report SPL. In Underwater Acoustics, we usually report dB ref 1 micro pascal in a 1 Hz band, so you need to correct your levels by the ENBW of the filter (fft bin)

There's a section in Kinsler and Fry's Fundamentals of Acoustics on converting levels for different frequency bands. Noise adds as 10log and signal adds as 20 log. I recommend you look at the book.


Most signal processing starts with time-sampling a signal: you have two signals, the vibration signal, $v$, and the tachometer signal, $t$. The first thing you need to be sure of is that these are sampled synchronously in time. Do you know if that's the case?

The reason order analysis is used in the analysis of rotating machinery is that the various frequencies present in the vibration signal often depend on the rotational speed of the main shaft --- which can vary with time.

Order analysis gets around this by resampling the vibration signal synchronously with rotation angle of the main (reference) shaft.

So what you need to do is to take the time-sampled signals you have and produce a single rotation-angle-sampled signal.

To generate the rotation angle signal, you need to take the tachometer signal, $t[n]$, and find the instantaneous phase of it. The way I'd start this is to find the Hilbert transform of $t$ and form the analytic signal:
$$
a[n] = t[n] + j H\left\{ t[n] \right\}
$$
where $H$ is the Hilbert transform and then the instantaneous phase is:
$$
\theta[n] = {\rm atan2}\left( {\tt Im}\{a[n]\} , {\tt Re}\{a[n]\} \right)
$$

From there you should be able to resample both $a$ and $v$ so that you get each evenly sampled in phase.

Once $v$ is evenly sampled in phase, you can then just use an FFT to do order analysis.

  yes they are sampled synchronously in time, does the hilbert transform still need to be applied ?

If they were not sampled synchronously, then it would be almost impossible to perform order analysis at all.

  What is the purpose of the sweep filter ?

The aim of the sweep filter is to turn the tacho signal into as pure a tone as possible.  If it's not a tone, then the Hilbert transform trick will not work correctly for finding the instantaneous phase.

  do i still need to apply FFT? Or does FFT only works on the broadband signal?

Once the vibration signal is resampled, it will be synchronous with the main shaft.  To get the orders (frequencies, in effect), you need to have some sort of frequency plot. The FFT is the easiest.

The effect of this analysis is illustrated in this Matlab application note. I believe the matlab function rpmordermap effectively computes the short-time Fourier transform (but for "orders" rather than "frequencies").


How can one put information within a signal waveform and make it to traverse through different kind of medium?

Maybe what you want is information theory. Information can be quantified by Shannon entropy. An analog signal channel always has some level of noise which limits the rate of information that can be sent through it. That limit is called the Shannon limit, or more specifically for a given channel bandwidth the Shannon–Hartley theorem.


Is there better than hamming code for correcting errors

I think the question hinges around the trade off of how heavily you code the information and the rate. Generally adding more parity bits will give better error detection/correction ability but it will reduce the communication rate since each each information word is being mapped into a longer code word. Because of this trade off, I am not sure you will find better codes than the mentioned code in terms of it having less parity bits while increasing error detection/correction ability.


Checking convolution property II of z transform

The term-by-term product of the sequences $(a_0, a_1, a_2)$ and
$(b_0, b_1, b_2)$ is $(c_0, c_1, c_2) = (a_0b_0, a_1b_1, a_2b_2)$. The corresponding $z$-transforms are
\begin{align}
A(z) &= a_0 + a_1z^{-1} + a_2z^{-2}\\
B(z) &= b_0 + b_1z^{-1} + b_2z^{-2}\\
C(z) &= a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}
\end{align}
According to the formula in Wikipedia pointed out to you by MattL, $C(z)$ can be obtained from $A(z)$ and $B(z)$ via a contour integral in
the complex plane.  The integrand of the contour integral (which
is with respect to the complex variable $v$) is
\begin{align}
v^{-1}A(v)B\left(\frac zv\right)
&= v^{-1}\left(a_0 + a_1v^{-1} + a_2v^{-2}\right)
\left(b_0 + b_1\left(\frac zv\right)^{-1}
+ b_2\left(\frac zv\right)^{-2}\right)\\
&= \frac {a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}}{v}
+ \sum_{m=0}^2\sum_{n=0, n\neq m}^2
\left(a_mb_nz^{-n}\right) v^{-m + n -1}\tag{1}
\end{align}
where in that last sum in $(1)$, the exponent of $v$ in each term is either $\geq 0$ or $\leq -2$. In other words, if we regard each of the $8$ terms in that sum as a separate function of $v$, then the function
either has no poles inside the unit circle ($-m+n-1 \geq 0$) or
has a pole of multiplicity
greater than $1$ at the origin ($-m+n-1 \leq -2$).
Now, write the contour integral of the right side of
$(1)$ as the sum of $9$ different contour integrals and
apply Cauchy's integral theorem to evaluate each of these $9$ contour integrals. The conclusion is that the integrals of all the terms
in that sum in $(1)$ integrate to $0$ while the contour integral
of $\displaystyle \frac {a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}}{v}$
(which has a single pole at the origin) gives the desired
$C(z) = a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}$.
Generalization of all this to sequences of arbitrary lengths is
straightforward.

If you don't know
enough complex-variable theory to figure out why it is that the contour
integral of the first term in $(1)$ is just the numerator and why the integral of the other terms is $0$, then you
will not be able to carry out the verification that you desire.


The easiest way (to my mind) to solve the problem is to

Use the identity $\sin(A\pm B) = \sin A \cos B \pm \cos A \sin B$, substituting the known numerical values of $\cos B$ and $\sin B$,
Gathering the results to express your sum of sinusoids in the form of $C \sin A + D \cos A$,
Expressing the resulting function as $\sqrt{C^2+D^2} \sin\left(\omega t + \theta\right)$


PSD in MATLAB, problem with the coefficient $\frac{1}{2\pi}$

I would simply say that it doesn't matter so much. The Fourier Transformation and its inverse are a pair and the two formulas are entangled and require a normalization factor, which is up to conventions. See eg enter link description here

If your theory/model does not require a special convention, it is basically a degree of freedom that you can choose to fit your needs.


Are pulse shaping and windowing functions the same in concept?

Windowing and pulse shaping are very similar, and are implemented in similar ways, but their purposes are different.

Windowing means multiplying a discrete signal by a pulse with certain spectral properties before calculating its Fourier transform. The purpose is reducing the effect of discontinuities at both ends of the discrete signal.
In QAM and related digital signals, the information is transmitted in the amplitudes of a sequence of pulses. Pulse shaping means choosing a pulse with the desired spectral and orthogonality properties.

Note that some digital communications textbooks assume (many times implicitly) a rectangular pulse shape, which is then converted to another pulse shape (using a filter, as in the example in your question). In my opinion, this formulation is unncecessarily complex -- it is better to explicity select  a pulse shape from the start.

Note also that pulse shaping may be understood as the convolution of a stream of deltas with a pulse shape; in that sense, pulse shaping is a sort of "dual" of windowing.


Where can I get passport size face photo database for my image classification experiment?

Here is an up-to-date index of face image databases commonly used by researchers.

http://www.face-rec.org/databases/

It contains links to the common datasets such as The Yale Face Database (regular and "B"), The Color FERET Database, The AT&T Database of Faces, etc.

It should cover your needs.


Converting a discrete real time domain to a discrete complex time domain signal

My updated answer to OP's update

To do this specifically for the case of the FFT of a product of two sinusoids, you must choose a sampling rate or modulation frequency such that the upper sideband of the "positive" frequency component aliases to be on top of the lower sideband of the "negative" frequency component. For instance using a carrier frequency of 1.2MHz and a signal frequency of 50KHz would give you one peak in the center with two sidepeaks. (Although to note with 4095 points this wouldn't match the plot position since the center bin would be 2047 which is closer to where the lower sideband is in the plot).

Otherwise for an arbitrary signal the fft of the following would provide this result:

$$s[n] = Ae^{j\omega_c n}B\cos(\omega_m n)$$

Where

$\omega_c$ is the normalized radian frequency of the carrier ($0$ to $2\pi$)

$\omega_n$ is the normalized radian frequency of the modulation ($0$ to $2\pi$)

$A$, $B$ are arbitrary gain constants for each term (which simply scales the result by $AB$)

The opening part of the question was you are doing this to enhance your understanding of the frequency domain of AM modulation. So I am assuming you really want help understanding what AM modulation is supposed to look like:

When you multiply two sinusoids to get AM modulation, the result is DSB-SC (Double Side-band Suppressed Carrier): The carrier frequency is suppressed and there will only be two sidebands. This is what a cosine wave looks like in the FFT: two sidebands centered on the DC bin (use fftshift command if you want to center DC in the middle of the plot, then showing a positive and negative frequency axis which may help intuition), but since there is no DC in the cosine, that bin is supressed. So a cosine on its own is a DSB-SC signal with the carrier being DC, which is suppressed (we call that the baseband signal).  If you multiply this baseband signal with a exponential carrier $e^{j\omega_c n}$ as I did above, this causes the frequency to shift by $+\omega_c$   , moving what was centered at DC to now be centered at $\omega_c$. If you muiltiply the baseband signal with a sinusoid, you shift in two directions (moving the signal to the positive and negative frequency carrier) since a sinusoid contains two exponential frequencies as further explained below.

What you show in the lower plot is Large Carrier AM, and specifically such a plot could be created from adding a DC offset prior to doing the upconversion:

$$s[n] = (C + A e^{j\omega_c n})B \cos(\omega_m n)$$

Converting the real signal from your product to a complex signal can of course be done, but the result of this would be the expected two sidebands only that you now have in either the positive frequency axis or negative frequency axis. This is the Hilbert Transform that will convert your real signal that has complex conjugate symmetric positive and negative frequencies to a complex signal that only has positive or negative frequencies depending on the sign of Q.

My earlier answer regarding why the upper and lower sidebands are not equal magnitude in the FFT result

Either window your time domain data before taking the FFT or increase the length of your time domain data (and therefore your FFT). You are seeing the effects of “scalloping loss” which occurs when the actual frequency terms are not an bin centers.  Zero padding your time domain data before taking the FFT will also help to make the result look closer to what you would expect since that will interpolate more samples in the frequency domain resulting in balanced sum and difference frequency magnitudes.

Also the command fftshift() will move “DC” to the center of the axis to give you a plot with positive and negative frequencies which may be more intuitive. Keep in mind that a real signal will have positive and negative frequency components in the FFT result since each bin represents an exponential frequency term $e^{j\omega t}$, and as is clear in Euler's identity, a sinusoidal term is composed on two exponential frequency terms. For example:

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

Your time domain repesentation would be a sine wave at 900 KHz and 1.1 MHz, so the result if you increase the time domain length by adding more samples or zero-padding will approximate the expected result showing two equally weighted frequency impulses in the positive frequency axis at 900 KHz and 1.1 MHz as well as two frequency impulses in the negative frequency axis at -900 KHz and -1.1 MHz.


Minimum number of bits required for each character

the formula is

$\lceil \log_2 210 \rceil = 8$


Can you average out focal length from different camera calibration sessions?

No, it is not a good idea.  Typically, you want to limit the volume of space relative to the camera in which you want to do your measurements, and you want to use a single calibration target appropriate for that volume.  In theory, you can calibrate with any number of different patterns.  In practice, however, most implementations assume that the pattern is the same in all calibration images.

In any case, if you insist on using multiple patterns, you still want to jointly optimize the parameters, not just to average them after the fact.  And you would most likely have to implement that yourself.


The matlab codes that implemented D. Ellis's algorithm are on their website:

http://labrosa.ee.columbia.edu/projects/beattrack/


How to get a real signal from complex ifft values?

You can plot the real and imaginary components from the FFT result separately (or the magnitude and the phase results separately, as two plots).  Then you can recombine the data from the two plots into a complex spectrum required to do a reconstruction IFFT.  Both components are required, otherwise you've thrown away half the information about your signal, and can't reconstruct it (unless the original signal was purely symmetric or purely anti-symmetric to begin with).

The plot data isn't really doubled when you do two plots because the result of an FFT of strictly real inputs in conjugate symmetric.  Thus you can create the input for your reconstruction IFFT from only the first half of the data for the real plot and for the imaginary plot (by conjugate mirroring it to the other half).


Mistake or not - Fourier Series of x(2t+3)

Ok, I see now: They did it straight from the computation. The 3/2 is because they divided T by two in $e^{j3k{\frac {2\pi} T}}$.
It’s important, because it tells us that the first intuition we have that we can apply scaling to Fourier Series or Fourier Transform and then apply transformation to the Fourier Transform  is wrong. We might mistakenly think that we can scale time, apply transform, and then apply shift by 3 to the transform to move by 3 units in time. As you can see it does not produce the same result, since the shift is by 3/scaling-factor and not straight 3.


If the watermark is additive and static and, the content is smooth, and the vide long enough you could compute the watermark that minimizes $\sum ||S(F_i - W)||^2$, for the frames $F_i$ and a fixed watermark $W$. The solution $W$ is expected to converge to the watermark.
But it is always better to use the original video.


I think a picture of the QPSK constellation would be helpful to intuitively see what is occuring. Consider two possible transitions in the QPSK constellation as shown below, from "B" to "C" (180° transition) and from "A" to "C" (90° transition):

Note the difference between a 180° transition as when we move from "B" to "C" compared to a 90° transition in moving from "A" to "C". The instantaneous amplitude will be the vector magnitude at any point along the dashed lines during the transition. If the signal and has infinite bandwidth (not realistic, but if), the transition will be instant and therefore the vector magnitude will never change (will be constant envelope). However with any filtering there will be a finite time to complete the transition, and we will observe the magnitude change according to the vector length along the path taken. I placed stars where the magnitude is minimum in each case, and it is clear how in the 180° case where we pass through the origin, the magnitude approaches zero while in the 90° case the reduction in magnitude is much less.
Based on actual filtering and prior symbol history the paths will deviate somewhat from what I shown in the picture above, but the concept still holds. For example below shows a constellation of a root raised cosine signal where we can make out the same effect with the different transitions.

Minimizing the peak to average ratio is important in power sensitive applications such as Satcom and battery operated devices. This is the motivation for modified waveforms such that no trajectories go through the origin such as $\pi/4$-QPSK and waveforms that are constant envelope such as MSK and and FQPSK.


Differential Equation and Transfer Function of Multiple Summer Integrator Block Diagram

First, transform the variables into Laplace domain for dealing with algebraic rather than differential equations, which greatly simplifies the labor. And then properly re-route those two feedback branches to simplify the block diagram yet still having the same overall transfer function.
All three sub-blocks of the system are integrators whose transfer function is $H(s) = 1/s$. And your feedback paths have unity gain.
Your initial system diagram looks like this:

The transfer function is defined like: $$ H(s) = \frac{Y(s)}{U(s)} $$
In the first step, lets move the upper feedback path, which is added to the output of the first integrator, to the left adder node.
Note that the feedback signal was $Y(s)$ . When you move this feedback signal to the front of the first integrator block, you should modify the feedback signal properly so that the two configurations produce exactly the same signals. Therefore, label the new feedback signal as $Y(s) s$, (more generally $Y(s)/H_1(s)$, where $H_1(s)=1/s$). This results in the following block diagram:

Observe that the middle adder is left free, and can be removed:

Then consider the lower feedback path which samples its signal from the input of the last integrator whose output is $Y(s)$. Therefore, the lower feedback path signal is $sY(s)$. Hence we can label the lower feedback path with the signal $sY(s)$.
This way, we can (apparently!) remove the feedback lines, as we have found their resulting (effective) input signals, and obtain the following final (simplified) block diagram in which you can merge the three cascaded blocks.

Finally, this is a simple input-output relationship (with feedback). : $$ Y(s) = \frac{1}{s^3} \{ 2sY(s) + U(s)\} $$
from which you find the transfer function as :
$$Y(s) \{ 1 - \frac{2}{s^2} \} = \frac{U(s)}{s^3} $$
$$H(s) = \frac{Y(s)}{U(s)} = \frac{1}{s^3 - 2s} $$


Interpolation (sampling frequency 44.1 kHz ➔ 88.2 kHz)

Your original 44.1 kHz sampled signal has frequencies up to 22.05 kHz, so you should lowpass filter at 22.05 kHz after dilution with zeros. Your filter should have a gain of 2. Otherwise the signal amplitude drops to half because you set half of the samples to zero. Like Jim Clay says, you can combine these things; there will be plenty of opportunity for optimization.

Downsampling (sampling frequency 88.2 kHz ➔ 44.1 kHz)

If your 88.2 kHz sampled signal is the one you lowpass filtered at 22.05 kHz, you can simply decimate by throwing away every second sample. Otherwise, first lowpass filter at 22.05 kHz to prevent aliasing, and then decimate. If you combine lowpass filtering and decimation, you can optimize again.


The Baum-Welch algorithm uses the EM (Expectation Maximization) algorithm to estimate the model parameters $(T, E, \pi)$, where:

$T$: the transition probabilities
$E$: the emition probabilities
$\pi$: probability distribution on the states

Some years ago, I made the following quick-and-dirty implementation (may be fairly broken now), for the discrete case.
Hope this helps.


How to extract the Freezing Index from Accelerometer data?

I found this article is very useful. There is a function in this article that computes the average power of the signal in a specific frequency band. And since the Freezing Index is the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively), we can use this function to get the power in each frequency band and divide them.

Here is the function:

def bandpower(data, sf, band, window_sec=None, relative=False):
    """Compute the average power of the signal x in a specific frequency band.

    Parameters
    ----------
    data : 1d-array
        Input signal in the time-domain.
    sf : float
        Sampling frequency of the data.
    band : list
        Lower and upper frequencies of the band of interest.
    window_sec : float
        Length of each window in seconds.
        If None, window_sec = (1 / min(band)) * 2
    relative : boolean
        If True, return the relative power (= divided by the total power of the signal).
        If False (default), return the absolute power.

    Return
    ------
    bp : float
        Absolute or relative band power.
    """
    from scipy.signal import welch
    from scipy.integrate import simps
    band = np.asarray(band)
    low, high = band

    # Define window length
    if window_sec is not None:
        nperseg = window_sec * sf
    else:
        nperseg = (2 / low) * sf

    # Compute the modified periodogram (Welch)
    freqs, psd = welch(data, sf, nperseg=nperseg)

    # Frequency resolution
    freq_res = freqs[1] - freqs[0]

    # Find closest indices of band in frequency vector
    idx_band = np.logical_and(freqs >= low, freqs <= high)

    # Integral approximation of the spectrum using Simpson's rule.
    bp = simps(psd[idx_band], dx=freq_res)

    if relative:
        bp /= simps(psd, dx=freq_res)
    return bp

Then, I created this simple function to return the FI:

def freeze_index(data, sf, band1, band2, win_sec=None, relative=False):
    fi = bandpower(data, sf, band1, win_sec, relative) / bandpower(data, sf, band2, win_sec, relative)
    return fi

And here is how I called it in the rolling window function:

df[col + '_fi'] = df[col].rolling(win_size).apply(freeze_index, args=(fs, [3, 8], [0.5, 3], win_size,))[step_size - 1::step_size]

I hope it's the right solution, and I wish it helps.


How mic setting will affect audio classification model?

We can't tell you what your classifier does, sorry.

But yes, you have a systematic bias in your data, and your classifier will cling to whatever is the strongest discriminator if it works as hoped. I will rename your classes to make this clearer:

100 audio samples : class 'Microphone 1 subclass 1'
100 audio samples : class 'Microphone 1 subclass 2'
100 audio samples : class 'Microphone 2'

If that strongest discriminator is the microphone used in the downsampled signals (and that is not an unreasonable assumption to make), then yes, your classifier performance will depend on the microphone used.

The "properest" way would be to go out and acquire some of class 'A' and 'B' with microphone 2; but assuming that is unreasonably much work: you could try playing back all sounds and record them through both microphones, and use random batches so that the you don't only see class C with the properties of microphone 2, but also classes A and B through 2, and C through 1. Maybe that is sufficient to reduce the role of the microphone characteristics for classification.

But really, you have a class that is really identical with being recorded via microphone 2, and your choice to call that class "C" instead of "recordings done with microphone 2" is your choice alone!


Bayesian Information criterion is independent of prior - What does 'prior' mean here?

"Prior" refers to the "already known" probability (or its probability distribution function) of an event happening


Choosing the right overlap for a window function

Please refer to this paper on the optimum overlap percentage for the Blackman-Harris window, which is derived to be 66.1%. It has a lot of useful information on spectral analysis and windows.


Mathematical operations with given discrete frequency response function

Is it just element-wise multiplication of the FFT array with the FRF
  array, followed by an IFFT

No. Filtering in the frequency domain is not trivial since it's equivalent to circular convolution (not linear) and will result in time domain aliasing. Google "overlap-add" or "overlap save" for a description of adequate algorithms

  To take it one step further, what if I only had a filtered signal y(k) and my discrete FRF, could I compute the original signal x(k)?

Sometimes. In essence what you are trying to do is an "inverse filter". This is only possible if the filter itself is invertible. In other words: information that's been lost in the filtering process cannot be recovered.  In practice that means you can't have any zero (or too small) amplitude in your filter and you need decent signal to noise at all relevant frequencies


What do you mean by noise vairance and mean in an image?

Noise is not just a sudden variation in intensity; that could also be an edge of an object, a texture pattern etc. DC bias variations could also be considered as noise.

Noise is mostly an unrelated (independent of the signal) and random signal, most typically, imposed on the true signal of interest. Note that distortion (which depends on the signal) is something else.

For practical reasons, noise is mostly assumed to be white, and uncorrelated for simplicty but there are various different noise types as well, such as bandlimited (colored) noise, correlated noise etc.

For a WSS, zero mean, noise source, the variance gives its power. Stated mathematically: The variance of a WSS random process (noise) $v[n]$ is
$$ \sigma_v^2 = \mathcal{E} \{ (v[n]-\mu_x[n])^2 \} = \mathcal{E} \{ v[n]^2 \} - \mu_x^2 $$

When the process is zero mean; i.e. $\mathcal{E} \{ \mu_x[n]\} = 0$,  then
$$ \sigma_v^2 = \mathcal{E} \{ v[n]^2 \}  $$

Which is the definition of power of the WSS random process (noise) $v[n]$

Mean of ergodic noise is, its time (sample) average in the most practical sense. Or it's exptected value as ensemble average, in a more theroretical sense.

It's difficult to separate noise from true signal. A class of techniques to reduce noise, relies on statistical methods to characterize and differentiate noise and true signal, and adjust their filter characteristics according to noise variance (noise power) versus signal variance (signal power).

Optimal Wiener filters exemplify a class of linear (LTI) filters that compute minimum MSE expectations of true signals in noise, for given signall and noise statistical characteristics.


Obtaining the frequency and phase spectra from fourier transform?

They simply are the magnitude and angle of the Fourier transform.


Channel Estimation / Equalization - Estimate Channel Inverse Using White Noise Statistics Only

A Regular Random Process is is the the result of White Noise Going through a Minimum Phase LTI System.

A Non Perfectly Predictable Random Process can be defined as (See Wold Theorem / Wold Decomposition):

$$ x \left[ n \right] = \sum_{k = 0}^{\infty} {g}_{k} u \left[ n - k \right] + z \left[ n \right] $$

Where $ {\left\{ {g}_{k} \right\}}_{k = 1}^{\infty} $ is a Minimum Phase LTI System, $ u \left[ n \right] $ is White Noise and $ z \left[ n \right] $ is a Perfectly Predictable Process.

Since the questions deals with colored noise, $ z \left[ n \right] = 0 $.

Hence the process can be defined as following:

The optimal linear predictor is given by:

$$ \hat{x} \left[ n \mid n - 1 \right] = \sum_{k = 0}^{\infty} {h}_{k} x \left[ n - 1 - k \right] $$

Where the filter $ H \left( f \right) $ is given by:

$$ H \left( f \right) = 1 - \frac{1}{G \left( f \right)} $$

Hence, given the Linear Prediction Filter one could define the Whitening Filter as:

$$ W \left( f \right) = 1 - H \left( f \right) = \frac{1}{G \left( f \right)} $$

The Linear Predictor Filter can be approximated by the LMS Filter:

Which now, the Whitening Filter can be easily defined and approximated using the LMS Framework.


This symbol denotes a rectangular pulse of length $N$:
$$\mathcal{R}_N(n)=\begin{cases}1,&0\le n\le N-1\\0,&\textrm{otherwise}\end{cases}$$
I'm not sure where it is defined for the first time, but this definition is clear from the equation above Eq. $(5.24)$ on page $130$ of the 3rd edition of Digital Signal Processing Using Matlab by V.K. Ingle and J.G. Proakis.


I'm afraid your statement isn't true. This can best be seen in a suitable choice of basis, one that simplifies the cross correlation. This basis if of course the shift invariant periodic Fourier basis on your support interval.

Let's label the basis vectors $F_n$ for integer $n$. The cross correlation of two different such basis vectors vanishes, because the vectors are orthogonal and shift invariant: $\langle F_n,F_m \rangle =0$ for $m\neq n$. If we slightly generalize the cross correlation to work for complex signals as in $$\langle A,B \rangle (\tau):=\int_{-\pi}^\pi A^*(t)B(t+\tau)dt$$ and assume normalization of the Fourier basis, then $\langle F_n,F_n \rangle (\tau)=F_n(\tau)$.

With this, we can construct your $x(t)$,$y(t)$ and $z(t)$ signals like follows:
$$x(t) := c_{xy} F_{n}(t) + c_{xz} F_{k}(t)$$
$$y(t) := c_{yx} F_{n}(t) + c_{yz} F_{m}(t)$$
$$z(t) := c_{zy} F_{m}(t) + c_{zx} F_{k}(t)$$

with mutually different integer indices $k,m,n$. The complex coefficients $c_{ab}$ for $a,b\in\{x,y,z\}$ can be chosen arbitrarily and independently.

The corresponding correlation functions are:

$$\langle x,y \rangle (\tau) = c^*_{xy} c_{yx} F_n(\tau)$$
$$\langle y,z \rangle (\tau) = c^*_{yz} c_{zy} F_m(\tau)$$
$$\langle x,z \rangle (\tau) = c^*_{xz} c_{zx} F_k(\tau)$$

So each cross correlation function is a complex sinusoid with a phase factor in front. The phase factors and the sinusoids of all three correlation functions can be chosen absolutely independently. If you want to talk about maxima instead of phase factors you can take the whole setup to the reals by either adding the conjugate basis functions or just take the real part of these correlation functions. Either way, you can chose the coefficients to give you any set of maximum lags $\tau_{xy,zy,xz}$ without the requirement of the connection between them you conjectured.

So as requested let's get more explicit. The construction principle that guarantees independence of the three cross correlation functions should be relatively obvious: The orthogonality of the Fourier basis allows us to construct the signals so that the signals $x,y,z$ only overlap pairwise. In other words, the example above demonstrates that there is no transitivity for cross correlations that would allow general statements about $\langle x,z \rangle$ from $\langle x,y \rangle$ and $\langle y,z \rangle$.

The simple construction from above makes them overlap with just one basis function, which is a complex sinusoid. This is probably why the result cannot be immediately seen to apply to the original question. If you want maxima we need to use real numbers, and single real basis functions also do not have a single maximum, so we need to have more than one basis function of overlap between the signals (except for $F_1$, which already only has one maximum per period).

The simplest full example that I can come up with right now and that follows the design principle described above is:
$$x(t)=\sin(t)+\cos(4t)+\cos(5t)$$
$$y(t)=\sin(t)+\sin(2t)+\sin(3t)$$
$$z(t)=\sin(2t)+\sin(3t)+\sin(4t)+\sin(5t)$$

The base indices that overlap between the signals are $1$ for $xy$, $2$ and $3$ for $yz$, and finally $4$ and $5$ for $xz$. The latter have a phase shift introduced to realise your requirement. The overlaps group adjacent basis functions so that we have a single guaranteed maximum.

At this point it's probably easier to inspect the resulting correlation functions graphically.

First, $\langle x,y \rangle$:

$\langle y,z \rangle$:

and finally the independent third pair $\langle x,z \rangle$:

The maxima of the first two cross correlations are clearly at $\tau=0$, while the third pair correlates best at some $\tau>0$.

This example works in both discrete and continuous time, as it only uses properties of the Fourier basis that are realised in both domains.


You seem to have figured out that the system is actually given by

$$y[n]=\frac14 y[n-1]+x[n],\quad \forall n\tag{1}$$

Since you're looking for the zero-state output, you can simply compute the system's impulse response $h[n]$ by applying $x[n]=\delta[n]$ as an input signal (with $y[n-1]=0$), and compute the response to any other input signal using this impulse response. For $x[n]=\delta[n+m]$, the response is simply $h[n+m]$.

The recursion for the impulse response is

$$\begin{align}h[0]&=1\\
h[1]&=\frac14\\
h[2]&=\left(\frac14\right)^2\\
&\vdots\end{align}\\$$

I'm sure you can take it from here and derive a general expression for $h[n]$.


Transforming Non-Uniformly Sampled Image onto Rectangular Grid

You are trying to achieve something called interpolation and this can be done in roughly one line using MATLAB.

Interpolation is well explained over the web so I won't detail the functioning in this answer (Wikipedia page is very complete and redirect to many interpolation methods, here is the MATLAB documentation for 2D interpolation).

In a nutshell you parse as input both old locations (your longitude and latitude) and the data (the color corresponding to each location) as well as new locations (the rectangular and linear grid you want) and interpolation returns the new values corresponding to these locations.

Old locations could be either linear or non-linear as well as for the output locations. Moreover You can either ask for fewer or more points as the input number.

You should read carefully the documentation to address potential drawbacks of interpolation and to take account of the different methods available.

Note: if you have scattered data instead of gridded data you should use this function instead of interp2.

Here is an example with a non-linear sampled sine on 50 points.

The red plot shows a linear resampling that keeps the original number of points (same resolution) but some information is lost near $\pi/2$ as the new sampling rate is lower than the original one at this location.

The green curve illustrates a linear resampling with a greater number of points (150).This time less information is lost but more points are extrapolated between original ones (near $0.5$ and $2.5$).

Note that it is a trivial example to illustrate the possibilities given by interp1. it does not shows properly why information can be lost or extrapolated.


As you correctly noted compressed sensing, compressive sampling, sparse sampling all mean the same thing. Some authors also call it sparse sensing.
The idea behind compressed sensing is that a sparse signal can be recovered from very few linear measurements. In symbols, if $\mathbf x$ is $N\times 1$ sparse$^\ddagger$ vector, and $\mathbf A$ is an $M\times N$ matrix with $M\ll N$, and we measure $\mathbf y=\mathbf{Ax}$, then compressed sensing theory tells us$^\dagger$ that we can exactly recover $\mathbf x$ from $\mathbf y$. This is remarkable because it says that we can recover the original signal from fewer measurements.

Dictionary learning on the other hand deals with an entirely different problem of representing a bunch of data vectors in a parsimonious way. Given a set of data vectors $\{\mathbf x_1, \mathbf x_2,\ldots, \mathbf x_K\}$, we would like to find another set of vectors $\{\mathbf v_1, \mathbf v_2,\ldots, \mathbf v_L\}$ (called "atoms") such that that each data vector $\mathbf x_i$ can be represented as a linear combination of these $\mathbf v_j$'s. The set of atoms is called a dictionary. The goal here is to learn a dictionary that is much smaller than the number of data vectors$^*$ i.e. $L < K$.

Given a set of atoms in a dictionary and a vector $\mathbf y$, the goal of sparse coding is to represent $\mathbf y$ as a linear combination of as few atoms as possible.

Finally, sparse dictionary learning is a combination of dictionary learning and sparse coding. The goal here is two-fold: finding a parsimonious representation of the set of data vectors and ensuring that each data vector can be written as a linear combination of as few of the atoms as possible.

Compressed Sensing v/s Sparse Coding
Both of these techniques deal with finding a sparse representation but there are subtle differences.

Compressed sensing deals specifically with the problem of solving an underdetermined system of linear equations i.e. fewer data points than the original signal. From an unknown sparse signal $\mathbf x$ and sensing matrix $\mathbf A$, we observe the data vector $\mathbf y = \mathbf{Ax}$. $\mathbf A$ has fewer rows than columns. Compressed sensing theory deals with the following kinds of questions questions:

Under what conditions is the under-determined set of linear equations solvable and how do we solve it in a noise-robust, computationally tractable manner?
How do we design sensing matrices $\mathbf A$ for various applications?

In contrast, sparse coding does not deal with the question of designing $\mathbf A$. Moreover you aren't interested in solving under-determined system of equations --- $\mathbf A$ is allowed to have more rows than columns.$^\%$

References:

Compressive Sensing [Lecture Notes]

Dictionary Learning

Online dictionary learning for sparse coding

Footnotes:

$^\ddagger$ Sparse means the vector has very few non-zero elements.

$^\dagger$ $\mathbf A$ and $M$ need to satisfy some technical conditions.

$^*$ Unlike standard transform methods such as Fourier transform, dictionary learning is data-adaptive. When taking a Fourier transform, the basis vectors $\mathbf v_j$'s are fixed ahead of time (complex exponentials). In dictionary learning, they are learned from data.

$^\%$ This is called an over-complete dictionary.


Is it the instantaneous sound pressure in pa?

How would a piece of software know about sound pressure? What it gets is the digitized voltage from a microphone, with zero information about the physical sensitivity of that microphone.

  or the just the instantaneous amplitude of the sound signal with no unit?

Exactly! :)

  How can I use this information to calculate the SPL(no-weigting) of the signal?

Calibrate your recording setup with a sound source of known power. That's the only way short of buying an expensive calibrated measurement microphone, amplifier, digitizer. But then, there will be little need for librosa.


When I use fft Matlab function, Do I have to choose NFFT>N and power
  of 2 (Y = fft(signal,NFFT)/N;), or can I just use Y = fft(signal); ?

Not necessarily, but it is better to have it in power of 2 for faster computations.
Also NFFT>N will make your response more clearer (just for visual purpose) as the frequency resolution doesnot changes with NFFT>N. You can also look into zero padding as it is what you will achieve with NFFT>N.

As also mentioned in the document of FFT:

The execution time for fft depends on the length of the transform. The
time is fastest for powers of two and almost as fast for lengths that
have only small prime factors. The time is typically several times
slower for lengths that are prime, or which have large prime factors.

  In my case, is signal = detrend(signal); necessary before doing FFT?

It depends on whether you want to measure the frequency component of 0-10 Hz range (specifically from your description)
If DC component is of no need then you can perform signal = detrend(signal);

  In order to have 10 Hz wide bins, using 4096 Hamming Window is the way
  to do so?

NO, first of all windowing is used to decrease the side loobes which occurs due to spectral leakage. There is a effect of response broadening due to side level suppression but I think this is not what you needed. Also do remember that windowing will be applied on raw time signal before zero padding.
You have to play with number of samples to have your desired resolution

  Do my plots need further processing in order to achieve my goal?

I cannot guess this, but might be from the above answers you can estimate what has to be done in order to achieve your specific goal.

Hope that this answer might help.


Is sampling a Fourier transformed signal and fourier transforming a sampled signal the same?

If I number your equations (1), (2) and (3) starting from the first, then your Eq (1) is correct as it represents a conceptual model of the ideal sampling process: multiplying the continuous signal with an impulse train. Consequently, Eq (3) is also correct since multiplication in time domain is convolution in frequency domain (recall that the Fourier transform of an impulse train in time domain is also an impulse train in frequency domain with a different scaling factor). Convolving the two is what causes the spectrum to repeat at multiples of the sample rate.

Your Eq (2) however is not entirely correct. Fourier transforming the signal first and sampling it later implies sampling in frequency domain. Interestingly, this should create 'aliases' in time domain, i.e., repetition of the same waveform at multiples of 'sample rate'.

Why it is not completely incorrect here is that your signal is given to be continuous and periodic. A periodic signal has a sampled spectrum due to the reason described above; technically, it is said to have a Fourier series. So when you compute $F(f)$, either you should have in time domain those complex sinusoids with Fourier series coefficients, or you should have in frequency domain a frequency domain impulses with those coefficients. This might have caused confusion to whoever made the assignment.

Update: A complex spectrum can be drawn in 3D with frequency, I and Q dimensions, but the most common practice (and that reveals the most information) is to draw two separate graphs: one for magnitude response against frequency and the other phase response against the frequency.


Why does this video's audio change depending on the playing device?

The 2 stereo channel might be out of phase (near opposites, thus cancelling each other out when mixed to mono).

The human ear-brain doesn't just sum stereo left-right to the 2 ears, but instead uses any phase differences to help determine directionality.


Finite and Infinite support in time and frequency domain

The sinc is just fine an example as a signal with infinite support:

Support being defined as the smallest interval in which the function has non-zero values, it's trivial to see that $\text{sinc}(x)=\sin(x)/x$ (for all $x$ but a single point) becomes arbitrarily small for large $x$, but never actually constantly $\text{sinc}(x>\xi)=0$ for any finite $\xi$. So, q.e.d., sinc has infinite support.

We use this exact property to illustrate to students all over the world that the perfect filter (i.e. something that looks like a rectangle in frequency domain) can never be implemented as a FIR filter, as that would need to be infinitely long.


Your first method assumes that the system is entirely characterized by its response to the input $x(t)=\delta(t)$. This would only be the case if the system were linear and time-invariant (LTI). In that case, the input-output relation would have the form

$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{1}$$

where $h(t)$ is the system's impulse response. However, the given input-output relation is not of the form $(1)$, and hence the system is not LTI. Consequently, the system's response to $x(t)=\delta(t)$ does not say much about the system's behavior in general. That's why your first method of determining the stability is not valid.

The second method is valid. You found a bounded input signal for which the output is unbounded, and in this way you've shown that the system is not BIBO-stable.


Well, there is sampling of singing voices.  it's not the latest.  one thing that you perhaps need to worry about is that, unless you get a sample root for each note going up the MIDI note scale, then you will likely be pitch shifting the notes up or down a little.  it turns out that this causes a noticeable distortion sometimes called "munchkinization" or the "Alvin the Chipmunk" effect.  to avoid that, each cycle of the quasi-periodic note needs to be isolated into little "wavelets" (not to be confused by the wavelet transform) or "grains" or "FOFs".  then you would output the grains or FOFs at a rate of the fundamental frequency of the note without stretching or scrunching the grain (this requires overlapping and adding).  there's an old paper by Keith Lent that discusses pitch shifting vocal notes without shifting formants, and another paper by me that analyzes the algorithm a little more rigorously.  the latter, i can send a copy.

but for synthesis, you need not perform the real-time pitch detection that you would have to for pitch-shifting the voice live.  what you can do is sample some notes from these sung voices and perform analysis on them off-line and generate a file of these grains or FOFs and store them in your synthesizer.


Consider the signal $$s[n]=\frac{a^{|n|}}{1-a^2},\, |a|<1$$

Assuming $$x[n]=a^n u[n]\stackrel{\text{DTFT}}\longleftrightarrow\frac{1}{1-a e^{-j\omega}},$$
 one can show that

$$s[n]=x[n]*x[-n]$$

The proof is simple. For example, for $n>0$ we have
$$\sum_{k=n}^{\infty}x[k]x[k-n]=\sum_{k=n}^{\infty}a^ka^{k-n}=\sum_{k=n}^{\infty}a^{2k-n}=a^n(1+a^2+a^4+\cdots)=\frac{a^n}{1-a^2}$$

Hence, the DTFT of the signal $s[n]$ is given by
$$S(e^{j\omega})=X(e^{j\omega})X(e^{-j\omega})=\frac{1}{1-a e^{-j\omega}}\frac{1}{1-a e^{j\omega}}=\frac{1}{1-2a\cos(\omega)+a^2}$$
In summary,
$$\boxed{\frac{a^{|n|}}{1-a^2}\stackrel{\text{DTFT}}\longleftrightarrow\frac{1}{1-2a\cos(\omega)+a^2},\, |a|<1}$$

Now consider

$$G(f)=\frac{1}{1+1.4\cos(2\pi f)+0.49}=\frac{1}{(1+0.7e^{-j2\pi f})(1+0.7e^{j2\pi f})}$$

whose IDTFT can be calculated using the above formula for $a=-0.7$ as

$$g[n]=\frac{(-0.7)^{|n|}}{1-(-0.7)^2}$$

Since
$$\frac{1+0.8(e^{-j 2\pi f}+e^{j 2\pi f})+0.64}{1+1.4\cos(2\pi f)+0.49}=G(f)+0.8e^{-j 2\pi f}G(f)+0.8e^{j 2\pi f}G(f)+0.64G(f)$$
The IDTFT is
$$g[n]+0.8g[n-1]+0.8g[n+1]+0.64g[n]$$
But obviously, the variable $\tau$ is used instead of $n$ in the give equation.


I don’t think that you can do multi-frame super resolution in the traditional sense if there is no camera rel scene movement, nor if the camera applies «proper» Nyquistian spatial prefiltering?

I think of good old interlacing as an example of potential super resolution. You get a 2-frame cadence consisting of separate time and spatial samples, usually moderately filtered. If there is no movement, you can get up to double spatial resolution. If there is movement that can not be resolved, you get high temporal resolution.


The problem was I've masked the shifted spectrum. After I've used ifftshift on my mask, I obtained the desired effect.

PS. geometrikal - thanks for putting me on the right track!


How to figure out number of multiplications in rational sampling

That depends a fair bit on how exactly you implement the conversion but if you follow the textbook approach it would probably look like this:

The generic order of sampling rate change is "upsample", "filter", "downsample".
That means that you need to run the filter at the upsampled rate, i.e 5000
Your IIR filter has 9 non-trivial coefficients, i.e. you need to do 9 multiplies per input/output sample
In total that's 9*5000 = 45000 multiplications per second.

That's a lot and that's whey rational sample rate converters often implemented as polyphase FIR filters. Using an IIR requires you to calculate all output samples of the filter although you will be throwing 6 out of 7 away.


Fusing motion and color in particle filter

If you can make the assumption that the two observations are independent then the the likelihood function you need is just the product of the two likelihood functions.

If for particle x, the likelihood of it in terms of the colour distribution is p(z_c|x) and the likelihood of it in terms of the motion contour is p(z_m|x) then the likeliood of both observations combined is

 p(z|x)
=p(z_c,z_m|x)
=p(z_m|z_c,x)p(z_c|x)
=p(z_c|z_m,x)p(z_m|x)

With the assumption of independence this becomes:

=p(z_c|x)p(z_m|x)

So just calculate the likelihoods as you have been doing and multiply them.

If you would rather not make this assumption then you need to specify p(z_m|z_c,x) or p(z_c|z_m,x).


How to find position, scale and rotation of a known shape/contour in an image?

A Hough transform can certainly work. However, with position, scale, and orientation, you would have to use a 5-D Hough transform (2 DOF for position, 1 for scale, 1 for orientation), which might be computationally quite demanding.

As an alternative, you might want to look at SIFT or one of its extensions like SURF. It sounds like they can be adapted to do the job. And OpenCV has them.


What is the relationship between passivity of impulse response and passivity of the system?

Let me clarify some definitions first. There is no such thing as a "passive impulse response", there are only passive systems. If you like, you can call the impulse response of a passive LTI system "passive", but that's not common usage. From what I understand from your question, you probably mean by a "passive impulse response" an impulse response $h[n]$ that satisfies

$$\sum_{n=-\infty}^{\infty}|h[n]|^2\le 1\tag{1}$$

i.e., an impulse response with energy not greater than $1$. However, as you've shown yourself in an example, this is not sufficient for the corresponding system to be passive. Let's use the following definition of passivity [1]:

  A discrete-time system is said to be passive if, for every finite energy input sequence $x[n]$, the output sequence $y[n]$ has, at most, the same energy, i.e.,
  $$\sum_{n=-\infty}^{\infty}|y[n]|^2\le\sum_{n=-\infty}^{\infty}|x[n]|^2\tag{2}$$

It is more straightforward to investigate the passivity of a linear and time-invariant discrete-time system in the frequency domain than in the time domain. Let $H(e^{j\omega})$ be the frequency response corresponding to the impulse response $h[n]$:

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{3}$$

Since the Fourier transforms of the input and output sequences are related by

$$Y(e^{j\omega})=H(e^{j\omega})X(e^{j\omega})\tag{4}$$

the condition

$$|Y(e^{j\omega})|^2\le |X(e^{j\omega})|^2\quad\text{for all }\omega\tag{5}$$

holds if and only if $H(e^{j\omega})$ satisfies

$$|H(e^{j\omega})|\le 1\quad\text{for all }\omega\tag{6}$$

Integrating $(5)$ gives

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}|Y(e^{j\omega})|^2\le \frac{1}{2\pi}\int_{-\pi}^{\pi}|X(e^{j\omega})|^2d\omega\tag{7}$$

from which it follows via Parseval's theorem that

$$\sum_{n=-\infty}^{\infty}|y[n]|^2\le\sum_{n=-\infty}^{\infty}|x[n]|^2\tag{8}$$

must hold. So it is not condition $(1)$ which makes a system passive, but condition $(6)$.

A causal and stable transfer function with real coefficients satisfying $(6)$ is called a bounded real transfer function [2]. This type of transfer function is important for the realization of low-sensitivity digital filters.

[1] S.K. Mitra, Digital Signal Processing, A Computer-Based Approach, 2nd ed., p. 70

[2] P.P. Vaidyanathan and Mitra, S.K., Low Passband Sensitivity Digital Filters: A Generalized Viewpoint and Synthesis Procedures, Proc. IEEE, vol. 72, no. 4, April 1984, pp. 404-423.


Shannon capacity with distribution different from Gaussian

Wikipedia gives the more general definition. Channel capacity is defined using Information Theory. If you are unfamiliar with that, it will be hard to explain. The channel capacity is always the same for a fixed communications channel (meaning the statistics are fixed).

The more general definition of the capacity $C$ is
$$ C = \sup_{p_X(x)} \mathbf{I}(\mathbf{X};\mathbf{Y}) $$ where $\mathbf{I}$ denotes the Mutual Information, $\mathbf{X}$ is the input to the channel and $\mathbf{Y}$ is the output of the channel, and the supremum is taken over all possible input distributions $p_X(x)$.

The equation that you have provided is specific to the AWGN channel. The plots are demonstrating the performance ceiling that exists because you have chosen not to use a Gaussian input sequence (for convenience usually). Notice that it follows the capacity curve to a certain limit but is then limited by the constellation size.

Codewords are relevant when trying to reach capacity, but not generally applicable in the statement of what capacity is. In other words, the channel capacity is independent of the codewords chosen.


First of all, SURF only uses it as a blob detector to find interest points.
The first order derivatives only give you directional information such as edge-ness. However, the Hessian (as the matrix of second derivatives) contains local structure information such as curvature / concavity etc. For instance, its eigenvectors can be used for cornerness measures. It is just more informative. Also, its determinant is used as an automatic scale selector.

Another common alternative would be a Laplacian detector, but Hessian shows better scale selection properties.


Beginner's book in signal processing with practical examples on fault detection in electrical motors

When I had a look at rotating machinery, the best reference I could find is Bob Randall's Frequency Analysis.  This was generated in conjunction with Brüel & Kjær as they sold lots of nice (and expensive) vibration measuring / monitoring equipment.  A short Google search shows that B&K have made the PDF of this book available for free download from their website.

The book is targeted at mechanical engineers, who (back in the day) were never exposed to the Fourier transform.


The $\omega$ in the frequency response of a discrete-time system $H(e^{j\omega})$ is indeed unitless. The frequency response $H(e^{j\omega})$ is periodic with period $2\pi$. If the discrete signal is obtained via sampling with sampling rate $f_s$, then the relation between $\omega$, the actual frequency $f$, and the sampling frequency $f_s$ is
$$\omega=2\pi\frac{f}{f_s}\tag{1}$$
Consequently, as a function of $f$, $H(e^{j2\pi f/f_s})$ is periodic with the sampling frequency $f_s$.


Does the Nyquist frequency of the Cochlear nerve impose the fundamental limit on human hearing?

No.
A quick run-through the human auditory system:

The outer ear (pinnae, ear canal), spatially "encodes" the sound direction of incidence and funnel the sound pressure towards the
ear drum, which converts sound into physical motions, i.e. mechanical energy
The middle ear (ossicles) is a mechanical transformer (with some protective limiting built-in) that impedance matches the air-loaded ear drum to the liquid-loaded oval window of the
Cochlea (inner ear). The vibration excites a bending wave on the basilar membrane. The membrane is highly resonant and transcodes frequency into location: for any given frequency the location of the resonance peak is in a different spot. High frequencies wiggle very close to the oval window, low frequencies towards the end of it. This motion is picked up by the
Cochlea neurons, which transmit the intensity of the excitation at their location to the brain. About 20% of the neurons are efferent (come out of the brain) and are used to actively tune the resonance with a feedback loop (which causes tinnitus if misadjusted)

So in essence the Basilar Membrane performs sort of a mechanical Fourier transform. The frequency selectivity of the Neurons is NOT determined by the firing pattern but simply by their location. A neuron at the beginning of the basilar membrane is sensitive to high frequencies and a neuron at the end detects low frequencies. But they are more or less the same type of Neurons.
The Nyquist criteria doesn't come into play at all since no neuron is trying to pick up the original time domain waveform. The couldn't anyway: human neurons have a maximum firing rate of less than 1000 Hz and average firing rates are way below that. The firing rate of a cochlea neuron represents "Intensity at a certain frequency" where that frequency is determine by the location of that specific neuron.
So you can think of it as a short term Fourier Transform. Instead of a single time domain signal you get a parallel stream of frequency domain signals where each individual signal has a much lower bandwidth.
A cochlea implant basically does the short term Fourier transform internally and then connects the output for each frequency range to the "matching" neurons in the cochlea nerve. Theoretically you can create ">20 kHz" hearing with an implant that can actually receive and process higher frequencies and simply routes them to existing neurons, i.e. you could feed 40 kHz activity to the 10 kHz Neuron. The human would have a sensation when exposed to 40 kHz but it's unclear what they could do with that: they would have "relearn" how to hear. Aside from the highly questionable practical and ethical issues, it probably wouldn't be useful. In order to get to 40 kHz you'd have to give some other frequencies, and chances are that evolution has chosen the current "normal" range for humans pretty carefully.


Frequency response of Low pass filter and High pass filter

In previous sections of the book, the fact that a discrete-time signal's spectrum is periodic may have been mentioned. It can be described formally as follows:
$$X(e^{j\omega})=\frac1{T} \sum_{k=-\infty}^{\infty}X_C\biggr(j\biggr(\frac{\omega}{T}-\frac{2\pi k}{T}\biggr)\biggr)$$
being $X_C(j\omega)$ the Fourier Transform of the continuous signal, and $T=1/f_s$ the sampling period. That expression is actually saying that the discrete-time spectrum can be computed from the continuous spectrum just by expanding its width $T$ times, and repeating this every $\frac{2\pi}{T}$. To keep things simple, we can normalize the frequency so that every signal's spectrum is periodic every $2\pi$, no matter how we sampled it.

Keeping this in mind, a low pass filter isn't just a "step" around $\omega=0$, but also around $\omega=2\pi,4\pi,...$ and, of course, their negative counterparts. I mean low frequencies repeat around $\omega=0, 2\pi,4\pi,...$, while high frequencies do around $\omega=-\pi,\pi,3\pi,...$

Shifting its frequency response by $\pi$, we get $H(e^{j(\omega-\pi)})$, which would look like the following:

Remember that $\omega=\pi$ corresponds to your Nyquist frequency (= half the sampling frequency).

Image credit: Oppenheim & Schafer, Discrete-time signal processing


First question is your DAC AC-coupled? In that case, the AC coupling will remove the DC offset caused by your DDS.

Second question, unless your DAC is perfect, your DAC will also have a DC offset. Do you know what's the typical DC offset level of your DAC? Again, if the DC offset caused by your DDS is less than say 10% of the intrinsic DC offset of your DAC, I don't think it's worth it to fix it.

Third question, is the ratio between the sampling frequency and the signal frequency a whole number? If it's an even number, you should not have any DC offset, but if you have an odd number of samples per period, you could get a DC offset, For example, with a ratio of 5 because you might get two round-ups and three round-downs giving your a negative DC offset.

Here's what you could try:

change the signal frequency so that the ratio between fs and f is not a whole number.
Try to add phase dithering, this could remove your DC "spur"


Impulse response of a system from the input-output relation

It looks like you answered your own question but to be complete..The impulse response $H(z)=\frac{Y(z)}{X(z)}$ can be found from the difference equation which you have found from the block diagram. So taking the z-transform of both sides of the difference equation we have: $Y(z)=X(z)(\frac{2}{16}+\frac{3}{16}z^{-1}+\frac{6}{16}z^{-2}+\frac{3}{16}z^{-3}+\frac{2}{16}z^{-4})$ and so we have $H(z)=\frac{2}{16}+\frac{3}{16}z^{-1}+\frac{6}{16}z^{-2}+\frac{3}{16}z^{-3}+\frac{2}{16}z^{-4}$.

Now taking the inverse z-transform you get the impulse response to be $h[n]=\frac{2}{16}+\frac{3}{16}\delta[n-1]+\frac{6}{16}\delta[n-2]+\frac{3}{16}\delta[n-3]+\frac{2}{16}\delta[n-4]$, which has the coefficients that you listed.


What is the advantage of MATLAB's filtfilt

You can best look at it in the frequency domain. If $x[n]$ is the input sequence and $h[n]$ is the filter's impulse response, then the result of the first filter pass is

$$X(e^{j\omega})H(e^{j\omega})$$

with $X(e^{j\omega})$ and $H(e^{j\omega})$ the Fourier transforms of $x[n]$ and $h[n]$, respectively. Time reversal corresponds to replacing $\omega$ by $-\omega$ in the frequency domain, so after time-reversal we get

$$X(e^{-j\omega})H(e^{-j\omega})$$

The second filter pass corresponds to another multiplication with $H(e^{j\omega})$:

$$X(e^{-j\omega})H(e^{j\omega})H(e^{-j\omega})$$

which after time-reversal finally gives for the spectrum of the output signal

$$Y(e^{j\omega})=X(e^{j\omega})H(e^{j\omega})H(e^{-j\omega})=
X(e^{j\omega})|H(e^{j\omega})|^2\tag{1}$$

because for real-valued filter coefficients we have $H(e^{-j\omega})=H^{*}(e^{j\omega})$. Equation (1) shows that the output spectrum is obtained by filtering with a filter with frequency response $|H(e^{j\omega})|^2$, which is purely real-valued, i.e. its phase is zero and consequently there are no phase distortions.

This is the theory. In real-time processing there is of course quite a large delay because time-reversal only works if you allow a latency corresponding to the length of the input block. But this does not change the fact that there are no phase distortions, it's just an additional delay of the output data. For FIR filtering, this approach is not especially useful because you might as well define a new filter $\hat{h}[n]=h[n]*h[-n]$ and get the same result with ordinary filtering. It is more interesting to use this method with IIR filters, because they cannot have zero-phase (or linear phase, i.e. a pure delay).

In sum:

if you have or need an IIR filter and you want zero phase distortion, AND processing delay is no problem then this method is useful
if processing delay is an issue you shouldn't use it
if you have an FIR filter, you can easily compute a new FIR filter response which is equivalent to using this method. Note that with FIR filters an exactly linear phase can always be realized.


The convolution process works by multiplying every single pixel on the image with the matrix/kernel. And for the outer pixel in top, bottom, right and left there are several ways to do it:

You can ignore it, so you can start from the second line if the kernel is 3x3.
You can put a 0 value as the outer border.
You can copy the closest neighbor pixel value in as the outer border.
You can pretend the image is circular by using the opposite pixel. Top with bottom and left with right border.

As far as I know the matrix goes from left to right direction. And for that matrix (the edge enhancement), it gets pixel that have huge difference with the left neighbor. The pixel is multiplied with 1 and the left neighbor is multiplied with -1. If the two pixels have huge difference, the result will be near to white (1 or 255) but if the two pixels have small difference, the result will be near black (0).

I hope this answers your questions.


There really is no practical difference between the two. The important thing to understand about aliasing is not the exact definition of the word, but rather the concept: when two frequency bands alias and there is a signal in one of these bands, after sampling you can't tell which band the signal came from. In a sense, by the time you're asking what is or isn't aliasing, it's already too late.
Consider this: you have a digital audio recording of perfect silence. However, somewhere during the recording process some RF-interference was mixed in. So the recording of silence now contains some aliased RF-noise, which isn't added to any existing spectrum, because there is none. After all the recording is silence, even though clearly the speaker cone is moving. Maybe the lack of spectrum was actually a feature of the spectrum?
If you understand the concept of aliasing, then it is almost always clear from the context what the exact definition in that particular context happens to be.
Edit:
I want to clarify further that what the text is really trying to say is simply if you violate the sampling condition, then the bands will alias and if the bands alias, then the signals can alias as well (but when they do, you really can't tell; at that point it's too late).


Isolating frequency-specific oscillation from transient

The notation and terminology in this and the closely related paper could have been more simplified I think.

To understand what's going on there, we need to start from page 1583. Specifically:

$y(t)$ is basically your measured signal (the output from each neuronal population). A scaled version of that is what you would sense via some modality (in an ideal world).
$s(\omega, t)$ is the frequency-time representation of a $y(t)$. This is complex.
Given that this is about Event Related Responses (ERR):

$g(\omega, t)_T$ is the average absolute value of $s(\omega,t)$ after a number of trials. This is real valued.
$g(\omega, t)_e$ is the average of $s$ times the average of the conjugate of $s$ (both of which will be complex, the result of which however will be real).
These look similar but their magnitude is not.

Then the explanation of the breakdown of that "power" follows. What this means is that, based on the assumptions in the paper, we think that the neuronal activity is due to three things: a) Background activity, b) Activity because of connections between neuronal populations and c) Induced activity by the stimulus. Furthermore, these three components are assumed orthogonal. This means that if one "happens", the rest do not. This is important for later on.
WHERE/WHEN this holds then we could estimate the background activity by looking at the spectrum BEFORE the stimulus onset, we could estimate the activity due to the stimulus (the part that is phase locked to the stimulus), because we know the "shape" of the stimulus and therefore what is left should be due to structural connections.

These are captured in $g(\omega, t))_i$ and $g(\omega,t)_T$.
Notice here that the baseline is only $g(\omega)$

So, when you arrive at equation 23, what you see is the explanation of that "adjustment" where:

$g(\omega,t)_T$ is the frequency-time representation of the averaged stimulus response after $T$ trials as before.
$g(\omega,t)_e$ is the frequency-time representation "power of average" as before
$\hat{\eta}$ now denotes the proportion of the frequency-time representation of $g(\omega,t)$ that is due to $g(\omega,t)_e$. The "trick" here is this, you have already assumed that the three components are orthogonal. WHERE/WHEN this holds, the "pseudoinverse" ($g(\omega, t)^+$) of a given $g(\omega, t)$ multiplied by the "shape" of the stimulus response $g(\omega,t)$ over the trials would be zero and you would be left with the "shape" of the stimulus response. So, that $\hat{\eta}$ business there (and the subtraction) is what allows you to do that "separation".

Perhaps it is now clear how to produce that $g(\omega,t)^+$. Provided that you have a given $g(\omega,t)_e$, invert it with pinv (or anything else equivalent to a pinv in another platform).

But, I wonder how effective would that be for what you are trying to do. Irrespectively of this specific technique, a "pulse" or a step or any abrupt change of a signal in the time domain usually has a broad representation in the frequency domain (whether that is Fourier or Wavelet, the point is that you need more components to synthesize the step). You are looking for "oscillations" in a 0-200ms interval. What is the frequency of the oscillations you are looking for?

Hope this helps


Using a Butterworth high-pass filter together with Fast Fourier transform

As the comments say: it really depends on your signals. If the energy in the "don't care" region is significantly larger than the energy at the frequency of interest, it's helpful to pre-filter it first. Otherwise you may see "spectral leakage" of the high energy components into your target area.

On the downside: the pre-filter will affect the group delay the transient response of your overall system.

This being said: if the filter is not too steep and the cutoff frequency not too close to your target frequencies, pre-filtering is generally a good idea. The pre-filtering is much cheaper than the FFT, so it have almost no impact on CPU or memory requirements.


Understanding why (at least) 2 points of the FFT of a real signal are real

First, "why (at least) 2 points of the FFT is real" is incorrect.  It is why one or two points of the FFT is real.
Each point of a discrete Fourier transform (which the FFT is just an algorithm to compute quickly) is the complex coefficient for its corresponding frequency.  In other words, the DFT starts with the assertion that the signal in question can be fully described by
$$x_k = \sum_N a_k e^{\frac{2\pi}{N}k} \tag 1$$
For $x_k \in \mathbb R$, (1) just carries the constraint that $a_k = (a_{N - k})^*$, where $x^*$ is complex conjugation.
From this, we can show that $a_0$ (and $a_{N/2}$ for $N$ even) are constrained to be real:
Start with the rule I mentioned that in order for $x_k$ to be real, $a_k = (a_{N - k})^*$.  Another property of the DFT is that it's only uniquely defined over an interval $N$ points long, and then it repeats.  So $a_{k + N} = a_k$.  This means that for $k=0$, $a_{N - k} = a_{N} = a_0$.  Applying the conjugate property, this means that $a_0 = a_0^*$ -- so $a_0$ must be real.
For $N$ even, $a_{N/2} = a_{N - N/2}$.  So, again, $a_{N/2} = a_{N/2}^*$, so $a_{N/2}$ must be real.
So -- one (or two, in the case of an even number of points) arguments must be real.


What should be the output of the pdf for multivariate normal distrubition?

The PDF $p$ of a, say, bivariate Normal distribution with mean $\mu$ and covariance $S$ is a function that takes two parameters and has a scalar-valued output:

$$ p: \mathbb{R}^2 \mapsto \mathbb{R}$$

and is given by

$$p(x_1,x_2) = \frac{1}{2\pi |S|}\exp\left(-\frac{1}{2}(\vec{x}-\vec{\mu})^TS^{-1}(\vec{x}-\vec{\mu})\right)$$

where $\vec{\mu}\in\mathbb{R}^2$ is the mean of the distribution and $S\in\mathbb{R}^{2\times 2}$ is its covariance matrix. $\vec{x}=[x_1, x_2]^T$ is the function input and $p(x_1, x_2)$ is the function output.

So, the output is a scalar (a single value). Furthermore, the PDF can become larger than one for some input $\vec{x}$.


How to compute a sweeping cutoff filter in python?

To create a smooth sweeping filter or envelope filter, you can process as in video by "frame". You plan a 30% overlap of each frame with a crossfade.
You need to generate two tracks and mix them together to do that.


Why use cross-spectral density to calculate frequency response?

in the audio biz, we call this the "Two-channel FFT".  the cool thing about it is that you can measure the magnitude response of a room or something using music (that is decently broadbanded) as the test signal (and divide the output spectrum by the input spectrum).  you don't need to pollute peoples' ears with an impulse train or a maximum-length-sequence or pink noise or with a swept-frequency chirp.


Interleaving is not specific to OFDM. It is an age old technique used for improving error correcting codes.

Basically, most error correction coding techniques are good at suppressing bit-wise errors which are independent and random. However, when channels are not memoryless, there could be burst errors. i.e. several bit errors are likely to see in close by given the first bit error occurs. Hence, in such cases typical ECCs gets much weaker.

Burst errors in general are more adverse not only because they get better of error correction codes but they tend to produce more profound local gaps (such as hicups in voice or patch in video data). And yet, in most wireless and satellite channels they are more common. This is also applicable in storage mediums.

Interleaving, a very simple technique resolves this very efficient way. Basically you take a block of bits and interleave it in a particular known order. Now when a burst error occurs, after the de-interleaving at receiver the effected bits are now speared apart much more than the original burst. And hence it improves the work of ECC significantly.

Many modern channel encoding uses interleaving built-in.

You can read a very detailed outline here:

Yun Q. Shi, Interleaving for Combating Bursts of Errors IEEE CIRCUITS AND SYSTEMS MAGAZINE FIRST QUARTER 2004 pp 29


Use wavelet transform to extract a waveform of certain frequency

The mother wavelet is a bandpass filter. This means we have to apply the CWT(in this case), calculate the correspondence between scales and frequencies, zero out everything outside the known frequency band for the alpha wave and invert the transform.

data = Import["your_sample_eeg_data"];

ListLinePlot[data]

cwd = ContinuousWaveletTransform[data, MorletWavelet[], SampleRate -> 200, Padding -> 0]

freq = (cwd["SampleRate"]/(#1 cwd["Wavelet"]["FourierFactor"])) & /@
   (Thread[{Range[cwd["Octaves"]], 1}] /. cwd["Scales"]);

ticks = Transpose[{Range[Length[freq]], freq}];

WaveletScalogram[cwd, Frame -> True, FrameTicks -> {{ticks, Automatic}, Automatic},
 FrameLabel -> {"Time", "Frequency(Hz)"}, ColorFunction -> "SunsetColors", ImageSize -> 500]

So, we are interested in the 5th octave(the region between 10.51 and 5.25, of course, roughly speaking)

mwd = WaveletMapIndexed[#1 0 &, cwd, Except[{5, _}]]

WaveletScalogram[mwd, ColorFunction -> "SunsetColors"]

It's like performing a brain surgery with a chainsaw ;__;

ListLinePlot[InverseContinuousWaveletTransform[mwd]]

Those are the basic steps I would perform if I have to extract a signal within a certain frequency band.


Calculating the total energy of a signal

The energy of a discrete-time signal is the sum of the squares:

$$
E_x = \sum_n |x[n]|^2
$$

If you need the average power of your signal, divide the energy by the number of samples in your vector. Power is the average energy per sample of your signal.

If your signal vector is a sample of a stationary process (i.e. a recording clip of someone talking), then average power probably has more meaning than energy (which will grow with the sample length).


Samples per symbol and number of symbols for QAM

Number of symbols

A baseband signal model that is used is: $x(t)=\sum_{k=1}^N a_k p(t-kT)$, where $a_k$ are the symbols (there are many types, PSK, QAM, PAM for example), $p(t)$ is the pulse shaping filter, $T$ is the symbol period, and $N$ is the number of symbols. The pulse shaping filter limits the bandwidth of $x(t)$.

  Samples per symbol

Having multiple samples per symbol is called oversampling. The samples per symbol, $\text{sps}$, is usually not explicitly chosen but is determined from the symbol rate of the signal, $F$, and sample rate, $F_s$, on a particular platform: $\text{sps} = \frac{F_s}{F}$. Your statement that "samples per symbol is sampling rate" is not correct. Visually, having more samples per symbol "fills in" the symbol and will make it look smoother. There are many reasons, besides the pulse shaping filter, to oversample. For example, some algorithms in the receiver are designed to operate on the oversampled received signal, while others expect to work on the received symbols.

Code: https://github.com/B-William/DSPSE/blob/master/spsScript.m


The system with impulse response given by $h[n] = \cos(\pi\sqrt{n})u[n]$ is BIBO-unstable because the sum $\sum_{n=-\infty}^\infty |h[n]]$ diverges instead of being convergent as is needed for BIBO-stability. Note that
for all positive integers $k$, $h[k^2]=\cos(\pi k)$ has value $\pm 1$ and so
$$\sum_{n=-\infty}^\infty |h[n]] = \sum_{n=0}^\infty |\cos((\pi\sqrt{n})|$$
is a sum that contains infinitely many $+1$ terms (and all the other terms are guaranteed to be positive too since $\cos(\pi r)=0$ if and only if $r = k+\frac 12$ where $k$ is an integer, and there is no integer $n$ whose square root is of the form $k+\frac 12$). So, $\sum_{n=-\infty}^\infty |h[n]]$ diverges, and the system is BIBO-unstable.


DTMF Tone Generation in MATLAB not sounding right

You have commented the line that you assign fs. Then the time duration t should be defined with respect to fs.

It should be as follows:

fs = 8000;
t = 0:1/fs:500;


Region of Convergence for pole at z=-1 ?

If the system is causal then the ROC is outside the pole with maximum magnitude. Since there is only one pole at $z=-1$, the ROC is $|z|>1$. Note that the system is not stable since the ROC does not include the unit circle.


Fourier Transforms, symmetry, real/imaginary

Your suspicion is correct. You can show this as follows:

$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{1}$$

$$X^*(\omega)=\int_{-\infty}^{\infty}x^*(t)e^{j\omega t}dt\tag{2}$$

Now if $x(t)$ is purely imaginary we have $x^*(t)=-x(t)$, and, consequently,

$$X^*(\omega)=-\int_{-\infty}^{\infty}x(t)e^{j\omega t}dt=-X(-\omega)\tag{3}$$


How to calculate bandwidth require for OFDM

The channel bandwidth and FFT size alone don't provide enough information to describe the entire structure of the OFDM signal. Recall the following relationship:

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

where $\Delta f$ is the subcarrier spacing, $f_s$ is the sample rate used at the modulator input, and $N$ is the FFT size.

In your LTE examples, the two channel bandwidths must use the same subcarrier spacing, as it takes four times as much bandwidth to carry four times the subcarriers (512 versus 128).

For 802.11a, however, the FFT size is typically 64; with a sample rate of 20 MHz, this yields a subcarrier spacing of 312.5 kHz. Of these 64 potential subcarriers, 12 of them are unused (sometimes referred to as "virtual carriers"). Of the remaining 52 subcarriers, 4 contain pilot tones, so only 48 subcarriers carry information at any given time.

For your notional OFDM system, you can choose the amount of bandwidth required by your system, as required based on your throughput requirements. Choose the sample rate that you will input to the modulator, then choose a convenient FFT size. This defines your number of subcarriers and therefore the subcarrier frequency spacing.


Most of recordings nowadays are being done in stereo - two channels. So when you do a recording in studio, then most of tracks are mono (i.e. microphones, guitars). The main job of sound engineer is to place them in panorama (Left - Centre - Right) to create a nice spatial feeling. For this purpose you have pan knob to set this mono track to the left or right speaker (make it more louder or quieter in either of them). Many tracks are useful when you want to apply different filtering to different instruments or add some reverb, change the volume, etc. In the end all waveforms are being summed and voila!

The main track is obviously a stereo track based on your adjustments of track's pan - that's why you have two channels loaded into your MATLAB.


How can I accurately represent an mp3's wave form in respects to time?

Storing the samples in a square array is weird. You can't represent all song lengths (in samples) as squares of integer array dimensions but yet the complete array is being plotted.

Your code seems to miss one x+=drawStep when transitioning between rows of the array. Neither does it "connect the dots" across the rows. If the last sample of each row is duplicated at the next row then it's fine as a row represents only waveForm.length - 1 samples.

Additionally, it just might be that x and drawStep are not accurate enough. Try making both doubles instead of floats. If there is an improvement but not yet good enough, calculate x by multiplication before each j loop. You should typecast songLengthInSeconds to double too.


How is Bayesian Estimation related to filtering?

In Wiener filtering, you filter a noisy signal to more closely resemble a desired signal that you have access to.  In Bayesian estimation, you take prior knowledge into account to estimate some state given noisy measurements.  In frequency filtering, you just remove frequency content from a signal.


SNR for complex valued signal and when to use SNR and when to use Eb/N0?

Question A:
For the considered system model
$$ y[n] = x[n] + w[n], $$
as already pointed out in your previous question, it holds:
$$ E\{ y[n]y[n]^* \} = E\{ (x[n]+w[n])(x[n]+w[n])^* \} \\= E\{ x[n]x[n]^* \} + E\{ x[n] \}\underbrace{E\{ w[n]^* \}}_{= 0} + \underbrace{E\{ w[n] \}}_{= 0}E\{ x[n]^* \} + E\{ w[n]w[n]^* \} \\= E\{ x[n]x[n]^* \} + E\{ w[n]w[n]^* \}.$$
As a result, your $SNR$ will always be$$ SNR = \frac{E\{ x[n]x[n]^* \}}{E\{ w[n]w[n]^*\}}. $$
So for both of your cases the $SNR$ you specify is correct. It just depends on your definition of $x[n]$ and $w[n]$.

Question B:
For your application I would recommend to use $\frac{E_b}{N_0}$, because the energy per bit depends on the modulation scheme (BPSK, QAM etc.). Therefore comparing the $BER$ in terms of $\frac{E_b}{N_0}$ will be more fair.

You can find further information:

https://en.wikipedia.org/wiki/Eb/N0

https://de.mathworks.com/help/comm/ug/awgn-channel.html

http://www.gaussianwaves.com/2008/11/relation-between-ebn0-and-snr-2/


Impulse property - needed mathematical proof

For any pair of functions $u(t)$ and $v(t)$ we have that

$$\frac{d}{dt}(u(t)v(t))=v(t)\frac{du(t)}{dt}+u(t)\frac{dv(t)}{dt}$$

If we integrate on both sides of the equation we get

$$\int\frac{d}{dt}(u(t)v(t)) \ dt=\int v(t)\frac{du(t)}{dt} \ dt+\int u(t)\frac{dv(t)}{dt} \ dt$$

In this case, $u(t)=\delta(t)$ and $v(t)=\phi(t)$, so

$$\int\frac{d}{dt}(\delta(t)\phi(t)) \ dt=\int \phi(t)\frac{d\delta(t)}{dt} \ dt+\int \delta(t)\frac{d\phi(t)}{dt} \ dt$$

The left side of the equation can be solved easily as

$$\int\frac{d}{dt}(\delta(t)\phi(t)) \ dt = \frac{d}{dt}\int\delta(t)\phi(t) \ dt=\frac{d}{dt}\phi(0)=0$$

as $\phi(0)$ is a constant. We also know that the delta function has the following property

$$\int \delta(t)\frac{d\phi(t)}{dt} \ dt = \left.\frac{d\phi(t)}{dt}\right|_{t=0}$$

So going back to the first equation and replacing with these results:

$$0=\int \phi(t)\frac{d\delta(t)}{dt} \ dt+\left.\frac{d\phi(t)}{dt}\right|_{t=0}$$

And that's the proof.


Using the flip-and-drag method to draw a convolution

Convolution is a linear operator, so let me explain just how to convolve $x_2(t)$ with $\delta(t-t_0)$. The only thing you would have to do afterwards is to convolve the same signal $x_2(t)$ with the other Dirac delta, multiply the function you get by some amplitude, and then add both convolution results.

The Dirac delta has what is known as sifting property:

$$\int_{-\infty}^\infty f(t)\delta(t-t_0)\,\mathrm{d}t=f(t_0)$$

This means that if we multiply a function by an impulse located at $t=t_0$, the area under that differential along the $t$-axis corresponds to the value of the function at $t_0$.

To solve your convolution graphically, imagine the following. Leave the impulses where they are. Then flip your convoluting signal, $x_2(t)$, and take it to start at $-\infty$. Start sliding this function to the right. At first, the convolution (the integral) returns $0$ for all $t$ because the product of $x_2(t)$ and $\delta(t+1)$ is $0$, as they are not intersecting anywhere. But, at some point, they both come across each other. This will continue to happen as long as the sliding version of $x_2(t)$ takes some non-zero value in $t=t_0$. When the triangle goes through entirely, then the convolution will be zero for every other $t$ again, as they will not encounter again.

I'll put some (extremely rudimentary) drawings to make this easier to see.

The first picture corresponds to the cases where the triangle and the impulse don't take any non-zero values at common points. The second one, where the triangle and the impulse are both non-zero at some value (in this case, always at $t_0$). The third picture corresponds to the moment where the triangle has passed through entirely, and thus the convolution is $0$ again.

Try to think about it, and do it on your own. As a rule of thumb, convoluting with deltas is by far the easiest convolution that there can be. In general, if we have something like this (with $A$ being some constant that represents the amplitude of the impulse): $$g(t)=f(t)*A\delta(t-t_0)$$ Then $g(t)$ will just be $f(t)$ but centered at $t_0$ and multiplied by $A$. That means, you just shift $f(t)$ and invert it or scale it depending on the value of $A$.


Signal vs. noise at a particular frequency

The usual scenario is where you want to find if the received signal is noise: $$r(t)=n(t),$$ or if it is a signal plus noise: $$r(t)=s(t)+n(t).$$ It is generally assumed that you know $s(t)$. For example, you could have a radar application where you transmit a pulse and need to find if it was reflected back; or a digital communications receiver, where you want to see if the signal received was $p(t)$ (say, indicating a bit 0) or $q(t)$ (a bit 1).

The way to do this is by calculating the correlation between $r(t)$ and $s(t)$. The idea is that if the signal is present, you have
\begin{align}
\text{corr}(r(t),s(t)) &= \text{corr}(s(t)+n(t),s(t))\\
&=\text{corr}(s(t),s(t))+\text{corr}(s(t),n(t))
\end{align}
which is large. If the signal is not present, then the correlation is small.

This scheme also works in the digital domain, assuming that you sample fast enough. In general you want to calculate the correlation of a large number of samples, so that (by the law of large numbers) the probability of being unlucky and getting "bad" noise samples is small.

In your question, you don't say if you know what the signal is. If you don't, but you know some of its properties, you may be able to look for it. For instance, if the noise is zero mean and the signal is not; or if you know that the signal's energy is concentrated in a certain band (the noise is flat). However, a specific answer is impossible without more information.

If you don't know anything about the signal, detecting it may well be impossible.


This isn't really a MATLAB-specific issue; I see a couple more general questions:

How do you implement a digital IIR filter?

You can apply any general digital filter by convolving its impulse response with the signal that you want to filter. That looks like:

$$
y[n] = \sum_{k=0}^{N-1} x[k] h[n-k]
$$

This works great for FIR filters, but you run into complications with IIR systems, because their impulse responses are infinitely long ($N \rightarrow \infty$ in the above sum). That makes the filter output difficult to calculate using the straightforward convolution sum above. To actually implement an IIR filter, we need a way to express the filter output in a way that is more computationally tractable.

The attribute of an IIR filter that makes its impulse response infinitely long is the system's recursive nature; there is feedback from the output back to the filter input. This means that the output of an IIR filter at any given instant of time can depend upon two separate information sources:

The input signal to the filter, past and present, and
Past values of the filter output signal.

This leads to the difference equation representation of the IIR system:

$$
a_0 y[n] = \sum_{k=1}^{M-1} a_k y[n-k] + \sum_{k=0}^{N-1} b_k x[n-k]
$$

This representation illustrates that we can implement an IIR filter by calculating a weighted sum of $M$ past outputs of the filter and $N$ past inputs of the filter (in practice, $N$ and $M$ could and often are the same; the order of the filter when defined as above is $max(M,N) - 1$). This is a closed-form expression that captures the behavior of the filter fully, suitable for automated implementation, and is the equation referenced in the MATLAB help excerpt that you showed in your question. Note that $a_0$ is almost exclusively assumed to be equal to $1$; I only included it because it is shown in the quoted help text, referenced as a(1).

So to specifically answer your question, as the text suggests, MATLAB implements IIR systems using the above equation. When you design an IIR filter with MATLAB or otherwise, you get two sets of coefficients ($a_k$ and $b_k$ in the difference equation) that define how to weight past filter inputs and outputs in the filter implementation. There are some other subtleties that are sometimes relevant; for instance, the text refers to using a "direct-form II transposed" realization of the filter. You then apply the filter by calculating the difference equation above for each desired output sample.

There are in fact a number of topologies that can be used to implement digital filters, each with their own tradeoffs. DF2T is often used because of its more efficient structure (it contains a minimum number of delay elements) and its improved robustness to roundoff error (since the feedforward coefficients $b_k$ are implemented first).

How does MATLAB design Butterworth filters?

The documentation for the butter function does not specify how it generates the discrete approximation to the analog Butterworth filter prototype. However, by peeking into butter.m, you find:

% step 5: Use Bilinear transformation to find discrete equivalent:
if ~analog,
    [a,b,c,d] = bilinear(a,b,c,d,fs);
end

So it seems that it uses the bilinear transform to map the Butterworth filter prototype to a digital filter realization.