Definition of Perfect Code (and example of Hamming Code) Let's consider the following sentences about coding a message before transmitting it in a communication system (here you find the complete pdf)

First question: can you explain me these two definitions of perfect code? I do not understand what "bit patterns" mean (I think they are the received bits but I am not sure). And I do not understand what do they mean practically.

Now let's consider the Hamming Code (7,4): here the generator matrix G (such that C = B*G, by referring on the upper blocks scheme) is shown:

Second question: I do not understand these two sentences. The rows of G do not appear to be one bit distant from a codeword and the second sentence is completely obscure for me. Can you explain me them?

I give an explanation that avoids the pejorative comments of Engineer.

  First question: can you explain me these two definitions of perfect code?

Those are not two definitions of perfect codes but rather a single somewhat poorly-phrased definition of a perfect code. It should read

  A perfect (binary) $t$-error-correcting code of block length $n$ is a set $\mathcal C$ of binary vectors of length $n$ that enjoys both properties listed below:

  Given any $\mathbf x$, one of the $2^n$ binary vectors of length $n$, there exists a codeword $\mathbf c \in \mathcal C$ such that $\mathbf x$ and $\mathbf c$ differ from each other in $t$ or fewer places, that is, every binary $n$-vector $\mathbf x$ is at (Hamming) distance $t$ or less from some (suitably chosen) codeword in $\mathcal C$.
  Each binary $n$-vector $\mathbf x$ is at distance $t$ or less from one and only one codeword in $\mathcal C$.

It is important to bear in mind that both properties must hold in order for $\mathcal C$ to be called perfect.

This definition is best understood in terms of the notion of Hamming balls in binary spaces. A Hamming ball of radius $t$ and centered at $\mathbf y$ is the set of all binary vectors $\mathbf x$ that are at Hamming distance $t$ or less from the center $\mathbf y$ of the ball. Note that a Hamming ball contains
$$\binom{n}{0}+\binom{n}{1}+ \cdots + \binom{n}{t}$$
binary vectors in it.
So, the first property says that if $|\mathcal C| = M$, then every binary $n$-vector $\mathbf x$ lies in one (or more) of the $M$ Hamming balls that are centered at the codewords, while the second property says that the Hamming ball that contains $\mathbf x$ does not contain any codewords other than the codeword at the center of the ball.  Thus, the $M$ Hamming balls centered at the codewords fill the entire binary $n$-space without any overlap, that is, each $\mathbf x$ belongs in exactly one of the  disjoint $M$ Hamming balls of radius $t$, and so we have a perfect packing of disjoint Hamming balls into the binary $n$-space. This is something that is impossible in our more familiar Euclidean spaces: non-overlapping balls don't fill entire spaces and there are always points in the space that don't belong in any of the balls.

  Second question: I do not understand these two sentences.

The first sentence uses the fact that the $[7,4]$ Hamming code is a $1$-error-correcting perfect code and since every binary vector $\mathbf x$ must belong to a Hamming ball of radius $1$ centered on a codeword, it must either be a codeword (the center of the ball) or be one of the $n$ binary vectors at distance $1$ from the codeword at the center. In binary $n$-space, a Hamming ball of radius $1$ has exactly $\binom{n}{0}+\binom{n}{1} = n+1$ binary vectors in it.

The second sentence points out an addendum to the unstated notion that if $\mathbf c$ is the transmitted codeword and if the channel makes at most one error (flips no more than one of the transmitted bits), then the maximum-likelihood decoding algorithm, which maps received vectors into the nearest codeword, will decode correctly.  The received vector $\mathbf x$ is at distance $0$ (no flips) or distance $1$ (exactly one flip) from $\mathbf c$ and so is in the Hamming ball of radius $1$ centered ar $\mathbf c$. The addendum, which is stated more explicitly, is that if two errors have occurred, then the received $\mathbf x$ lies in a Hamming ball centered at some other codeword $\mathbf c^\prime$ and so will get decoded into $\mathbf c^\prime \neq \mathbf c$ by the maximum-likelihood decoding algorithm, resulting in decoding error. Left unsaid (since these are presumably not things that are dreamt of in the book author's philosophy) are the facts that decoding errors will also occur if the channel makes three, four, five, six, or seven errors.


What does Cooley–Tukey FFT algorithm calculates in order to compute DFT 0
I'm relative new to this subject, I've watched many videos explaining FFT and DFT and read some articles.
I wanted to see how I could implement FFT in C++ and then I found this code, it works but I don't fully understand it, for example, I'm not sure what is he using the nested-for loops for, I thought it's for Matrix multiplication but maybe it's just using the DFT formula for each element of the series, in addition why does he needs to reverse the bits?
if someone could explain in "simple words" the idea behind this code I'll appriciate it because I just started learning it and I've come crossed super complicated stuff which I don't understand.
Here is the FFT Function I'm using:
void fft(Iter_T a, Iter_T b, int log2n)
{
    typedef typename iterator_traits<Iter_T>::value_type complex;
    const complex J(0, 1);
    int n = 1 << log2n;
    for (unsigned int i = 0; i < n; ++i) {
        b[ReverseBits(i, log2n)] = a[i];
    }
    for (int s = 1; s <= log2n; ++s) {
        int m = 1 << s;
        int m2 = m >> 1;
        complex w(1, 0);
        complex wm = exp(-J * (PI / m2));
        for (int j = 0; j < m2; ++j) {
            for (int k = j; k < n; k += m) {
                complex t = w * b[k + m2];
                complex u = b[k];
                b[k] = u + t;
                b[k + m2] = u - t;
            }
            w *= wm;
        }
    }
}

This is the iterative version of the algorithm. Most of the difficulty there is tracing the indices of the subarrays in each stage of the algorithm. For instance bit reversing is used so that the even inputs are moved to the first half of the subarray and odd inputs are moved to the end of the subarray at all levels, to improve memory access locality.
Conceptually the Cooley-Tukey is a recursive algorithm in which each step you split the input in even/odd indices subarrays, and compute the first and second half of the DFT.
$$\begin{eqnarray}
X_k &=& \sum_{n = 0}^{N-1} x_n e^{-2j\pi k n/N} \\
&=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, 2n/N}}_{\textrm{even term}} + \underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, (2n+1)/N}}_{\textrm{odd term }} \\
 &=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, n/(N/2)}}_{E_k} + e^{2j\pi\,k/N}\underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, n/(N/2)}}_{O_k} \\
\end{eqnarray} $$
Then $E_k$ and $O_k$ are themselves the coefficients of a DFT of the decimated even/odd indexed terms of $x$ (of size N/2).
The loop computes two terms one for $X_k = E_k + O_k$ and, (using the periodicity $E_k = E_{k+N/2}$ and $O_k=O_{k+N/2}$), $X_{k+N/2} = E_k + (-1)O_k$ (the $-1$ comes from the exponential mulitplying the $O_k$).
Here is a more readable C++ implementation (from rosetta code)
#include <complex>
#include <iostream>
#include <valarray>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray  odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }
}


Calculate THD From DFT Plot Consider the following DFT plot below, showing the output spectrum of an 9-bit ADC.

I want to calculate the THD (Total Harmonic Distortion). I am wondering how to do it correctly. I am aware of the following formula.

$$
THD=10\log\Bigg(\frac{P_{HD}}{P_S}\Bigg)=10\log (P_{HD}) - 10\log (P_S)
$$

So from my understanding $P_{HD}$ is the sum of all harmonics. But the signal power is shown in dBc, so I can only say how far the power of the harmonics are below the carrier, but I do not have the information about the absolute power levels. Is this correct? Then is it actually possible to calculate the THD, if yes, how?

What is shown in the plot is the PSD (Power Spectral density), which is measured in W/Hz (or dBc). So, in order to get the actual power in a given frequency range, you'd need to multiply with the bandwidth of each bin (which is determined by the length of the DFT window). Let the bandwidth of each bin be $B$.

Then, your formula becomes

$$
THD = 10\log\left(\frac{P_{HD}}{P_S}\right)=10\log\left(\frac{\sum_i B\cdot S_i}{B S_{HD}}\right)
$$

where $S_i$ denotes the Density at the $i$th harmonic and $S_HD$ denotes the density at the fundamental. As you see, B cancels in the expression, so you can ignore it. Note, that $S$ are in linear scale, (i.e. you'd need to transform the measurements in the diagram to linear scale before summing them up etc.)


If you had supersonic hearing, would you be able to hear the carrier wave of an AM radio station? Would you be able to hear the carrier wave of an amplitude modulated station, if you could hear in the MHz range, or would the high frequency be filtered out after rectification?

The AM or medium wave band extends from 526.5 kHz to 1606.5 kHz in Europe or 535 kHz to 1705 kHz in the US.

Ultrasound (not supersound! :-) frequencies range from 20kHz and up:

However, a bigger problem is that AM stations broadcast using electromagnetic radiation, whereas sound travels by acoustic waves in the air.


What is beat-synchronous chroma? I have a pretty good understanding of what chroma is when represented against time, but I don't understand what it means when one represents chroma against beats

Short answer

Normally a chromagram constits of chroma feature vectors per time frame. If you want to have a musical time axis in beats instead of a physical time axis in seconds, you have to get beat information (e.g. by a beat tracker) and average all chroma features that belong to the time span of a beat.

Long answer

Since the comments asked for a longer answer, here it is: A chromagram presents the energy that is contained in the twelve pitch classes over time. Starting from the STFT, you could compute a chromagram version by summing up all frequency bins that correspond to a particular pitch class, individually for each time frame. Since the STFT uses a constant hopsize, you end up with a chromagram which is constant in physical time. Since chromagrams are used for music applications, it is sometimes desirable to have a musical time axis, e.g. in beats. Since the tempo of music can change, the physical time span of an beat can vary. For converting physical to musical time, you could annotate the beat positions within the music or apply an automated procedure, kown as beat tracking. Then you could simply average the chroma features that belong to the time span of a beat.

In more technical wordings: $\mathcal{C} \in \mathbb{R}^{N \times 12}$ is your normal chromagram matrix with $N \in \mathbb{N}$ time frames. $\mathcal{B} : \mathbb{N} \rightarrow \{s| s \subset \{1,\dots,N\}\}$ is a function that gets the index of a beat as input and outputs a set of all time frame indices belonging to this beat index. You want to compute the beat-synchronous chromagram $\mathcal{C}^\text{beat} \in \mathbb{R}^{B \times 12}$ for $B \in \mathbb{N}$ beats. Usually $N \gg B$. You can do this simply like this: $\mathcal{C}^\text{beat}_{b, c} = \frac{1}{|\mathcal{B}(b)|} \sum_{n \in \mathcal{B}(b)}  \mathcal{C}_{n, c}$ with beat index $b \in \{1,\dots,B\}$ and pitch class index $c \in  \{1,\dots,12\}$.

This is not the place for comprehensive literature, but only two references for general chromagrams and beat-synchronous ones:

Meinard Müller: Fundamentals of Music Processing, Springer 2015. (Chapter 3.1)

Daniel P. W. Ellis and Graham E. Poliner: Identifying Cover Songs with Chroma Features and Dynamic Programming Beat Tracking, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Honolulu, 2007, pp. IV-1429–IV-1432.


Difference Between Correlation / Convolution and Matrix Multiplication Can anyone please clarify the difference between correlation/convolution and matrix multiplication? As I thought either convolution or correlation is similar to matrix multiplication. I read this question difference between convolution and multiplication, that is about simple multiplication.I am asking about matrix multiplication.
Thanks in advance.

Well I will try to explain.

Let us first discission in time domain:

1) Let us say you have two signals, x and y. By Convolution in time domain, you mean that you flip(invert) one of the signals(lets say ,y) and then slide it over the other(x). By corerlation, you just don't flip them, rest is the same.

2) Now when we talk about correlation, we first discuss 1-D case, where you have just one signal(x) and then for multiple signals(x,y...)

One Variable:

The variance for x will tell you how spread or varied that signal is, i.e. a small will correspond to signal which is changing very much.

As as engineers, we like that just at looking at value we can make prediction about the nature of the signal. For that reason, variance is not a good estimate, as you cannot compare signals easily by just looking at their variance. For that reason we introduce Standard Deviation, which is infact the underroot of variance, but by looking at it you can easily compare the nature of two signals.

Multi Variable:

As variance is in one signal case, in mutli variable, where you have X,Y ...etc signals, we replace it by Covariance. A covariance will tell you how much X and Y vary with each other. But again like variance, we cannot compare just based on co-variance, so we introduce correlation(Standard deviation in one signal case). Now by looking at correlation, you can compare the signals. And then a correlation coefficient,which is just the normalization of correlation, and it ranges from 1 to -1. If the correlation coefficient between 2 signals(x and y) in this case, near to 1, means these both signals are same or same behavior, whenever one signal goes up, the other also and vice versa. 0 means that they have no relationship at all and -1 means inverse relationship.

Computing:
Now to compute the convolution/correlation, we use matrix multiplication as a tool, just like  you use variable multiplication for calculating variance or mean in 1 signal case. In multi dimension case(x,y,...) you use matrix usually instead of variables to make your computations easy(especially on computers).

Now in frequency domain:

For calculating, convolution, in time domain, you had to flip and then slide(slide means that each and every element of variable needs to be slided) which is computationally very expensive. Instead, in frequency domain, the same result as that of convolution is achieved just by multiplication and yes when you signal is in freq. domain, you simply multiply the signal(X)-which is matrix with Signal(Y), which is also a matrix.

So, now you will be able to understand that, Yes convolution is same as matrix multiplication(where matrix X and Y matrix of signal) but ONLY IN FREQUENCY DOMAIN.

I hope it helps.

P.S: Simple mulitplication is done in theory, when you are writing equation, but when you have a real signal , then you sample it and it gives a matrix, which you further multiply in matlab.


Problems with IFFT not being symmetrical I have two signals, a measurement and a reference which I have performed an FFT on.  They have both been windowed with a Hanning window, and now I would like to deconvolve these to get the impulse response.

I have followed these steps:

Divide my reference FFT by my measurement FFT. (Point division of complex numbers).
I then take the conjugate of this resultant array.
I then perform a forward FFT on this.
I then take the conjugate again.
I then divide by the length to normalise the signal.

If I am to understand correctly, this should give me the symmetrical deconvolution of the two signals.  However, for my input of 1024 samples, I get the first 512 samples equal to 0, and then the second half of the samples are a little bit noisy with a definite impulse at the expected delay time.

My measurements are actually completely electronic so should be noise free.  I can confirm that my two input FFTs are symmetrical (I am using the full FFT, not just the positive values), I am not doing any 'FFT Shift'.

If anybody has any suggestions or hints, I would greatly appreciate it.

Thanks

Here is an example using Matlab/Octave code. There is no need to conjugate.

x = rand(1000,1);
y = rand(1000,1);
z = conv(x,y,'full');
y2 = ifft( fft(z,1999) ./ fft(x,1999) );
norm(y - y2(1:1000)) / norm(y)

ans =

   1.5654e-13


FFT in MATLAB is shifted in frequency I have a signal in MATLAB on which I'd like to perform an FFT. My signal is stored in s, and I use the code below (inspired by the MATLAB help):

L = length(s);
nfft = 2^nextpow2(L);
S = fft(s,nfft)/L;
fftf = 1/Ts/2*linspace(0,1,nfft/2+1);
ffts = (2*abs(S(1:nfft/2+1)));

s is a GMSK-modulated bit sequence, that is, it varies between $f_c-2400$ Hz and $f_c+2400$ Hz in my case when transmitting 0 or 1, respectively. $f_c$ is set to 100kHz. For a long input, say, 250 bits worth of 1's and 250 bits worth of 0's, I get what I expect, see the first image below.

If I instead choose a low number of bits, say 10 1 bits followed by 10 0 bits, I get as expected, but it is shifted down to ~90kHz instead of centered at 100kHz. This is something I can't quite understand - it seems changing the sample rate and length of the FFT changes absolutely nothing.

Can anyone explain to me why? Thanks in advance!

Long data:

Short data:

The code used to generate the signal and FFT:

%% Configuration
clear; clf;

DataRate = 9600;          % 9600 kbps for AIS
N = 100;                  % Samples per bit
Tb = 1/DataRate;          % Bit period
Ts = Tb/N;                % Sampling period
BT = 0.5;                 % AIS spec, time-bandwidth product
Ftrans = 100e3;           % Frequency of "transmitted" signal

num = 200;
Bits = zeros(1,num)+1;
Bits = [Bits zeros(1,num)-1];
clear num;

%% Modulation

% Prep a time axis from -2Tb to 2Tb
t_g = -2*Tb:Ts:2*Tb;

% Gaussian response to rectangular pulse [Haykin4th, p. 397]
x = pi*sqrt(2/log(2))*BT;
gr = 1/2*(erfc(x*(t_g/Tb-1/2))-erfc(x*(t_g/Tb+1/2)));

% Truncate to 3Tb, pulse centered at 1.5Tb
gr = gr(0.5*N+2:3.5*N+1);

% Normalize
% when integrating, we want to end at 0.5 (phase changes by 0.5pi)
% so, we want sum(y)=0.5 -> normalize by sum(y) and divide by two.
gr = (gr)./(2*sum(gr));

% Generate the Gaussian filtered pulse train by centering a "Gaussian
% rectangle" on each bit, and adding inter-symbol interference
f = zeros(1,(length(Bits)+2)*N);
for n = 1:length(Bits)
    f((n-1)*N+1:(n+2)*N) = f((n-1)*N+1:(n+2)*N) + Bits(n).*gr;
end

% Since gr corresponds to changing the phase 0.5, multiplying by pi and
% integrating gives the resulting phase.
theta = pi*cumsum(f);

% Prep I,Q
I = cos(theta);
Q = sin(theta);

% Transmitted signal, shifted to ftrans
t = linspace(0,length(Bits)*N,length(I))*Ts;
s = -sin(2*pi*Ftrans.*t).*Q+cos(2*pi*Ftrans.*t).*I;

%% FFT

L = length(s);
% faster w/ a pow2 length, signal padded with zeros
nfft = 2^nextpow2(L);
% do the nfft-point fft and normalize
S = fft(s,nfft)/L;
% x-axis from 0 to fs/2, nfft/2+1 points
fftf = 1/Ts/2*linspace(0,1,nfft/2+1);
% only plotting the first half since its mirrored, thus 1:nfft/2+1
% why multiplied with 2?
ffts = (2*abs(S(1:nfft/2+1)));

%% Plotting

% FFT PLOT
plot(fftf/1e3,ffts);
title('FFT of transmitted signal S(f)');
set(gca,'xlim',[Ftrans/1e3-20 Ftrans/1e3+20]);
ylabel('|S(f)|');
xlabel('Frequency [kHz]');
grid;

Adjusting the sample frequency by changing N seems to have no effect - but changing num from e.g. 10 to 100 (changing the number of bits) clearly shifts the plotted spectrum closer to 100kHz.

Unfortunately I don't have Matlab now but I'm thinking your problem has to do with linspace. Remember linspace gives the actual first and last values always. When you're doing the FFT binning I don't believe you want the boundaries to be pinned that way. I'd try using the colon operator instead and appropriate starting and delta values.

EDIT:
If you look at the definition of the DFT it is:
$$X(k)=\sum_{n=0}^{N−1}x(n)e^{j\frac{2πkn}{N}}$$

This shows that the first DFT output at $k=0$ is at baseband. Every subsequent frequency is at an interval of $F_s/N$. So if you had 10 points the bin centers would be from $0...(\frac{N-1}{N}F_s=0.9F_s)$. For 1000 points it would span $0...0.999F_s$. Linspace, on the other hand always leaves the first and last points exact. You need the start and spacing exact, or to adjust the ending point with linspace.


Power spectrum estimate from FFT I'm plotting the FFT power-spectrum of a signal in MATLAB.
I uploaded the 8000 samples time-series signal in a text file here:
http://wikisend.com/download/896484/signal.txt

I'm using the following code:

Fs = 500; % sampling frequency
L = length(y); %number of samples

complex = fft(y)/L; % complex signals
f = 0 : Fs/L : Fs/2; % frequency bins
amplitude = 2*abs(complex(1:L/2+1)); % amplitudes
pow = (amplitude).^2/2; % power

semilogy(f,pow,'-ro');
grid on;

And I'm having the following plot:

The confusion is, where I work they are using their own modified FFT function and if I run that code I obtain the following plot:

If you check the power in above figures the linear power ratio is 16 between my code and their's. And for other measurement I'm getting the ratio 60, 120 or even 1.

I'm confused and need help if one can plot FFT for me and see if my code and plot correct.

The difference is in the scaling of the power spectrum. I suppose you will find that the difference in scaling always equals $L/F_s$, which for the numbers in your question is indeed $8000/500=16$.

Take a look at Power Spectral Density Estimates Using FFT for the correct scaling. If you normalize the FFT result by the FFT length, you need to scale the squared magnitude of the normalized FFT by $L/F_s$. Furthermore, you need a factor of $2$ if you throw away the negative frequencies (I see you did that). Note, however, that you shouldn't apply a factor of $2$ for the bins at DC and at Nyquist, because they are not mirrored at negative frequencies.

Here is a simple Matlab code from the above quoted mathworks page for computing a periodogram-based one-sided power spectrum estimate using the FFT (my comments):

Fs = 1000;                          % sampling frequency (Hz)
N = length(x);                      % even! (=> bin N/2+1 is Nyquist)
xdft = fft(x);
xdft = xdft(1:N/2+1);               % DC to Nyquist (one-sided)
psdx = (1/(Fs*N)) * abs(xdft).^2;   % periodogram scaling
psdx(2:end-1) = 2*psdx(2:end-1);    % scaling for one-sided PSD
                                    % note: no factor 2 at DC and Nyquist


Accuracy of time domain signal when one frequency is absent Suppose the known frequency spectrum of a signal is:
$$
G(n\Delta f)=9,4,3,5,1,6,...
$$
where $n=0,1,2...$ and $\Delta f=1$ Hz

If, for some reason, we cannot sample $G$ at $n=0$ (i.e. $G(0)=9$ is not available), but can sample the frequencies everywhere else, including at finer $\Delta f$'s.

What can we do to minimize the error of the recovered time-domain signal $F(t)$?

For example, how fine $\Delta f$ should be? How many frequencies needed to be sampled? Can we use the discrete Fourier transform or should we use the continuous Fourier transform by interpolating the shape of $G$? etc.

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

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

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

Another more recent area is known as "Matrix Completion"

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

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

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

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

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


How does transition width of a low-pass filter affect a message which is well within the pass-band? I'm taking a communication systems course and we use GNURadio in the lab. In an exam, we had a simple question where we were asked to demodulate a Dual Side Band (with carrier) AM signal in GNURadio. Carrier frequency was known to be precisely 500KHz and message was audio (known to be < 15KHz). The incoming signal has sampling rate of 1.92MHz.

EDIT: i.e incoming signal is given to be of the form $(1 + k_am(t))cos(2{\pi}{f_c}t)$ with $|k_am(t)| < 1$

To do that, I took the absolute value of the incoming signal, and used a Low pass Filter with cutoff frequency 17KHz and transition width 100Hz to obtain the audio clip. I was listening to the audio sampled at 48KHz (with proper decimation of 40 set in the Low Pass filter block). However the received audio clip was very jittery. I could listen and identify the song, but there were periodic higher frequency spiky noises.

The TA observed this, changed the transition width of the LPF from 100Hz to 1KHz and it was  completely fine now - I could hear it neat and clear! I argued that I only tried to make the filter better by making it more ideal - brick-wall type, and also, this transition is happening where I don't have my message (transition occurs at 17KHz and my message < 15KHz). He replies that I pay a price when I make a filter transit faster, and added I would understand that when I do the digital communications course.

Could someone kindly explain to me, what price I pay (in digital domain) when I decrease the transition width of the Low Pass filter? Thanks in advance.

For an FIR filter there are 3 main components that determine the filter length for equiripple designs:

Passband ripple
Stopband attenuation level
Transition width (Width from the edge of the passband to the edge of the stopband)

For other filter designs the filter order may be related to flatness of the passband and the rate of fall-off in the stopband (and initial stopband attenuation). The transition width is often expressed as a fraction of the original sampling frequency.

In your case the original fraction is $\frac{100}{1.92\times10^6}\approx 52.1\times10^{-6}$ which is very small and is going to require a very long filter. Depending on the algorithm used to design the filter, the algorithm may run into numerical difficulties. As an example, designing a 2000 tap FIR equiripple filter using Matlab's Remez type routines is often problematic.

Increasing the transition width by a factor of 10 will roughly reduce the filter length by a factor 10 as well. This is just a rough approximation and results will vary based on the design routine and the filter specs being used.

If you really need such a small transition width, it is better to perform you decimation in steps, e.g, break it down to decimate by 10 and then by 4. In the first decimate by 10 you can use a wider transition width. The transition width in the last filter will be the same as the original requirement, but because of the reduced sampling rate the filter length will be much smaller. Crochiere and Rabiner's "Multirate Signal Processing" goes into more details.


Convolution between two unequal length vectors Considering an example, where the model is
\begin{align}
x[n] = \mathbf{h}^\mathsf{T}\mathbf{u}[n] + w[n].
\label{Eq1}
\end{align}
$\mathbf{h} = [h_0,h_1,\ldots,h_{M-1}]^\mathsf{T}$ of length $M$ which represents the coefficients of the LTI SISO channel and $\mathbf{u}[n] = [u[n], u[n-1], u[n-2],\ldots,u[n-M+1]$ is the input to the channel. Let $\mathbf{v} = [v_0,v_1,\ldots,v_{M-1}]^\mathsf{T}$ of length $M$ be the coefficients of the equalizer.

If the length of the equalizer is different, say the channel length is $M$ and equalizer length is $Q=2M-1$ then how would the convolution happen since there will be a mismatch in the number of elements of $\mathbf{h}$ and $\mathbf{v}$

I can do conv(h,w) in Matlab where h is the coefficients of the channel and w denotes the coefficients of the equalizer. Assuming known values of h = [1,0.2,0.6] and the estimated coefficients of the equalizer as `w = [0.8,0.1,0.2,0.3,0.5]'. How would the convolution work? Can somebody please show first few steps. The first few sample will be zero based on my calculation. But i am not sure.

The discrete convolution sum operation is not restricted to equal length vectors. You can, and most of the time you do, convolve two different signals of arbitary lengths. Your confusion is probably with something else. The equalizer length can be different than that of the channel model length. That should not pose a problem but it would of course effect the equalization performance.


Given input and output time series, obtaining filter design Say I have 5 minutes of input and output audio. One method I know is to do FFT on windows from input and output. Then divide FFT output into bins and find average energy in various bands for all the windows. After this, take ratio of energy bands from output windows with corresponding energy bands from input windows. And then take average of all these ratios. This will give an approximation to FFT(Y)/FFT(X) i.e. transfer function of the filter.

Are there any other methods to find transfer function given input and output time series?

If your system is linear, then you can find the impulse response $h(n)$ of the system using the convolution sum $y(n)=\sum_{k=0}^{N} h(k)x(n-k)$, which can be written in matrix-vector form as $y=Xh$, where $X$ is a convolution matrix of the input data and $N$ is the order of the system, say $N=10$ for example. Since you know $X$ and $y$, you can find $h$ via the Pseudo inverse, assuming you have created an over-determined system $y=Ah$. The pseudo-inverse is $h=(X^T X)^{-1}X^T y$. After finding $h$, find the FFT of $h$ or take the Laplace transform of $h$ to get the transfer function. I hope this helps you.


Comparing an audio signal with a previous 'taught/profiled' signal I'm working on a audio sample recognition project, where first you profile a sound and then run the sound detection function.

Currently during the 'profile' stage I'm grabbing copies of the magnitude frequency bin from the sound you want to detect (for a changeable pre-set amount of FFTs). During the detection phase I'm trying to use these to compare the similarity of the incoming identical length FFT, to see if they 'match'.

What would be some good methods to compare the two FFTs given that I have a 'learn' phase which can be used. I've tried comparing the maximum frequency bin, which works well to an extent. Voice recognition is easily fooled, with a tone generator (Which makes sense)'.

I'm sampling at 16000Hz with each FFT at 2048 samples long, equating to 128ms of data. Which I believe gives me a frequency range of ~8000hz.

FFT won't work because it captures too many details which are rather irrelevant to the task. You need to extract a representation of sound which is more or less invariant by pitch and amplitude changes (properties of sound which are irrelevant to your problem). So you need a "blurred" version of the spectrum - and MFCCs are exactly that.

The textbook-ish, baseline system for speaker verification is to extract MFCCs (first coefficient excluded) on the "profile" recording; train a Gaussian Mixture Model on this data. The more data we have in the profile, the better (ideally, the most common phonemes should all be represented).

To score a new recording (during the "recognition" phase), extract the MFCCs, and compute the likelihood of the model given this data (you'll have to normalize by the length of the recording). Compare the score to a threshold - set according to the false positive rate you want to achieve.

For speaker recognition, you train one model per speaker to recognize; and select the speaker associated with the model with the higher score.


What is the entropy for these cases? This question stems from an article,

"Entropy estimation of symbol sequences"  download link

where the abstract mentions the need for using symbols in information theory.

  Random chains of symbols (zeros and ones or characters) $s_1, s_2,\ldots,$ drawn
  from some finite alphabet
  appear in practically all
  sciences. Examples include spins in one-dimensional mag-
  nets, written texts, DNA sequences (here symbols are alphabets), geological records of
  the orientation of the magnetic field of the earth, and bits
  in the storage and transmission of digital data.

The paper presents technique for finding the entropy of symbols. My Question is that if entropy of symbols is important then what about the entropy of real valued numbers?
From the perspective of entropy, is there a proof or some information in text where it provides information that the entropy of random symbols is greater or lesser than the  entropy of real numbers.

Q1: In other words, what would be the entropy of random variables having binary values and entropy of real valued random variables?  A fair coin takes 2 values -0/1 with equal probability and its entropy is 1 bit. So, is 1 bit the maximum entropy? Would entropy of binary random variables be less or more than entropy of real valued random variables and how to show?

So, I would like to know if there is a proof or material from which it can be inferred if the entropy of characters is lesser/ greater or same with  the entropy of real numbers.

Q2: I want to know why symbols are used in coding and other applications such as storage from the viewpoint of entropy and information theory.

The entropy $H(X)$ of a continuous random variable $X$ is infinite. Proof is trivial (note that we can, without loss of generality, use the natural logarithm, since any other logarithm is the same but with a factor: $\log_b(x) = \frac 1{\ln(b)}\ln(x)$):

\begin{array}
\\H(X) &= E\{I(X)\}\\
&= \int\limits_{-\infty}^{\infty}{-\ln\left(P(X=x)\right) dx}
\end{array}

The point here is that for any really continuously distributed random variable, the probability that a specific real value gets "hit" is $P(X=x)\equiv 0 \,\forall x$, and hence, the information contained in that event happening would be infinite. You hence cannot directly compare discrete alphabets with dense intervals (such as the whole $\mathbb R$).

Hence, people invented differential Entropy. For your real-valued (=1-dimensional) case, this is nothing but the obvious:

$$h(x) = \int\limits_{-\infty}^{\infty} f(x)\log_b\left(f(x)\right)dx \text{,}$$ with $f$ being the continuous probability density function of $X$.

Now, all you need to compare your discrete case with the continuous is model the discrete PDF for the alphabet $\mathbb X$

$$
f_\text{discrete}(x) = \sum\limits_{x_i\in\mathbb X} \delta(x - x_i)P_\text{discrete}(x_i) $$

$\delta$ is the dirac function.

Caveat not everything you know about the discrete entropy applies to differential entropy. Read the Wikipedia Article!


Determining if a system involving integrals is time varying or time-invariant I've been studying for a signals and systems class coming this fall and am trying to understand the following problem.
$$$$
Show whether the following system is time invariant or time varying:
$$S_1:y = T_1[x],\;\; y(t) =  \int_0^t{[x(\tau)]^3\,d\tau},\;\;t\geq 0$$
solution:
$$z(x)=t[x(t-\sigma)] =  \int_0^t{[x(\tau-\sigma)]^3\,d\tau}$$
$$=  \int_{-\sigma}^{t-\sigma}{[x(\tau')]^3\,d\tau'}$$
Where the change of variable $\tau-\sigma \rightarrow \tau'$ is applied. and
$$y(t-\sigma) =  \int_0^{t-\sigma}{[x(\tau)]^3\,d\tau}$$
So I know that the answer is that $z(x) \neq y(t-\sigma)$ and thus the system is time varying. So my question is why doesn't 0 get subtracted by $\sigma$ in integral expressed in $y(t-\sigma)$?

Furthermore if I switch the limits of integration like so:
$$S_2:y = T_2[x],\;\; y(t) =  \int_{t}^{t+1}{[x(\tau)]^3\,d\tau},\;\; t\geq 0$$
Would $S_2$ be time invariant?

The fixed bound at 0 in $\int_0^t$ indeed suggests a time-variance.You only need a counter-example for $S_1$, like the constant $x(t) = 1$.  $y(t+1) = t+1$, while $S_1(x(t+1)) = S_1(x(t)) = y(t)$.

For $S_2$, the two bounds moving at the same time are promising, and a standard variable change shows the time invariance.


What does mean when we say mel scale reflects how humans hear musical tone According to the standard mel frequency conversion formula, the mel frequency corresponding to 10000 Hz is 3073.22.

Does it mean human ears perceive 10000 Hz frequency of any sound as 3073.22 Hz?

The Mel scale is not in Hz units but rather Pitch level.

The point of the Mel scale is to return an estimate of the perceived pitch and it is more useful when looking at differences rather than single values.

So, if you notice the Mel scale as it extends past 3500 Hz, it tends to a straight line. This means that a human being would rate differences in frequency as expected (doubling of frequency results in proportional increase in the Mel scale).

But below 3500 Hz, the curve is more sensitive to frequency variation.

This is part of more general observations in perception where perception and stimulus are in logarithmic relationship.

Hope this helps.

EDIT:

In terms of an example, you can think of the siren of an ambulance. The Mel scale "predicts" that if you were to transpose the siren's centre frequency to a higher one, you might not be perceiving it as going up and down that much while the modulation level is exactly the same. And that's something to think about. Objectively same quantities evoking different perceptions. Taken to the extreme, there could be a frequency range that the Mel "line" becomes so shallow that if asked a listener would say "all I hear is a constant tone" although the siren still goes on at its regular Neeee-Noooo-Neeee-Nooo sound.


DFT of discrete signals, why do we only analyze frequency bins equal to number of input samples? If we have a signal $x[n]$ such that we have $N$ samples i.e. $n=0, \ldots, N-1$, then when we analyze the DFT $X[k]$ we only analyze for $k=0,\dots,N-1$ as well.

Why is the range of $k$ tied to the range of $n$?
Why don't we analyze more, or less number of frequencies?

A DFT is an invertable (square matrix) transform, thus the complex result vector is the same size as the complex input vector.  But you can analyze as many (additional) frequencies as you want by interpolating the results after the DFT (which requires additional computation that doesn't really increase separation resolution, but can produce a smoother looking plot).

You can analyze less if you want by using the Goertzel algorithm, which can be used to compute any single frequency result bin of a DFT.  (But calculating less than N orthogonal bins won't be invertable or "complete".  And any more than log(N) bins will likely be a lot slower than just doing a full FFT and throwing away what you don't need.)


How do I filter out this horizon from my contour? I am trying to identify the T-rex and obstacles (cactus and bird) from Google's no-internet game.

What I have done is:

Apply GaussianBlur filter
Apply Canny
Apply Dilate
Apply erode
findContours that returns all external contours, RETR_EXTERNAL

My approach generally works, except that when there is a small bump on the ground next to a cactus, my script will consider the T-rex, the bump on the ground and the cactus as a giant horizontal contour, how do I filter out this horizon? The blue box is the outline of the contour I found.

I tried reducing the size of the image I am taking vertically so that the horizon would be left out, but it introduced a different problem: T-rex legs would be cut off and occasionally, T-rex would not be considered as a contour as it is not closed.

When there is no bump on the ground next to a cactus, my script produces the expected result.

I have found a function provided by opencv2, threshold().

With this function, I can apply different combinations of threshold values so that this horizon can be filtered out.

The current values I am using is:

cv2.threshold(img, 100, 255, cv2.THRESH_BINARY)


Should we upsample the channel when upsample the signal I have an OFDM system with  Number of sub-carriers $N= 1024$, modulated data using $QAM$ modulation is transmitted via those subcarriers as follows:

modulation --> ifft --> adding CP --> upsampling --> conv channel --> adding noise

Then at the receiving side, the received signal is processed as follows, discarding the channel

delay --> downsampling --> CP removal --> fft --> MMSE equalizer --> demodulation

The issue which I am facing in that system is in the equzlizer step,  When using the original channel, I can not get the performance back. Howenver, when I estimate that channel after the step of fft and then use the estimated channel, The performance is OK !!

Why does that issue happens? I think because of the upsampling step, because when I don't upsample the signal, that becomes fine but when I use the upsampline and use then equalize using the original channel, I can't get the performance back.

In simulation, we do upsample to simulate the "analog" signal which is a continuous signal in reality but we represent it as a high-time-resolution in our digital/discrete simulation. When we upsample our data before going through the channel, we should also upsample the "impulse response" (IR) of the channel. It means that the timing-base of both signal and IR of channel must be same. For example, without upsampling, take for example of a channel with an IR of {4,0,0,3,0,1.5,0,0,0}. First, We may apply our signal (without upsampling) to this channel. Then, we upsample our signal with a factor of 2 and want to apply that to this channel again. this time we have to upsample our channel's IR with a factor of 2 just like the signal. So, the upsampled channel'IR would be: {4,0,0,0,0,3,0,0,1.5,0,0 ....}. So, the time-base or time-resolution of both signal and channel should be equal.


Fourier series representation Conceptually when we want to represent a peroidic series, for example a pulse train, we find the Fourier coefficients and get a representation in the time domain.

However what is conceptually wrong with using say an infinite sum of time shifted rect functions to represent it?

Sorry I had a formatting problem, so i am putting this in the bottom post...

My method is as such:

Assuming we have a periodic pulse $x(t)$ such that $x(t) = 1$ for $0 < t< T$ and $0$ for $T < t < T_p$; thus $x(t)$ has a period of $T_p$.

Finding the fourier coefficients Ck via:

$$Ck = \frac{1}{T_p} \int_{-\infty}^{\infty} x(t) * e^\tfrac{-j2\pi kt}{T_p}$$

over 1 period, and thus we can represent x(t) as:

$$x(t) = \sum_k C_k e^\tfrac{j2\pi kt}{T_p}$$

and doing the Fourier transform we will get this form:

$$X(f) = \sum_k C_k \delta(f-kf_p)$$

which is discrete.

However if we consider $x(t)$ to be of this form:

$$x(t) = \sum_n \text{rect} \left(\frac{t-nT_p}{T} \right) $$

Applying fourier transform of $x(t)$ to get (in the form):

$$X(f) = \sum \text{sinc} * e^{-j*W}$$ -- (2)

where sinc() is due to the FT of rect and the e^(-j*W) comes out due to the time shifting property of FT.

Comparing X(f) in (1) and (2), we see that 1 is discrete and the other continuous.

However they come from the same x(t), so isn't this a contradiction?

Sorry for the long post.

However what is conceptually wrong with using say an infinite sum of time shifted rect functions to represent it?

There is nothing conceptually wrong with it.  Fourier transforms decompose a signal into a sum of complex sinusoids, but you can also decompose a signal into many other things, which may be more useful in certain applications.  The Haar wavelet transform, for instance, breaks up a signal a sum of rectangular pulses:

 source

We use sinusoids in lots of applications because it makes the most sense in those applications.  For instance, why do we almost always decompose audio signals into sinusoids?  Because our cochleas do the same thing:


What is the amplitude of the impulse resonse I'm learning how to construct a discrete signal as a sum of impulse responses of the system. The final formula is then:

$$
y[n]=  \sum_{k= -\infty}^\infty x[n] h_{k}[n]
$$

My question is, what is the amplitude of the impulse response, because i thought an impulse response was already a scaled response to the unit impulse inputted into the system? But here they are multiplying again with a scaling factor of x[n]? So is the unit impulse's amplitude just 1 then?

In case of discrete signals, the amplitude of unit impulse is 1.

It does not mean that your impulse response should have unit amplitude. The impulse response can have any real valued amplitude in it.

Note: The area under unit impulse in continuous time domain is 1 not the amplitude. The amplitude of a unit impulse approaches infinity as the area under remains constant (=1).


How to Apply Gaussian Filter to Raw Audio Files in Python? I am recording raw files in python to later break them into phonemes, but the noise in the surrounding environment is hampering the result. So, there is a need to apply a filter to the recorded raw audio files. How can this be done ??

Do it using the following steps:

Look for the optimal frequency for the Low Pass Filter. Usually close as possible to the bandwidth of your desired signal.
Design an LPF filter according to the frequency above. If you're after Gaussian based filter you need to optimize the STD ($ \sigma $) parameter.
Apply the filter either using convolution, Using Numpy's convolve() function (Only in case of FIR Filter) or Scipy's lfilter() function (Which, in case of FIR Filter does convolution as well yet can also handle IIR Filters).

Remark
Gaussian based filters aren't optimal for the task you are after (Their passband isn't flat).
I'd would design a filter using Remez Method (See Scipy's remez()) which is Equiripple Method or using Chebyshev Filter Type II (cheby2()) and Butterworth Filter (butter()) which has flat passband (Namely they won't do any harm to the audio content).


Haar Wavelets vs STFT for graceful audio quality degradation (newbie question) I'm playing around with FFT (using Naudio library) and with Wavelets (using TarsosDSP library) trying to get better understanding of both of them.

I've noticed that when I zero out FFT/Wavelet coefficients (for the same window size - and then perform inverse transform and listen to the result) starting from low to high energies and compare the loss of audio quality Wavelets seems to suffer less. They have less noise for the same number of zeroed out coefficients (or require less non-zero coefficients to have roughly the same amount of noise) and require less non-zero coefficients for audio still to be recognizable.

Is it inherent property of Haar Wavelets or have I messed up somewhere and FFT should perform equally well in this regard?

For some non-stationary 1D signals, containing sudden changes, or a lot of noise, wavelets can have a close to optimal performance. Here, by wavelets I mean short oscillating shapes spanning several frequency ranges or scales.
Wavelets are somehow better than Fourier to that respect, as they can deal with peicewise regular + jumping signals.

However, this covers a huge range of possibilities for continuous and discrete tools.

Hence, I am confused with your "Haar Wavelets vs STFT" title. Haar wavelets can be orthogonal, STFT very redundant. So one cannot compare them directly in terms of sparsity and quality of approximation.

I will thus hypothesize that you performed a Hamming window of some audio signal frame of length $L$. Then applied either an FFT or an orthogonal Haar discrete wavelet transform on the result. Then performed a nonlinear approximation by keeping the $N$th largest coefficients, with $L \ll N$.

The result you evoke is a bit surprising to me, and I am a wavelet supporter: it seems a little too good for Haar wavelets. Haar wavelets are often  not the best fit for audio, and few audio coders use wavelets. What makes me wondering about your results is:

the length of your frame, or is the signal sufficiently nonstationary to win by KO on Fourier?
"starting from low to high energies": I have problem understanding what it means with wavelets. Are you performing a fair treatment on Fourier ?
magnitude and complex symmetry: Fourier coefficients are usually "doubled" on real signals. Did you take that into account? Are you sure you have reconstructed Fourier coefficient fully, and not from magnitude only?

All in all, a sketch of your code could help us deciphering the results.


understanding FFT2 function of Matlab I am trying to understand what happen when we take FFT2 of a matrix in a matlab.

first have a look at this simple example,

a=ones(8); %8x8 Matrix of Ones
fft2(a)

it will generate the following output,

    64     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0

in above matrix , 64 is the sum of all the matrix values and it is the DC coefficient.
Does that mean the 64 is calculated as 8*8 = 64 ?

Now the actual question,
i generated a matrix of 4,4,

   a = [zeros(4,2) ones(4,2)]

% 0     0     1     1
% 0     0     1     1
% 0     0     1     1
% 0     0     1     1

First 4 rows and Two columns are zeros while the remaining are ones as shows above.

Now if i take its FFT2,

    fft2(a)

It generates the following output,

  8.0000            -4.0000 + 4.0000i         0            -4.0000 - 4.0000i
        0                  0                  0                  0
        0                  0                  0                  0
        0                  0                  0                  0

According to the book, the values shown in the output are the DC coefficient and it is indeed the sum of all the matrix values.

In previous example it was easy to calculate the DC coefficient but in this example how these values are generated ?

  8.0000            -4.0000 + 4.0000i         0            -4.0000 - 4.0000i

If i calculate the fft instead of fft2 it will give the following output

%FFT -> DFT of Vector
%FFT2 -> DFT of Matrix
 0     0     4     4
 0     0     0     0
 0     0     0     0
 0     0     0     0

The fft() function calculates the one-dimensional FFT of its input argument. If the input argument is a vector, then the operation is pretty simple to understand; the output is just the result of efficiently calculating a discrete Fourier transform on the input. If the input is a matrix, then, as many MATLAB functions do, each column is treated as a separate one-dimensional signal and are transformed separately. Therefore, the $n$-th column of the output matrix corresponds to the DFT of the $n$-th column of the input matrix.

The fft2() function calculates the two-dimensional FFT of a matrix. One way of envisioning this operation is to take standard one-dimensional FFTs down each column of the input matrix $A$ to yield a temporary result $B$. Then, take one-dimensional FFTs across each row of $B$ to yield the result $C$. The matrix $C$ would then be the result of the two-dimensional FFT of A. 2-D transforms (although not typically the FFT) are often used in image and video compression algorithms.


FFT/PSD/IFFT analysis on single axis piezoelectric accelerometer signals for curb impacts I'm trying to denoise the signal by performing PSD analysis and followed by IFFT. Ultimately, I want to generate Force and Displacement plots from the denoised acceleration signal.
Noisy Acceleration Signal($a_z$ vs t):

PSD analysis of the signal:

Setting a PSD > 0.001 in the code to filter out frequencies having less power than 0.001.
After IFFT($a_z$ vs t):

The denoised signal makes sense since I'm recording z-acceleration on curb impacts which comes out to be a series of impulses.
I'm a novice in signal processing and I don't know whether a windowing would have given a better result or not?
Further questions: Is it possible to find force distribution from the acceleration signal? I've been searching to find answers but none have given me a good idea.

Filtering by zeroing out bins is not a recommended approach as it will introduce significantly more time domain ringing, as detailed here
Why is it a bad idea to filter by zeroing out FFT bins?
Consider using the firls function available in MATLAB / Octave or Python scipy.signal to design optimized multiband filters around your frequencies of interest.
However the frequency content will be driven by repetition in the data; if the actual application will repeat (or not) at unknown and variable intervals, any such filtering techniques will not be useful.


Remove noise on edge (compression deffects) I have the images of cartoons, that were compressed. The example:

They have such noise, that is not easy to delete.

Even though the pixels are on the gray background, the noise pixels can be of very different colors and if I look at the histogram (of grayscale image) of such image I see a lot of bins around the bin of main color. But I can't just delete the colors of value, that are not of main color, because there is a risk to delete some important features (as eye, for example). Also I've tried to posterize image (example presents 8 colors), but some pixels still remain.

Also I've tried median filter, but it is not able to vanish such heavy noise ( I use 3*3 filter).

Please, give some methods that can effectively remove the noise in this case. I would be grateful for any help!

The Total Variation (TV, a.k.a. Rudin-Osher-Fatemi functional) denoising wold help here.
You can apply it on every color layer or look for a vectorized version if th eresults ar enot good enough for you.

Good references of tsable and easy to implememnt algorithms for TV denoising are for example:

Chambolle's algorithm
the Chambolle-Pock primal dual algorithm

Both references have good performance and are easy to implement from the papers.


What do cloud-based speech-to-text services upload to the cloud? Services like Siri and Google Now perform speech recognition (speech-to-text) on mobile phones by doing part of the computation on a server in the cloud.

What do they upload to the cloud?  Do they upload the raw audio?  Or do they do some feature extraction and only upload features (e.g., MFCC coefficients or something) to the cloud?  Is uploading the full audio important/useful for obtaining the best results (e.g., best accuracy of speech recognition)?

I'm curious about understanding the technical details of what is uploaded to Google/Apple, e.g., to better understand the privacy implications.  Note that I'm not asking about their privacy policy, or what they do with the data, or whether I should be concerned about privacy, or any of that stuff.  My interest is in the technical details of how speech recognition computation is split between your phone and their servers.

The research I've done, and what I've found:

Based on public statements that Apple has made, it appears that Siri works by uploading full audio clips (not just features) to Apple's server.  The article also discusses some speculation about whether doing local processing would be feasible.  (See also here.)
Back in 2011, Applidium reverse-engineered the Siri protocol.  They reported that the phone is uploading a full audio clip to Apple's servers: the only computation done on the phone is to compress the audio using Speex.  The scripts they used to help them are available on Github.
I haven't been able to find anything on Google Now.

Its probably safe to assume that all of the major companies send enough information to reconstruct the audio. This is because having that audio for training is such a valuable resource. A certain percentage of the audio segments will be listened to and transcribed by a human annotator.

Also features in these systems are more complicated than MFCCs. You obviously want to offload computing to the cloud so that's another reason to expect that they send the whole audio.


Why DFT is used for approximating CTFT when you can approximate CTFT-integral itself? I was using MATLAB for approximating FTs. Why DFT is used if we can approximate the transform-integration using summation.

Assuming that future vistors won't take the time to read all the comments, I'd like to give a very simple and straightforward interpretation of the discrete Fourier transform (DFT) as an approximation of the continuous-time Fourier transform (CTFT).
First, we need to assume that the signal to be transformed decays sufficiently for large value of $|t|$. (If it doesn't, we can still use the DFT to approximate the CTFT, but the truncation error will be large.) Furthermore, we choose the origin $t=0$ such that the relevant part of the signal appears for $t\ge 0$. Then we can approximate the Fourier integral in the following way:
$$X(j\Omega)=\int_{-\infty}^{\infty}x(t)e^{-j\Omega t}dt\approx \int_{0}^{T}x(t)e^{-j\Omega t}dt\tag{1}$$
with some sufficiently large $T$. The approximation $(1)$ generally results in a truncation error, which can be made small by an appropriate choice of $T$.
Now we can approximate the right-hand side of $(1)$ by a Riemann sum:
$$\int_{0}^{T}x(t)e^{-j\Omega t}dt\approx\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\Omega n\Delta t}\Delta t\tag{2}$$
where $\Delta t$ is the sampling interval and $T=N\Delta t$. The approximation $(2)$ introduces an aliasing error, which can be made small by choosing $\Delta t$ sufficiently small.
Computing the right-hand side of $(2)$ at discrete frequencies $\Omega_k=2\pi k f_s/N$, $k=0,1,\ldots,N-1$, where $f_s=1/\Delta t$ is the sampling frequency,
finally results in
$$\begin{align}X(j\Omega_k)&\approx\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\Omega_k n\Delta t}\Delta t\\&=\frac{1}{f_s}\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\frac{2\pi}{N}nk}\\&=\frac{1}{f_s}\textrm{DFT}\big\{x_d[n]\big\},\qquad x_d[n]=x(n\Delta t)\tag{3}\end{align}$$
Eq. $(3)$ says that the CTFT of $x(t)$ can be approximated on a set of discrete frequencies by the DFT of a sampled version of $x(t)$. Two types of error have to be taken into account: 1. the truncation error, from truncating the signal, and 2. the aliasing error, from sampling the signal.
Note that we can never avoid both errors, because for finite length signals (with zero truncation error) we always get an aliasing error due to the infinite bandwidth of finite length signals, and for perfectly band-limited signals (with zero aliasing error), we always have a truncation error because they can't have finite length. In practice we always have to deal with both types of error.
In sum, answering the question in your question title: under the aforementioned conditions, the DFT is an approximation of the CTFT integral, and it has the advantage that there exist very efficient algorithms (FFT) for its computation.


Difficulty absorbing Idea of interpolation? I am trying to develop my understand of interpolation and signal reconstruction and uptill now i have understood that there are 3 commonly studied types of interpolation
1)Zero order hold interpolation
2)First order hold interpolation
3)Ideal interpolation
As far as i am able to understand,interpolation is processing of filling the empty sample space of sampled signal and zero order and first order tell us how samples are connected to each other.
But still i have two main confusions:
1)What is meant by "order" here in interpolation?Is it referring to derivative??
2)Many figures literature sources that i come across don't explicity contain block of interpolation,although some include block of low pass filter or reconstruction filter as shown in block diagram in attached snap of science direct page. In such scenario,how we can differentiate between low pass filter,reconstruction filter and interpolation filter(zero order hold or first order hold etc)

Order refers to the order of the polynomial that is effectively used in the upsampling process. A zero-order hold will hold the previous value as a constant to insert all new interpolated samples, until the next sample, while a first-order interpolation will insert samples linearly between the existing samples, and so on.
Observe these two polynomial forms below using $m$ as the new higher rate samples and $n$ as the original lower rate samples, such that for $N$ original samples with  $n = 0 \ldots N-1$ there are $M$ higher rate samples with $M=In$ from $m = 0 \ldots M-1$, where $I$ is the upsampling factor:
Zero order hold:
$$y[m] = y[n-1]$$
First order hold:
$$y[m] = y[n-1] + \frac{m-In}{I}(y[n]-y[n-1])   $$
Often interpolation is done by upsampling by inserting zeros between the existing samples to have samples at the higher rate (and to be more complete in the mathematical forms above we would also operate on a zero-inserted signal at the higher rate to then show the two polynomial interpolations given with a consistent index of m), and this is followed by an interpolation filter which "grows" the zeros to the ideal (or close to ideal as limited by the filter design complexity) values for the interpolated signal.
To see how this works, and the additional distortion that a zero-order hold approach would introduce (known as passband droop) please consider how a zero-order hold mathematically can be viewed as the convolution of our original signal with a (sampled) pulse as depicted in the figure below:

Convolution in the time domain is multiplication in the frequency domain. The Fourier Transform of a pulse is a Sinc function (and similarly the Discrete-Time Fourier Transform, DTFT, of a sampled pulse is the Dirichlet Kernel which is an aliased Sinc). This explains the "pass-band droop" distortion that a zero-hold upsampling would introduce as shown in the magnitude of the frequency response for the sampled pulse in the figure below (freqz is a Matlab/Octave/Python function which can produce the DTFT):

In contrast, a zero-insert approach would introduce no distortion within the passband of the signal (but replicates the spectrum in frequency which are removed with the interpolation filter):

As further detailed in other posts, the upsampling process will introduce spectral images that must be filtered out to complete the interpolation process. The ideal interpolation filter (ideal is not realizable but can be approached at the expense of filter complexity and can be designed to comply with any distortion requirements). For example, below shows the new images (in dashed lines) that would be created from an upsampling by 4 (introduce 3 zeros in between each sample), and the desired response of a multi-band interpolation FIR filter (IFIR) that ideally would pass the pass band with no distortion, and completely reject the images.


A special case of 2 jointly Weak-Sense Stationary (WSS) stochastic processes I know that 3 conditions must be met in order a pair of stochastic processes $X(t)$ and $Y(t)$ to be characterized as jointly WSS:

1. $X(t)\;\; WSS$

2. $Y(t)\;\; WSS$

3. $R_{xy}(t_1,t_2) = R_{xy}(t_1 - t_2) = R_{xy}(t_2 - t_1)$ which means that their cross-correlation function should depend only on the difference $T = t_1 - t_2$.

My question is: if the first 2 conditions are met but the cross-correlation function does not depend on any variable (I mean that it's a constant e.g. $R_{xy}(t_1,t_2) = 0$) can we assume that the two stochastic processes $X(t)$ and $Y(t)$ are jointly WSS?

Yes, when $R_{X,Y}(t_1, t_2) = E[X(t_1)Y(t_2)] = 0~ \forall ~t_1, t_2$, $\{X(t)\}$ and $\{Y(t)\}$
are often called uncorrelated WSS processes,
and they are considered to be jointly WSS processes.  More nitpicky folks
agree that the processes are indeed jointly WSS when $R_{X,Y}(t_1, t_2) = 0~ \forall ~t_1, t_2$, but insist on reserving the adjective uncorrelated
for processes for which the cross-covariance function
$$C_{X,Y}(t_1,t_2) = E\left[\left(X(t_1)-E[X(t_1)]\right)\left(Y(t_2)-E[Y(t_2]\right)\right] =E[X(t_1)Y(t_2)]-E[X(t_1)]E[Y(t_2)] = R_{X,Y}(t_1,t_2)-E[X(t_1)]E[Y(t_2)]$$
equals $0$ for all $t_1, t_2$, that is, every pair of random variables $X(t_1)$, $Y(t_2)$ from the processes is a pair of uncorrelated random variables.  Of course, when at least one of the processes has zero mean
(as is usually the case in signal processing applications), then
$C_{X,Y}(t_1,t_2) = R_{X,Y}(t_1,t_2)$ and so the nit gets picked without
requiring any further consideration.


Why is a bandpass filtered analog PM signal still PM and not AM I know, that in theory a sinusoidally phase modulated (PM) signal with an expression like $\Re(e^{\mathrm{j}k\sin(\omega_s t)} \cdot e^{\mathrm{j}\omega_ct})$ with the signal frequency $\omega_s$, the modulation index $k$ and the carrier frequency $\omega_c$ has frequency components at the frequencies $\omega_c \pm n\cdot\omega_s$ with amplitudes calculated from the bessel functions of the first kind.
This means, that for no noise and for an infinite timespan, the fourier transform of such a PM signal consists only of a carrier and (an infinite number of) symmetric spurs around the carrier.
Now when I filter this PM signal using a bandpass with center frequency $\omega_c$ and bandwidth slightly larger than $2\omega_s$, I get a signal with a carrier and a single pair of sidelobes/spurs that are symmetric around the carrier. When I try to demodulate this signal (either using a dsp simulation or using a signal generator and a signal analyzer), I get almost no amplitude modulation and a roughly sinusoidal phase modulation.
I am wondering about two things, I don't understand:

Why is the bandpass filtered signal still a PM signal and not an AM signal? AFAIK, a sinusoidally modulated AM signal consists of the (optional) carrier at $\omega_c$ and two sidelobes at $\omega_c\pm\omega_s$. But this seems to be exactly the spectrum of the filtered PM signal. What am I missing here?
How are AM sidebands/spurs different from PM sidebands/spurs? Given I have a pair of spurs around the carrier at $\omega_c$ with offset $\pm\omega_s$, how can I determine if those spurs are a sinusoidal amplitude modulation or a sinusoidal phase modulation?

Edit 1:
Based on Fat32's answer, I played a little bit with different phase angles between spurs. In particular, I created one-sided frequency vectors with a (single-bin wide) carrier at one frequency and two (single-bin) sidebands with fixed and equal distances left and right to the carrier. The carrier has a mangitude of 1 and an angle of 0. The spurs have a magnitude of 0.5 each and can have independent angles.
If I am not mistaken, a $\cos(\omega t) - \sin(\omega t)$ type of expression should have spurs with phase angles of $\pm \frac\pi2$ relative the carrier (coming from the identity $\sin(\omega t) = -\frac{\mathrm{j}}{2}(e^{\mathrm j \omega t} - e^{-\mathrm j \omega t})$. When I configure my spurs with said angles of $\pm \frac\pi2$ and take the inverse fft of my frequency vector, the resulting analytic signal has zero imaginary component and a sinusoidal magnitude. Which means, it is pure AM modulation.
Did I get the phase angles for the $\cos - \sin$ expression wrong?
And what are the general rules here? When I have a carrier (with an angle of 0 w.l.o.g.) and a pair of symmetric spurs with relative angles $\phi_1$ and $\phi_2$, when exactly is the resulting (analytic) signal only AM or only PM?
Edit 2:
Another attempt of mine was to approach this mathematically. Looking at the analytic baseband signal of the demodulated signal from the first edit, it should have an expression of the form $A_1 + A_2 \left(e^{\mathrm j(\omega_s t + \phi_1)} + e^{\mathrm j(-\omega_s t + \phi_2)}\right)$.
As $A_1$ and $A_2$ are constants, I thought I can simply focus on the sum of the exponentials and determine

For which $\phi_1, \phi_2$ is $\frac{\partial}{\partial t}\left|e^{\mathrm j(\omega_s t + \phi_1)} + e^{\mathrm j(-\omega_s t + \phi_2)}\right| = 0$ (no AM)
For which $\phi_1, \phi_2$ is $\frac{\partial}{\partial t}\angle\left(e^{\mathrm j(\omega_s t + \phi_1)} + e^{\mathrm j(-\omega_s t + \phi_2)}\right) = 0$ (no PM)

Unfortunately, for general angles $\phi_1$ and $\phi_2$, I was not able to derive any general conditions for purely AM or purely PM modulation from a single pair of spurs with equal amplitude and distance and independent angles.

Conventional AM and narrowband PM signals look quite similar:
$$x_{AM}(t) = A \cos(\omega_c t) + m(t) \cos(\omega_c t) $$
$$x_{NBPM}(t) = A \cos(\omega_c t) - A k_p m(t) \sin(\omega_c t) $$
And based on this, their spectrum also look quite similar. However, there's an important distinction between them:
AM modulation changes (modulates) the amplitude of the carrier wave, but keeps carrier phase unaltered. Whereas, PM modulation, NW or WB, keeps the amplitude of the carrier constant, but alters (or modulates) its phase according to the message signal.
Ther reason for this different behaviour is about the way the sidelobes (spurs) are added to the main carrier. In AM, the sidelobes are added in-phase with the carrier (using two cosines for that), while in NBPM they are added in quadrature (sine added to cosine). This quadrature additon, leaves the amplitude (almost) unaltered but changes the phase-angle of the carrier.
So if you treat NBPM as an AM and try to demodulate it that way, you will get constant output (or nothing after the DC blocking), as the amplitude of the NBPM carrier is not changing...


How to check if a signal is power signal or energy signal? What is the procedure to check the signal type?

example:

$ x(t) = A  \sin (\omega t) $

$ y(t) = A e^ {-\lambda |t|} $

let's assume the signal, $x(t)$ is not identically zero for all $t$.
an "energy signal" (what i would prefer to call a "finite energy signal") is such a signal, $x(t)$ with a finite energy:
$$ 0 \ < \ \int_{-\infty}^{\infty} |x(t)|^2 \ dt \ < \ +\infty $$
BTW, sometimes for mathematical ease, we require a stricter sense of finite "energy":
$$ 0 \ < \ \int_{-\infty}^{\infty} |x(t)| \ dt \ < \ +\infty $$
and a "power signal" (what i would prefer to call a "finite power signal") is such a signal, $x(t)$ with finite power:
$$ 0 \ < \ \lim_{T \to +\infty} \frac{1}{T}\int_{-\frac{T}2}^{\frac{T}2} |x(t)|^2 \ dt \ < \ +\infty $$
i think that is the most fundamental definitions of the two classes of continuous-time signals.  you can do a very similar definitions for discrete-time signals, $x[n]$.


How can i record 1 MHz signal and do a real-time processing in a PC? I have a 1 MHz signal from transducer that needs to be analyzed in real-time condition, so the sampling needs to be more than 2 MHz. It consists of digital filtering, turn into spectogram and feed into a model for prediction. I am still stuck at what the device that sufficient for that task. i have tried to use high speed ADC and arduino, but the arduino cannot send the data to the PC rapidly. is there aother method that i could use?

Thanks

Stanley is right and what you need is a data acquisition board. Depending on your application concern, you can either buy it from an online store or build your own if you can do so.

If you would like to make it on your own, then as MBaz suggested you can select from mini pc boards such as a sufficient Raspperry pi that can either transmit your grasped data to your PC at ethernet 100T speeds or the Pi itself can also be used to do the whole analysis as MBaz indicated.

Unfortunately, here you cannot find help on the hardware design of such a board....


Bandlimiting and OFDM spectrum - what is considered "bandlimited"? I have been reading about OFDM modulation, it's not clear to me how you can have band-limited OFDM but given its prevalence in wireless communication there must be some designation of approximately band-limited that OFDM satisfies.

So my question is how is the bandwidth of OFDM systems calculated and if it's an approximation based upon the largest spectrum lobe, why is that good enough?

It doesn't even look like OFDM systems are pulse shaped to reduce the side-lobe strengths so what's going on?

This isn't a question specific to OFDM:

All transmissions need to be bandlimited by law. Hence, law defines band-limited here.

Typically, this happens in the shape of spectrum masks (often expressed in textual form in regulations); these say something like (this is just an example, actual legislation varies across regions and is definitely different):

"not more than -40 dBm/Hz within 10 MHz of the carrier frequency, not more than -90 dBm/Hz between 10 and 20 MHz of the carrier frequency, and -120 dBm/Hz more than 20 MHz away"

Regarding the specialties of OFDM:

You indeed can't pulse shape OFDM without losing its benefits. (You can do a bit of a tradeoff, and trade a bit of orthogonality for out-of-band radiation suppression, but you'd typically avoid that)

That's why you typically have "empty" carriers at the edges. You simply leave room for the sinc to decrease.


Pitch detection of two following music notes (C0 and C0#) Hi signal processing community,

currently I am working on a pitch detection tool which should work on audio input files (or later also instrumental or voice from microphone) to say which music note is playing at the moment.

At the moment I am getting samples in a window of 1024. Then I do FFT-Hamming on this array. Now it is easy to detect notes like C1, C3, C4, C5, since I search for the peaks in the resulting array. But I don`t know how to do it for following notes. (also for example C0, D0)

Is my window to small? Or am I doing something wrong?

I hope you can help me. Thank you.

Best wishes,

MusicMagician

There are many ways to do pitch detection (seee e.g. https://ccrma.stanford.edu/~pdelac/154/m154paper.htm ) and which one is the most suitable depends on your application requirements and constraints. FFT isn't great if you need a lot of accuracy of accuracy (in tune, out of tune?) or good resolution at low frequencies.

You also need to be aware that for many instruments, the harmonic have higher energy than the fundamental, so the peak in the FFT isn't always the pitch.


Finding the pitch value of a note I am trying to implement a karaoke scoring algorithm and, for a particular song, have the note information (done by a professional) that indicates what note occurs in which part of the lyric.

When I use a pitch detector and plot the pitch values for the portion of the lyric where a single note is held, the pitch varies a lot in the beginning (for the first few pitch values), hits a steady value (that roughly corresponds with the note that needs to be sung) and again shows a lot of variation before it transitions to the next note. I want to derive a single pitch value that is representative of the pitch value when the note is held steadily.

How do I go about it? Is this a problem of pitch smoothing (median)?

I am using a pitch algorithm based on the G.729 pitch detector with a pitch smoothing algorithm that seems to work really well (so far) on speech and on male singers singing in the presence of background music (very low volume).

You may need to ignore the transients between stable pitches (for example, most consonants when singing lyrics) as these may be un-pitched sounds or noises, and thus confusing your pitch estimation algorithm.

The singer(s) may also be producing some vibrato, which should correctly show as pitch variations towards the end of some held notes.


Up-converting a sinusoidal signal I want to up-convert a baseband sinusoidal signal of frequency $4000 \, \rm{Hz}$ using a carrier frequency of $20\, \rm{KHz}$. I read that one can up-convert a baseband signal to a passband signal using the following.

But after convolution with a carrier cosine, I get a $16\, \rm{KHz}$ signal. What could be the problem?

I used following Python code to do convolution with carrier frequency: samples * cosineSignal.

A 48k sample rate is too low to reliably reproduce your upper side of 20k+4k = 24kHz.  Thus you end up with a LSB (lower sideband SSB) signal, not an AM signal.


DFT practice question This is probably basic but as I am new to the field it confuses me a bit.
While looking at some solutions provided to a problem in the final step following happens:

$$\begin{aligned} \frac{1}{10}\sum_{l=-\infty}^{\infty}x[l]\sum_{k=0}^9 e^{-i2\pi(n-l)k\frac{1}{10}}
&=\frac{1}{10}\sum_{l=-\infty}^{\infty} x[l]10\delta_{10}[n-l]\\
&=\sum_{s=-\infty}^{\infty}x[n-10s]\end{aligned}$$

Now, according to the DFT table we have provided for the course, the transformation in the 1st step is defined like this

$$e^{2i\pi(k-k_0)/N}\longmapsto N\delta[k-k_0]$$

The delta with subscript was never defined before appearing in this task so searching online I found this:

$$\delta_N[k]=\sum_{s=-\infty}^{\infty}\delta[k-sN]$$
My questions are:
What exactly describes the $\delta_N[k]$?

In the 1st step regarding the delta, is one of the 2 sources wrong or could it be that in the transformation table it should be ''understood'' that the delta has to have the subscript depending on its period?
Is the last step just the result of applying the properties of the delta(with subscript) function on the x signal?

You have to realize that all sequences in DFT relations are inherently periodic. This is usually understood, but it can be made explicit, for instance by that subscript. So both sources are correct, but one explicitly notes the periodicity of the sequence.

If you write down the DFT relation

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

then you see that the expression on the right-hand side of $(1)$ is periodic in $n$ with period $N$. So $\delta[n]$ on the left-hand side is actually continued periodically. If we just consider one period, we don't care about that periodic continuation. However, if we are concerned about the behavior for all $n$, then we do need to take that periodicity into account, and we could have written $\delta_N[n]$ on the left-hand side of $(1)$ to be more explicit.

And in the given example we do care what happens outside that period because the index $l$ in the first sum extends from $-\infty$ to $\infty$.

The last step in the given equation is just a convolution of the sequence $x[n]$ with a unit impulse comb, so $x[n]$ is shifted by integer multiples of $10$ and all those shifted versions are summed up (hoping that the sum converges).


Linear phase FIR filter for impulse responses that don't appear symmetric I would like to clarify some confusion I have about linear phase FIR filters of which do not seem to have symmetric impulse responses.

Starting with a simple case, a delay, $ h[n] =\delta(n-n_0) $ does not seem to have  have a symmetric impulse response. For example,

$$h[n] = [0, 0 ,0,1,0] $$ is not a symmetric impulse response. How is it that the property of symmetry for a linear phase FIR filter holds true, in this case?

For another example, say I have a filter defined as

$$h[n] = [0,0,0,0,1,0,0,0,0,0,1]$$

I believe this filter is linear phase too. But it's impulse response is not symmetric.

How does the symmetric property hold true? Does it have to do with zero padding? But even if we don't zeropad the filter, it is still linear phase? If we zeropad it by a large factor, is it also linear phase?

For one more example, say I have this filter:

$$ h[n] = [1, 2, 3, 4, 3, 2, 1]$$

This is obviously linear phase. When I zeropad it by a factor N, the impulse response is no longer symmetric, but it stays linear phase for all values of N. Why is this the case?

Your confusion is understandable.

If you consider the definition of linear phase FIR filter and the associated symmetry conditions on their impulse responses, then you can arrive the conclusion that the first two cases

$$ h_1[n] = [0,0,0,1,0] $$

and

$$ h_2[n] = [0,0,0,0,1,0,0,0,0,0,1] $$

are non-symmetric. However, as you use zeros and ones in those impulse responses, it can be seen that the following two new augmented impulse responses are equivalent to those non-symmetric looking ones:

$$ h_3[n] = [0,0,0,1,0,0,0] $$

and

$$ h_4[n] = [0,0,0,0,1,0,0,0,0,0,1,0,0,0,0] $$

And they are symmetric and linear phase according to the definition as well. This happens because of the zeros and use of ones in a specific way.

However, the following impulse response

$$h_5[n] = [0,0,1,0,0,0,0,1,0,1,0]$$

is not symmetric, is not linear phase and cannot be made linear phase by augmenting with zeros...


Noise floor of a FM duplex If these 2 signals, that form this duplex, have ranges of 15-30 [kHz] and 30-45 [kHz] using:

$NoiseFloor_{dBm} =10\log _{10}(k\times T_{0}\times 1000)+NF+10\log _{10}(BW)$ (formula from wikipedia)

and if $k,T,NF$ and $\Delta f$ are given, only thing left to figure out is the $BW$ (bandwidth). Is the $ BW=2(\Delta f+f_{m})$ where $f_{m}=2\times 15=30 [kHz]$ correct formula for the bandwidth of this duplex?

Yes, the 30 kHz are correct, but only assuming your receiver has filters that filter out any noise outside these bands. That might or might not be the case – if your receiver in fact observes 0 – 45 kHz, for example, then the noise energy it'll see is defined by that 45 kHz bandwidth.

I'd also recommend you not just peruse the formula from Wikipedia (which I'd argue is in terrible shape), but try to understand how this noise floor calculation comes to be, and then, after that, apply the conversion to dBm.

Namely, you get Johnson-Nyquist thermal noise. This noise has a power density of $k_B T$, that is, Boltzmann's constant times the temperature of your system. This constant represents Energy per temperature above absolute zero (eg. Joules per Kelvin), so $k_B T$ is Energy.

In signal-processing systems we commonly care about Power more than energy (because, for example, if we know how much power gets dissipated by a resistor, we know the voltage across it), so we divide the Energy by time, or, multiply with bandwidth/frequency (which is the inverse of time):

$$P_\text{N, thermal}=k_B T \Delta f$$

Real-world active systems, especially receiver's amplifiers, add more than just the noise present due to thermal noise – the factor between the noise going in and coming out is the Noise Figure $\text{NF}$. Therefore:

$$P_\text{N, thermal}=k_B T \Delta f\, \text{NF}$$

Because we like our $\text{dB}$-based units, let's cancel out all units (I'm doing this in-line here just for clarity) and then do the dB-conversion. Then we get

$$P_\text{N, thermal}[\text{dBW}]=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)(\Delta f \cdot s) \, \text{NF}\right)$$

You wanted $\text{dBm}$, but a milliwatt is simply $\frac1{1000}$ of a Watt, so

$$\begin{align}
P_\text{N, thermal}[\text{dBm}]=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)(\Delta f \cdot s) \, \text{NF}\,\cdot 1000\right)
\end{align}$$

Let's apply a few logarithm rules here:

$$\begin{align}
P_\text{N, thermal}[\text{dBm}]&=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)\right)&+&10 \log_{10}\left(\Delta f \cdot s \right)&+&10 \log_{10}\left( \text{NF}\right) &+& 10\log_{10}\left(1000\right)\\
&=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)\right)&+&10 \log_{10}\left(\Delta f \cdot s \right)&+&\text{NF}_\text{dB} &+& 30
\end{align}$$

If you drop in the values of Boltzmann's constant and room temperature, you'll end up with the classical -174dBm/Hz; simply add your bandwidth (in dBHz) and your Noise Figure (in dB), and you get your observed noise power.


Zero Frequency Basis Equation for Discrete Sine Transform So, I'm evaluating the Discrete Sine transform Basis functions and I found on the internet the following picture:

So, I may be getting the point wrong, but, if the 2D DCT equation is:

Then the Discrete Sine transform shall be:

Then, it is clear that for the lower frequency (u=0,v=0) all cosines will be 1, but in the case of the sines all will be 0 (for frequency=0), therefore, the basis function for the lower frequency, in the sine case, will be zero, right?

Is the upper left image (u=0,v=0) correct? It should be constant, right? (btw: It shall also be zero, right?)

I think I got it:

While in the DCT there are definitions of the DCT that make cos(0) when u and v are zero, this does not apply to DST (see: https://en.wikipedia.org/wiki/Discrete_sine_transform).

In the DST definition, there is always a +1 or a +1/2 adding to the frequency.

This explains why it is not zero and also the graph above. It varies
(for u=0 and v=0) because the content of the sine is not zero.


Why does Fourier zero-fill interpolation pad the middle frequencies? In this: http://dspguru.com/dsp/howtos/how-to-interpolate-in-time-domain-by-zero-padding-in-frequency-domain
the author says that we can interpolate a signal by zero-padding the middle frequencies in the DFT. But then the high frequency components are pushed even higher, so shouldn't the result look more wiggly? Why is that not so?

Usually a fast Fourier transform (FFT) returns a vector that starts with the zero-frequency bin. This means that the Nyquist frequency (sampling frequency / 2) bin is in the middle of the vector. So you are not pushing apart the low and high frequencies by zero padding at the middle. The zeros go between the highest positive frequency bin and the highest-absolute-frequency negative frequency bin.

Rick omits the details of how to deal with the Nyquist frequency bin (which is there with even-length DFT): If you have a non-zero Nyquist frequency bin then you need to split it into a half-amplitude positive Nyquist frequency bin and a half-amplitude negative Nyquist frequency bin, and insert the zero padding between those. ...I think.


Help with equations for exponential ADSR envelope With application code, I have implemented a linear ADSR envelope for shaping the amplitude of an oscillator's output.  The parameters for attack, decay and release duration as well as sustain level can be set on the envelope and everything works as expected.

However, I would like to tweak the ramp shapes of the envelope to something that resembles what most synthesizers use for a more natural response: inverse exponential for the attack and exponential for the decay and release.  I am having trouble getting my formulas right for calculating the envelope output values for these types of ramp shapes.  To calculate the linear ramps, I am using the two-point form, plugging in the start/end $x$/$y$ values which are derived from the attack/decay/sustain/release input parameter values.  I cannot seem to work out the correct formula for exponential (standard and inverse) ramps using the same start/end $x$/$y$ point values.

I have saved a Desmos Graphing Calculator session that demonstrates the approach for linear ramps that I described above.

If anyone can help point me in the right direction, it would be much appreciated.

I think that what confuses you is that a decreasing exponential ($e^{-x}$) does never reach 0, so an ADSR generator with truly exponential segments would stay stuck ; because it would never reach the target value. For example, if the generator is at the height of the attack phase (say $y = 1$) and has to land to a sustain value at $y = 0.5$, it can't go there with a true exponential, because the true exponential won't decay to 0.5, it'll only asymptotically go to 0.5!

If you look at an analog envelope generator (for example the 7555 based circuit everybody seems to use), you can see that during the attack phase, when the capacitor is charging, it is "aiming higher" than the threshold used to indicate the end of the attack phase. On a (7)555 based circuit powered by +15V, During the attack stage, the capacitor is charged with a +15V step, but the attack stage ends when a threshold of +10V has been reached. This is a design choice, though 2/3 is the "magic number" found in many classic envelope generator, and this might be the one musicians are familiar with.

Thus, the functions you might want to deal with are not exponentials, but shifted/truncated/scales versions of it, and you'll have to make some choices as to how "squashed" you want them to be.

I am curious anyway as to why you are trying to get such formulas - maybe it's because of the limits of the tool you are using for synthesis ; but if you are trying to implement those using a general purpose programming language (C, java, python) with some code running for each sample of the envelope, and a notion of "state", read on... Because it's always easier to express things as "such segment will go from whatever value it has just reached to 0".

My two pieces of advice on implementing envelopes.

The first one is not to try to scale all the slopes/increments so that the envelope exactly reach start and end values. For example you want an envelope that goes from 0.8 to 0.2 in 2 seconds, so you might be tempted to compute an increment of -0.3 / second. Don't do that. Instead, break it down into two steps: getting a ramp that goes from 0 to 1.0 in 2 seconds ; and then applying a linear transform that maps 0 to 0.8 and 1.0 to 0.2. There are two advantages to work this way - the first is that it simplifies any computation you'll have relative to envelope times to a ramp from 0 to 1 ; the second is that if you change the envelope parameters (increments and start/end times) midway everything will remain well-behaved. Good if you're working on a synth, since people will ask to have envelope time parameters as modulation destinations.

The second is to use pre-computed lookup table with envelope shapes. It is computationally lighter, it takes out many dirty details (for example you don't have to bother with an exponential not reaching 0 exactly - truncate it at your whim and rescale it so that it is mapped to [0, 1]), and it's dead easy to provide an option to alter envelope shapes, for each stage.

Here is the pseudo-code for the approach I describe.

render:
  counter += increment[stage]
  if counter > 1.0:
    stage = stage + 1
    start_value = value
    counter = 0
  position = interpolated_lookup(envelope_shape[stage], counter)
  value = start_value + (target_level[stage] - start_value) * position

trigger(state):
  if state = ON:
    stage = ATTACK
    value = 0  # for mono-style envelopes that are reset to 0 on new notes
    counter = 0
  else:
    counter = 0
    stage = RELEASE

initialization:
  target_level[ATTACK] = 1.0
  target_level[RELEASE] = 0.0
  target_level[END_OF_RELEASE] = 0.0
  increment[SUSTAIN] = 0.0
  increment[END_OF_RELEASE] = 0.0

configuration:
  increment[ATTACK] = ...
  increment[DECAY] = ...
  target_level[DECAY] = target_level[SUSTAIN] = ...
  increment[RELEASE] = ...
  envelope_shape[ATTACK] = lookup_table_exponential
  envelope_shape[DECAY] = lookup_table_exponential
  envelope_shape[RELEASE] = lookup_table_exponential


Linear Regression Filter Properties I have a filter to detect when a slow moving signal (oven temperature) reaches steady state. I am using a linear regression of a moving window and looking at the slope. Based on the equations for linear regression including $\sum xy$, $\sum x^2$ and $(\sum x)^2$ I am thinking this is a non linear filter, is this true? Is it still true to say this is a FIR filter though?

Also if I am trying to get an estimate for reaching steady state is there a better filter, or can similar quality be achieved with a linear filter given Gaussian noise? Basically I am looking for a good theoretical analysis of this type of filter to way its pros and cons?

This is a linear filter because the signal values ($y_i$ in your notation) enter the equation linearly. It can be implemented as an FIR filter, and then it's called "Slope Filtering". The method was described in the IEEE Magazine on Signal Processing ("DSP Tips & Tricks"), Nov. 2008. You can read the article here.


Concepts / Algorithms for recognising words in the Time-Domain I'm currently working on a project that involves in Speech and recognising Phonemes (in kind of real-time).

The current algorithm works by always listening for input from the Mic, and input that is above a specific threshold then this is recorded to a .wav file and sent to the server for processing. E.g.

listening...
listening...
listening.. (max > threshold)
recording..
recording..
etc..

What I want to do is therefore have it recording when a specific name has been said. Someone basically speaks to it.. I.e. (let's assume it's called iris) so the algorithm would therefore work like the following:

listening...
listening..
listening.. (max > threshold)
  // determine if the input is "iris"
  // if true:
     begin recording the sentence
  // else:
listening..

I want to do this in the time-domain, then, process the sentences in the frequency domain. I've read an article, and, a paper that describes this process using Total number of Zero-crossings (using a specific threshold) my questions is whether or not this would be a good measure? I only need to identify one word, that is "Iris" and then everything else can be done using algorithms such as: MFCC, HMM's etc..

Any ideas? :)

I highly doubt that something as complex as recognizing a word can be done directly in the time domain using features as basic as zero-crossing rate.


Why is aliasing inherently non-linear? E.g. when having 2 signals $x(t)$ and $y(t)$ and transforming each one with adding their spectra, the operation is linear, as the result would be the same as the transformation of $(x + y)(t)$.

Even when looking at the sampling series, each of the infinite elements adds up linearly.

But how is aliasing a non-linear operation and how is it proven mathematically?

I think you mean "images", not "aliases".  They become aliases if there is foldover from resampling.

It's because you are not adding two signals, $x(t)$ and $\operatorname{III}(t)$, you are multiplying them that these images appear.

$$\begin{align}
x_\text{s}(t) & \triangleq x(t) \cdot \operatorname{III}(t/T)  \\
&= x(t) \cdot \sum\limits_{n=-\infty}^{+\infty} T \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x(t) \cdot \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x(nT) \cdot \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \delta(t-nT) \\
\end{align}$$

where $x[n] \triangleq x(nT)$ and

$$ \operatorname{III}(u) \triangleq \sum\limits_{n=-\infty}^{+\infty} \delta(u-n) $$

$$\begin{align}
\operatorname{III}(t/T) &= \sum\limits_{n=-\infty}^{+\infty} \delta\left(\tfrac{t}{T}-n\right) \\
&= \sum\limits_{n=-\infty}^{+\infty} \delta\left(\tfrac{t-nT}{T}\right) \\
&= \sum\limits_{n=-\infty}^{+\infty} T \delta(t-nT) \\
&= \sum\limits_{k=-\infty}^{+\infty} e^{j k \frac{2 \pi}{T} t} \\
\end{align}$$

The last line is from doing fourier series.  Now, if you use the shifting property of the Fourier Transform, then the Fourier Transform of $x_\text{s}(t)$ is

$$\begin{align}
X_\text{s}(f) & \triangleq \mathscr{F}\{x_\text{s}(t)\} \\
&= \mathscr{F}\{x(t) \operatorname{III}(t/T) \} \\
&= \mathscr{F}\left\{x(t) \sum\limits_{k=-\infty}^{+\infty} e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \mathscr{F}\left\{\sum\limits_{k=-\infty}^{+\infty} x(t) \, e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \sum\limits_{k=-\infty}^{+\infty} \mathscr{F}\left\{ x(t) \, e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \sum\limits_{k=-\infty}^{+\infty}  X\left(f-\tfrac{k}{T}\right) \\
\end{align}$$

where $ X(f) \triangleq \mathscr{F}\{ x(t) \} $ .

That non-linear process of multiplication generates frequency components than did not previously exist in $x(t)$.  Those new components are simply shifted versions of $X(f)$ and are called "images".


Generalizing the chain rule Given the real-valued functions $f_1$ and $f_2$ with $x\in\mathbb{R}$, then
$$ \frac{df_2(f_1(x))}{dx} = \frac{df_2(f_1(x))}{d f_1(x)}\frac{df_1(x)}{dx}$$

Is it then the case that if we also have a real-valued function $f_3$, then
$$ \frac{df_3(f_2(f_1(x)))}{dx} = \frac{df_3(f_2(f_1(x)))}{df_2(f_1(x))}\frac{df_2(f_1(x))}{df_1(x)}\frac{df_1(x)}{dx}$$

Is this correct?

If so, does the chain rule generalize to the following statement:
$$ \frac{df_{N}(f_{N-1}(...f_2(f_1(x))))}{dx} $$
$$= \frac{df_{N}(f_{N-1}(...f_2(f_1(x))))}{df_{N-1}(...f_2(f_1(x)))}    \frac{df_{N-1}(...f_2(f_1(x)))}{df_{N-2}(...f_2(f_1(x)))}...\frac{df_2(f_1(x))}{df_1(x)}\frac{df_1(x)}{dx}$$

My apologies about this looking like a straight calculus question. I am asking to figure out the backpropagation algorithm in its general form.

I'm not a mathematician but if the multiple composition $ f_1(f_2(f_3 \dots)))$
is equivalent to a sequence of single compositions $g_2(g_1)$, and another $g_3(g_2)$ and so forth.

The generalization looks reasonable for  $x, f_{1}, f_{2} \dots f_N \in \cal{R}$ and $f_i$ differentiable with no discontinuities.