Why are analytic signals so important in Time-frequency analysis? I am little confused about why we need analytic signals so bad in time-frequency analysis. What might happen if I use non-analytic signals to do time-frequency analysis?

Assuming time-frequency aims a providing a separation (at least visual) between signal components, the main reasons could be:

for quadratic distributions, which tend to yield interference  between components, "cancelling" at negative frequencies reduce the quantity of components that can interfere.
for linear distributions, the filter bank formalism, and especially the down-sampling operators, is simplified, reducing the impact of aliasing errors.

For real signals, the Hermitian symmetry yield that "no information" is lost in the analytic form, a complex-valued function that has no negative frequency components, and it is easy to go back to real.

However, this is not so simple in practice.

If the analytic signal constructed from a wide-sense stationary (WSS) real signal is always proper, this may not be the case for non-stationary signals, as underlined in Stochastic time-frequency analysis using the analytic signal: why the complementary distribution matters.
Moreover, computing the analytic signal on discrete, finite-length data is often only approximate, especially for real-time applications. Hence the energy in the negative frequencies is rarely zero.
Extending the concept of analyticity beyond 1D is not evident, and several design still exist.


Dependence of correlation peak height on GPS C/A signal's SNR, noise floor level (dBFS), correlation period length and ADC resolution Assuming correlation at 2X the chip rate (2046 chips per period) and 500 Hz frequency search bins.

There are probably additional variables here that I'm missing and I'm not sure how to quantify a good peak but.

I'm trying to determine how much amplification I need in my acquisition setup to get good correlation peaks when correlating over 1 C/A code period.

I'm using a ceramic patch + 33dB gain GPS antenna with a 3m RG174 coax + RTL-SDR v3 (7 bit ADC) with 49dB gain (I'm not sure both are probably actually less gain, I've no way to measure).
The original data was sampled at 2.048MSps complex baseband.

I'm correlating in Octave with the signal re-sampled to 2.046MSps and using circular cross-correlation with the C/A code extended to be 2046 chips per period.

I'm only getting good correlation peaks over at least 5 periods for the strongest of signals and 10 or 20 periods for weaker ones. I also had to use 100Hz frequency bins, 500Hz was not enough to get good peaks.

I think my Octave script is fine, I've cross checked it with this script https://www.mathworks.com/matlabcentral/fileexchange/26035-joint-frequency-and-delay-correlation (feeding it my data) and the results are similar, I see a weak peak for the strongest signal in the recording.

Even if your amplified signal matched the quantization noise floor of the ADC (representing a 3 dB degradation) you would still achieve a significantly larger and clearer peak from what you are showing.

The search bin size for 1 C/A code duration of 500 Hz is a good choice for the largest step (since the correlation versus frequency offset would follow a Sinc function response with the first null at $1/T$ where here $T$ is 1 ms, so the correlation would go to zero with a 1 KHz offset from true frequency.

With a 1 C/A code duration you can expect a processing gain of 30 dB ($10Log_{10}(1023)$) and the unspread GPS signal itself can be anywhere from -23 dB to -15 dB below the noise floor typically within the 1 MHz bandwidth of the main lobe of the GPS Signal. (You can measure your post correlation SNR by to see if it is in the ballpark of these levels plus the processing gain).

As for gain, the quantization noise of the ADC can be approximated by being 6 dB/bit + 1.76 dB below the full scale level of a sine wave, and that quantization noise is the total noise power spread evenly across the Nyquist bandwidth. I typically amplify the signal so it is at least 10 dB above this noise level at the input to the ADC, resulting in a 0.4 dB noise figure from the quantization noise. The 0.4 dB quantity comes from the power sum of the two noises (amplified thermal noise as the received signal and quantization noise): $10Log_{10}(1+10^{d/10})$ where d is the delta in dB that the received signal is above the quantization noise.

So if you ADC starts clipping with a sine wave that is at 10 dBm (for example), then for a 7 bit ADC, the quantization noise will be at 6dB x 7 bits + 1.76 = 43.8 dB lower or -33.8 dBm, spread across +/-1.023 MHz in your case with a sampling rate of 2.046 MHz. So within the 1 MHz main lobe of the GPS signal, the quantization noise will be 3 dB lower (since 1 MHz is half of this ADC bandwdith), so that puts the quantization noise in band at -36.8 dBm.

The received signal assuming a front-end noise figure of 3 dB is -174 dBm/Hz (thermal noise floor) + 60 dB (1 MHz bandwidth) + 3 dB (noise figure) = -111 dBm (these are typical numbers, your mileage may vary).

So if we want to place the recieved signal to be 10 dB above the ADC quantization noise floor in order to minimize further SNR degradation while maintain as high of a dynamic range as possible (the typical trade space and why we don't want to apply too much gain--we get little in SNR improvement but lose our ability to handle stronger signals or interference) we would need to have the received signal that starts at -111 dBm amplified to be -36.8 dBm (quantization noise floor) + 10 dB = -26.8 dBm.

With all those assumptions, your recommended gain is 84 dB. This hopefully shows you how to determine it for your specific case.


Expressing DTFTs in terms of one another with similar time domain signals I cam across a question in my DSP book asking this:
1). Express $X_2(e^{j\omega})$ in terms of $X_1(e^j\omega)$ without explicitly computing $X_1(e^{j\omega})$. ($X_1(e^j\omega)$ represents the DTFT of $x_1[n]$)
\begin{align}
x_1[n] &= \left[1 \ 2 \ 2 \ 1\right]\\
x_2[n] &= \begin{cases}
x_1[n] & 0 \leq n \leq 3\\
x_1[n-4] & 4 \leq n \leq 7\\
0   & \text{otherwise}\end{cases}
\end{align}
Here is my code and answer, but Im dont think my answer is correct. What should I be looking for in the time domain graphs that can give clues to the DTFT? ( i know my dtft function needs some work but I will do that next)
Thanks
%% P3.2
% 1). The magnitude plot of X2(jw) should have a similar shape to X1(jw) because
%     they contain the similar frequency components, but part of X2 is just phase
%     shifted, therefore the angle plots might be a bit different.
% 2).

n1 = 0:3;
x1 = [1 2 2 1];
stem(n1,x1); axis([0 7 -4 4]);

n2 = 0:7;
x2 = [1 2 2 1 1 2 2 1];
%    ( x1[n])(x1[n-4])
figure;
stem(n2,x2); axis([0 7 -4 4]);

k = 0:200;

X1 = dtft(n1,x1,k);

X2 = dtft(n2,x2,k);

function X = dtft(n,x,k)
% DESCRIPTION
%  Performs an approximate Discrete Time Fourier Transform
%  the approximation is due to the fact that the output is
%  not continous. ex. X = dtft(n,x,k);
%
% NOTE : Only for finite duration signals.
%
% INPUT VARIABLES
%  n = number of samples
%  x = array containing time domain signal samples
%  k = evenly spaced divisions of omega (frequency)
%
% OUTPUT VARIABLES
%  X = the transformed output
%  Various graphs of signal aspects
%
% REFERENCES
%  Adapted from "Digital Signal Processing Using MATLAB 3rd ed." - Ingle, V.,
%  Proakis, J. (pg. 65).
%
% DOCUMENTATION
%  ver 1.0 by Dominic Meads  5/8/2020
%  filename: dtft.m
%
% ENGINEER'S COMMENTS
%  The authors reccommend to use this more as an
%  excersize rather than a full function. As stated above, use the DFT for
%  better results.
%

w = (pi/100)*k;  % calculates the evenly spaced frequencies
X = x * (exp(-j*pi/100)) .^(n'*k); % calculates the DTFT

magX = abs(X);  % for graphs
angX = angle(X);
realX = real(X);  % divide into real and imaginary parts
imagX = imag(X);

figure('Color', [1 1 1]);
subplot(2,2,1); plot(w/pi,magX); grid off;
xlabel('frequency in units of pi'); title('Magnitude part');
subplot(2,2,2); plot(w/pi,realX); grid off;
xlabel('frequency in units of pi'); title('Real part');
subplot(2,2,3); plot(w/pi,angX); grid off;
xlabel('frequency in units of pi'); title('Angle part');
subplot(2,2,4); plot(w/pi,imagX); grid off;
xlabel('frequency in units of pi'); title('Imaginary part');
end

The problem is asking you to find an analytical expression given $x_1[n]$ and $x_2[n]$. There is no need to compute anything, let alone code it up!
The problem wants you to identify that given the discrete-time domain functions, you use properties of the DTFT to write $X_2(e^{j\omega})$ in terms of $X_1(e^{j\omega})$. You've already identified that $x_2[n]$ contains a shifted version of $x_1[n]$.
You can use the time-shift property which says that if the DTFT of $x[n]$ is $X(e^{j\omega})$ then
$$\mathcal{F}(x[n - m]) = e^{-j\omega m}X(e^{j\omega})$$
Keeping in mind how the domain of $x_2[n]$ is defined, you can use the DTFT properties to express $X_2(e^{j\omega})$ in terms of $X_1(e^{j\omega})$.


How to do the Fourier Transform of bounded function? I was trying to solve a Fourier transform of a function using the properties of Fourier transforms. The function is given as:

$\frac{At}{2}$ for $-2<t<2$ and $0$ for all other $t$. Doing the integral from the definition of of a Fourier transform gives the result: -$Ai\left(\frac{\sin(2\omega)-2\omega\cos(2\omega)}{\omega^2}\right)$. While when I try to do it with derivative property of Fourier transforms I get a different answer as $\frac{2A\sin(2\omega)}{2i\omega^2}$.

The idea is that I express the above bounded function as the product of the function with a rectangular pulse that has a period of 4. Like this: $\frac{At}{2}p(t)$. Then I take the derivative of this as just $\frac{A}{2}p(t)$ and then take the Fourier Transform of $\frac{A}{2}p(t)$. For the Fourier transform of the main function I simply divide by $iw$.

The two results are completely different from each other. What's wrong with taking my approach?

P.S I am sorry for not being able to write the formulas in a very mathematical manner.

Well, imho, one good way to do it is to take the signal's derivative in time and make use of the property $$\frac{d}{dt}x(t) \longleftrightarrow j2\pi f X(f)$$

I'll try to show it analytically.

Since the signal is discontinuous at $t=\pm 2$ you should include Delta function in the derivative, that is, $$\frac{d}{dt}\Big[\frac{A}{2}t\:\mathrm{rect}\Big(\frac{t}{4}\Big)\Big]= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big)$$

Given that a unitary rectangular pulse can be written as a sum of two unit step functions $$\mathrm{rect}\Big(\frac{t}{T}\Big) = u\Big(t+\frac{T}{2}\Big) - u\Big(t-\frac{T}{2}\Big)$$ then the pulse's derivative is $$u'\Big(t+\frac{T}{2}\Big) - u'\Big(t-\frac{T}{2}\Big) = \delta\Big(t+\frac{T}{2}\Big) - \delta\Big(t-\frac{T}{2}\Big)$$

Going back, we get $$\frac{d}{dt}x(t)= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big) = \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\Big(\delta(t+2) - \delta(t-2)\Big)$$

The sampling property of a Delta function yields

$$t\delta(t \pm t_0) = t\Big|_{t=\pm t_0}\delta(t \pm t_0) = (\mp t_0)\delta(t\pm t_0)$$

and in your case

$$\frac{A}{2}t\Big(\delta(t+2) - \delta(t-2)\Big)=-A\delta(t + 2)-A\delta(t - 2) $$

and in total

$$\frac{d}{dt}x(t)= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big) = \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) -A\delta(t + 2)-A\delta(t - 2)$$

The Fourier Transform of the latter is

$$F\Big\{\frac{d}{dt}x(t)\Big\} = 2A\mathrm{sinc}(4f) - Ae^{j2\pi 2f} - Ae^{-j2\pi 2f} = 2A\mathrm{sinc}(4f) - 2A\cos(2\pi 2f)$$

and from the derivative property $$X(f) = \frac{1}{j2\pi f}\Big(2A\mathrm{sinc}(4f) - 2A\cos(2\pi 2f)\Big) = \frac{jA}{2\pi f}\Big(2\cos(2\pi 2f) - 2\mathrm{sinc}(4f) \Big)$$

which gives

$$X(f) = jA\Big(\frac{2\cos(2\pi 2f)}{2\pi f} - \frac{\sin(4\pi f)}{(2\pi f)^2} \Big)$$

Finally, with the $j\omega$ notation

$$X(j\omega) = jA\Big(\frac{2\omega\cos(2\omega)}{\omega^2} - \frac{\sin(2\omega)}{\omega^2} \Big) =jA\Big(\frac{2\omega\cos(2\omega) - \sin(2\omega)}{\omega^2} \Big)$$


matched filter for doppler To remove the Doppler shift, we can use a bank of matched filter. But, to use this method, we need an interval $[-f_{max};f_{max}]$.

How can we approximate this interval?

you sample the interval and build a filter for each frequency shift, unless you don’t need to, like in typical carrier/phase recovery schemes.

You make a bank of filters, particularly if Doppler is of itself, of interest.  In SONAR, tones have better Doppler resolution, and the DFT provide a bank of filters.

The ambiguity function function for the waveform will show how much loss you can tradeoff for the frequency sampling you choose.

Each waveform will have different loss for frequency miss-match.  The ambiguity function also accounts for delay resolution.

A waveform like hyperbolic FM is robust to Doppler so a single waveform covers a wide range of possible shifts.

Both miss-matches can couple.

In comms, Doppler compensation is usually a matter of compensation of the frequency shift, not full time dilation.  Shifting a matched filter is a matter of frequency shifting and can be accommodated in carrier/ phase recovery.


Find time lag between two time series I am trying to find the time lag between two time series over t = [0,1000] using MATLAB (not that it matters). The first time series is simply t^2. The second is (t-15)^2 which is, of course, shifted to the right 15 units (e.g., seconds). My approach has been to find the cross correlation (computed using FFT) and then use the maximum of these to determine the appropriate time shift. For some reason, the answer I keep getting is 0, which doesn't make any sense. Can anyone tell me what I'm doing wrong? Thanks!

I've calculated the cross correlation using FFT as shown (all code is in MATLAB; note fft is the Fast Fourier Transform function and ifft is the inverse Fast Fourier Transform function):

t = [1:1000];
a = t.^2;
b = (t-15).^2;

if (length(a) < length(b))
  c = [zeros(1,length(b)-1) a zeros(1,length(b)-length(a))];
  d = [b zeros(1,length(b)-1)];
else
  c = [zeros(1,length(a)-1) a];
  d = [b zeros(1,length(a)-length(b)+length(a)-1)];
end

c_fft = fft(c);
d_fft = fft(d);

cross_spectral_density = c_fft.*conj(d_fft);
cross_correlations = ifft(cross_spectral_density);

[max_xcorr shift_idx] = max(abs(cross_correlations));
shift = shift_idx - max(length(a),length(b));

I have no good explanation for why this happens, except that perhaps quadratic signals are nonstationary in that their mean is changing (or unbounded).

If I do something to ensure the mean is bounded (and not changing) such as:

a = sin(t.^2);
b = sin((t-15).^2);

or

a = rem(a,100);
b = rem(b,100);

then I get the right answer: 15 samples.

This is what happens when I use the first:

This is what happens when I use the second:

both modified signals generate the same (correct) estimate of delay.

Scilab code to generate it below.  I also modified your zero padding etc. because it seems overly complicated.

Code Only Below

// 25867
t = [1:1000];
a = (t.^2);
b = ((t-15).^2);
a = a - fix(a/200)*200;
b = b - fix(b/ 200)*200;

c = [a zeros(1,length(b) - 1)];
d = [b zeros(1,length(a) - 1)];

mu = mean([mean(c) mean(d)])

c_fft = fft(c - mu);
d_fft = fft(d - mu);

cross_spectral_density = c_fft.*conj(d_fft);
cross_correlations = fftshift(ifft(cross_spectral_density));
xc = xcorr(c- mu,d- mu,length(t)-1);

[max_xcorr shift_idx] = max(abs(cross_correlations));
shift = shift_idx - max(length(a),length(b));

clf;
subplot(211);
plot(a);
plot(b,'r');
subplot(212)
plot(cross_correlations);
plot(xc,'g:')
plot(shift_idx,cross_correlations(shift_idx),'r.')
xtitle('Delay is ' + string(shift))


DFT with geometrically-spaced bins? The traditional Discrete Fourier Transform (DFT) and its cousin, the FFT, produce bins that are spaced equally. In other words, you get something like the first 10 hertz in the first bin, 10.1 through 20 in the second, etc. However, I need something a little different. I want the range of frequency covered by each bin to increase geometrically. Suppose I select a multiplier of 1.5. Then we have 0 through 10 in the first bin, I want 11 through 25 in the second bin, 26 through 48 in the third, etc. Is it possible to modify the DFT algorithm to behave in this fashion?

To quote my dissertation:

  A collection of transforms are given the name constant Q and are
  similar to the Fourier transform.

  Computation of the discrete Fourier transform can be very efficient
  when employing the use of the fast Fourier transform. However we
  notice that energy of a signal is divided into uniformly sized
  frequency buckets across the spectrum. While in many cases this is
  useful, we notice situations where this uniform distribution is
  sub-optimal. An important example of such a case is observed with the
  analysis of musical frequencies. In Western music, the frequencies
  that make up the musical scales are geometrically spaced. We therefore
  see that the map between frequency bins of the discrete Fourier
  transform and the frequencies of musical scales is insufficient in the
  sense that the bins match poorly. The constant Q transform addresses
  this issue.

  The aim of the constant Q is to produce a set of logarithmically
  spaced frequency bins in which the width of the frequency bin is a
  product of the previous. As a result we may produce an identical
  number of bins per musical note across the audible spectrum, thus
  maintaining a constant level of accuracy for each musical note. The
  frequency bins become wider towards the higher frequencies and
  narrower towards the lower frequencies. This spread in the accuracy of
  frequency detection closely imitates the manner in which the
  human-auditory system responds to frequencies.

  Additionally, the close matching of notes in western scales renders
  the constant-Q particularly useful in note detection; identifying a
  musical note value rather than an explicit frequency value.
  Furthermore the constant Q simplifies the process of timbre analysis.
  The frequencies of a musical note played by an instrument are often
  comprised of harmonically related partials. The timbre of the
  instrument can be characterised by the ratios of the harmonics. With
  the constant Q transform, the harmonics are equally spaced across the
  bins regardless of the fundamental frequency. This greatly simplifies
  the process of identifying an instrument playing a note anywhere in
  the scale simply by shifting the characterisation across the bins. A
  potential downside to using the constant Q transform is that it
  demands more computation than the Fourier transform.

  An efficient algorithm for transforming a discrete Fourier transform
  (which may be computed with the FFT) into a Constant Q transform is
  detailed in Brown and Puckette (1992).


How can I enforce aliasing with Audacity? I would like to hear how aliasing sounds like. How can I get this with Audacity?

(I was told that I should adjust the Fourier coefficients in Audacity. I have no idea where I can do so / how I should adjust them.)

Here's a Python script that does what Matt L said:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Enforce aliasing by reducing an audio-file to every n'th frame."""

from scipy.io.wavfile import read, write

def main(inputFilename, outputFilename, n):
    rate, data = read(inputFilename)
    write(outputFilename, rate/n, data[::n])

if __name__ == "__main__":
    from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
    parser = ArgumentParser(description=__doc__,
                            formatter_class=ArgumentDefaultsHelpFormatter)
    parser.add_argument("-i", "--input", dest="inputfile",
                        help="input WAV file", metavar="FILE")
    parser.add_argument("-o", "--output", dest="outputfile",
                        help="output WAV file", metavar="FILE")
    parser.add_argument("-n",
                        dest="n", default=2, type=int,
                        help="every n'th frame gets played")
    args = parser.parse_args()
    main(args.inputfile, args.outputfile, args.n)

It might not be audible at factor = 2.  Try other powers of 2.


How to predict the spectrum of an aliased/folded signal without looking at the FFT? This seems like a trivial concept but I just can't get the hold of it. I have a signal with two complex sinusoids having frequencies -25MHz and 17MHz. Currently it is being sampled at 64MHz so I get two peaks of fft at -25MHz and 17 MHz. If I throw every other sample away the component at 17MHz shifts to -30MHz and the one at -25MHz to 14 MHz.

I know it is due to the aliasing/folding but why exactly are these frequency components present in the aliased signal? Is there any way to predict the frequency content of the aliased signal without looking at the fft? What spectrum should I get if every fourth sample is dropped (i.e. sampling frequency now becomes 16MHz)?

Your description of the situation is faulty.  A complex signal sampled at 64 MHz can represent frequencies from -32 MHz to +32 MHz, so the presence of tones at -25 MHz and +17MHz is fine.  However, when you decimate by two the sample rate drops to 32 MHz, and the frequency range becomes -16 MHz to +16 MHz.  In other words, you can't have a signal at -30 MHz afterwards.

The aliasing shifts the signals down to the new Nyquist zone (the region between the negative and positive Nyquist frequencies, -16 MHz and +16 MHz respectively).  -25 MHz is 9 MHz below the negative Nyquist frequency so it is shifted up by 32 MHz (the width of the Nyquist zone) to +7 MHz.  +17 MHz is 1 MHz above the positive Nyquist frequency, so it is shifted down by 32 MHz to -15 MHz.

The picture above shows the before (top plot) and after (bottom plot) frequency spectrum.  It was created with the following Matlab commands:

sig = exp(j*(1:1000)*2*pi*-25E6/64E6) + exp(j*(1:1000)*2*pi*17E6/64E6);
freq = -32E6:64E6/1000:32E6;
dec = sig(1:2:end);
freqDec = -16E6:64E6/1000:(16E6-1);
subplot(2,1,1)
plot(freq, 20*log10(abs(fftshift(fft(sig)))))
subplot(2,1,2)
plot(freqDec, 20*log10(abs(fftshift(fft(dec)))))

EDIT: I see that you largely had the right idea when you calculated the aliased frequencies.  The only mistake that you made was to not divide the frequencies by 2 due to the decimation.  Thus, -30 MHz should have been -15 MHz, and 14 MHz should have been 7 MHz.


Designing an IIR to generate a specific data sequence I know that it's possible to design an IIR with specific poles and zeroes to create specific frequency responses.

Is it also possible to design an IIR such that when you give it an impulse, it generates values of a desired sequence? Or at least a sequence within some tolerance of those values?

Yes, there are time-domain design methods for IIR filters. One of the best-known ones is Prony's method. It is well described in the book Digital Filter Design by T.W. Parks and C.S.Burrus (ch. 7.5).

If the desired sequence is (right-sided) periodic, then there's a very simple and exact solution possible (which is usually an exercise for DSP students).


sound source localization - finding source vector This may sound a stupid question - but I am working on some simple microphone array beamforming sound source localization stuff.  I think I understand how to use beamforming to find the TOA of the sound source between 2 microphones (positioned, say on the x axis), and therefore how to find the angle, phi. If the source is positioned normal to the z axis then this makes theta (the azimuth) equal to pi/2. However how do I find the r coordinate of the source vector?  As I understand it , I cant localize the source without this coordinate.

Any help much appreciated, thanks in advance

In general, 2 microphones are not enough.  In 3 dimensions, you actually need 4 (non co-linear) for the 3 polar coordinates.

There are cases where if you can bring more information to your problem, like knowing the signal of the source and the transmit and receive times of the signal, you can estimate the range for simple known geometries.

A good reference is:

Strang, Gilbert, and Kai Borre. Linear algebra, geodesy, and GPS. Siam, 1997


How to segment femoral and tibial knee cartilage individualy? I am researching in medical image processing area. I want to segment knee cartilage into femoral cartilage and tibial cartilage. But two part have same intensity and only different position(femoral cartilage is above and tibial cartilage is below). I am finding the techniques or method to segment them individualy.Could you help me show to me that techniques? Thank you so much![enter image description here][1]

Looks like you have to do some motion tracking of points to separate the two. Is it possible to take images before and after some motion? That way you could model how the cartilage moved as a function of tibia or femoral motion, and you might be able to separate what portions of the cartilage belong to the tibia and femur.


Power of a periodic sequence I am trying to find the way to reduce the standard expression to compute the power of a generic sequence $x(n)$:

$$P_{\text{x}}= \lim\limits_{N \to \infty}\frac{1}{2N + 1}\sum\limits_{n=-N}^{N}|x(n)|^2\text,$$

when $x(n) = x(n + N_0)$ ($N_0$-periodic sequence), into the simplified formula

$$P_{\text{x}}= \frac{1}{N_0}\sum\limits_{n=0}^{N_0-1}|x(n)|^2\text.$$

How can I demonstrate that second one could be obtained from the first one, using the hypothesis of periodicity? I've tried something like this, but I get stuck in after $\text{(iii)}$:

$$\begin{align}
N &= N_0\tag{i}\\
P_{\text{x}}&= \lim\limits_{N_0 \to \infty}\frac{1}{2N_0 + 1}\sum\limits_{n=-N_0}^{N_0}|x(n)|^2\\
&= \lim\limits_{N_0 \to \infty}\frac{1}{2N_0 + 1}\sum\limits_{n=-N_0}^{N_0}|x(n + N_0)|^2\tag{ii}\\
&m = n + N_0\text{, so }\\
P_{\text{x}}&= \lim\limits_{N_0 \to \infty}\frac{1}{2N_0 + 1}\sum_{m=0}^{2N_0}|x(m)|^2\tag{iii}
\end{align}
$$

The basic trick is to bound the series above and below. Let us do it on one side, for positive indices.

For any $N> 0$, you can write $N=kN_0+r_N$, with $0\le r_N< N_0$. Then if $a_n$ (here $a_n = |x_n|^2$) is positive, $\sum_{n=0}^{N-1} a_n$ is increasing. Now  $kN_0 \le N< (k+1)N_0$, hence you have:

$$ \sum_{n=0}^{kN_0-1} a_n  \le \sum_{n=0}^{N-1} a_n  \le   \sum_{n=0}^{(k+1)N_0-1} a_n   $$
or

$$ k\sum_{n=0}^{N_0-1} a_n  \le \sum_{n=0}^{N-1} a_n  \le   (k+1)\sum_{n=0}^{N_0-1} a_n   $$

since  $a_n$ is  $N_0$-periodic. Now since

$$\frac{1}{(k+1)N_0} \le \frac{1}{N} \le  \frac{1}{kN_0}$$
by a simple point-wise product of the threefold inequality (valid because positive):

$$ \frac{k}{(k+1)N_0} \sum_{n=0}^{N_0-1} a_n  \le \frac{1}{N}\sum_{n=0}^{N-1} a_n  \le    \frac{(k+1)}{kN_0}\sum_{n=0}^{N_0-1} a_n  \,. $$

Since $k/(k+1)\to 1$ and $(k+1)/k\to 1$ as $N\to\infty$, you have your result, mostly. You can split the original series in $\sum_{-N}^{N}a_n = \sum_{0}^{N} a_n+ \sum_{-N}^{0}a_n - a_0$,  take care of the positive and negative indices separately, and leave $a_0/(2N+1)$ tend to $0$.


Does a RHP zero imposes some limitation in the time domain? I am studying control systems and I have seen that a non-minimum phase zero imposes a limitation on the bandwidth, which is that the bandwidth has to be lower than the frequency which is half the frequency of the RHP pole.
I have also seen that it imposes in the step response an undershoot, but as far as I know there are no other influences of RHP zero in the time domain.
Does a RHP zero impose any other limitation in the time domain?

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

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

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


Get an input signal from LTI system output, find impulse response I just study signal processing and I have some questions.

How to find impulse response $h[n]$ sequence of the system if I have
$x[n]$ and $y[n]$ (input and output) sequences? Hope there is a formula
to do that.
How to find $x[n]$ input sequence if I have $h[n]$ and $y[n]$?
Is it possible to find $x[n]$ if I just have set of $y[n]$ values, but I
don't know impulse response $h[n]$? How to do that if so?

Maybe you can also suggest something for further reading about this topics.

Ok.  You are touching on the topic of deconvolution.  You can do this most easily by transforming the data into the frequency domain, doing some basic math and then converting back to the time domain.

For example for case 1, if you transform the signals you have
$ H(f) = Y(f)/X(f) $ so $h[n] = IDFT[Y(f)/X(f)]$

You can do similar manipulations to handle case 2.

I'm not sure case 3 is possible.  I think you need more information to produce x[n]

Here is a site that discusses deconvolution:  http://terpconnect.umd.edu/~toh/spectrum/Deconvolution.html

This is a VERY simplified treatment of this question.  I'm not sure what level of detail is being requested in the posted question.


Odd-order Butterworth filter gain at cutoff frequency Wikipedia gives the gain of an n-order Butterworth filter as

$$G^2(\omega)=\left |H(j\omega)\right|^2 = \frac {{G_0}^2}{1+\left(\frac{j\omega}{\omega_c}\right)^{2n}}$$

here. If you look at the gain at the cutoff frequency, the denominator becomes $(1+j^{2n})$, which blows up to infinity for odd $n$. Am I doing something wrong? Are Butterworth filters, or that formula, restricted to even $n$?

Congratulations, you've found an error on the wikipedia page on Butterworth filters. The squared magnitude of the frequency response of an $n^{th}$-order Butterworth low pass filter is

$$|H(j\omega)|^2=\frac{1}{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}\tag{1}$$

There should be no imaginary unit $j$ in that formula.


ROC of transfer function given:

$$ H(z) = \frac{4z(z-1)}{z-0.5} $$

I would say, when all poles are in the unit circles, the impulse response is right sided and causal.

For that transfer function, you have 2 regions of convergence (ROCs): $R_1 = \{|z|>0.5\}$ and $R_2=\{|z|<0.5\}$.

$R_1$ includes the unit circle, so the corresponding impulse response corresponds to a BIBO-stable filter. Since it is an outer disc, it is also right-sided.

Conversely, $R_2$ corresponds to a left-sided impulse response, and an unstable filter.

The extra $z$ in the numerator means that the impulse responses will be shifted 1 sample to the left, so the impulse responses are not causal.


Apply LTI filter to incoming signal stream Typically, when I want to filter a signal with an LTI filter (in this case a low-pass filter), I just convolve the signal with the filter.

However, what do I do if I have a stream of data that I want to filter as it arrives? Do I have to save a certain amount of the signal and filter it? Is there some other approach or theory that applies to this scenario?

As a quick draft of an answer, it turns out the things I was looking for was:

Digital Filter Realization
Digital Low Pass Filter Realization


About combination between QAM and PAM I have confused question regarding the possibility of combination between QAM modulation and PAM in order to improve the spectral efficiency.  Is that possible?

As you see in the below Figure, it's the QAM modulation, where the four status of 1+j, 1+j, -1+j and -1-j are existed. where J represents the imaginary part.

My question, what's about using the amplitude with that? For example, in case of combination with PAM modulation, we can have eight status instead of four, which are:
1+j, 1+j, -1+j and -1-j  at the same time the four added status will be 1/2(1+j), 1/2(1+j), 1/2(-1+j) and 1/2(-1-j), where the amplitude is used to represent additional bits to transmit too.

Is that reasonable and logical to be used ? Is there any articles or books which explain that ?

It is quite common to combine amplitude and phase modulation. QAM is actually a good example because generally you will have symbols with different phases and different amplitudes. Just in the special case of $4$-QAM (which is actually the same as QPSK), all symbols have the same amplitude. But if you take as an example $16$-QAM then you'll see that the information is encoded in the phase as well as in the amplitude:

[Chris Watts [CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0)], from Wikimedia Commons]

Note that in general the exact position of the complex symbols is a compromise between ease of implementation and power efficiency. QAM is especially easy to implement and the penalty compared to more power efficient constellations is relatively small.


Control theory: how do you initialize input for a model predictive controller? Although I think this is a pretty popular algorithm, I posted this question on electronics stacksexchange but received no response. I am posting it here to see if there are many signal/control experts who can answer my very simple question

On these slides, there is a description of MPC algorithm. I am confused as to how the input $u$ is initialized and how its values are known in the future.

How is $u$ initialized at the beginning?

In the  very beginning, we need all the future values of $u$ in order to predict the value of $x$. I am fully aware of how $x$ is predicted $N$ steps into the future.

But where do you obtain $u_k, u_{k+1}, \ldots, u_{k+N-1}$ to predict the states in the first place?
For instance how do you know $u_{k+1}$?

How to get the next $u$ values after receding horizon step?

After obtaining the optimal $u^*$ from the optimization program, we keep the first element of $u^*$ and discards the rest. This allows us to produce our output $y$.

But how are the subsequent $u$ generated in order to feed the next iteration of the program?
Again, we need every single value $u_k, u_{k+1}, \ldots, u_{k+N-1}$, but for instance how do you know $u_{k+1}$? (This is the same question as question 1 I suppose)

All in all, I am not sure how we know the future of our input $u$ in order to predict $x$. Can someone please enlighten me as to how this is done in simulation or in practice?

This algorithm in general tries to solve an optimization problem each step, defined as,

$$
\begin{aligned}
& \underset{\textbf{u}}{\text{minimize}}
& & \sum_{i\ =\ 0}^{N-1}\left[x^T\!(k+1+i)\, Q\, x(k+1+i) + u^T\!(k+i)\, R\, u(k+i)\right] \\
& \text{subject to}
& & x(k+1+i) = f(x(k+i),\ u(k+i)) \qquad \forall\ i = 0, \ldots, N \\
& & & g(x(k+i),\ u(k+i)) \leq 0 \qquad \qquad \qquad \ \ \ \forall\ i = 0, \ldots, N
\end{aligned} \tag{1}
$$

where,

$$
\textbf{u} = \begin{bmatrix}
u(k) & u(k+1) & \cdots & u(k+N-1)
\end{bmatrix}^T. \tag{2}
$$

This means that the algorithm tries to find the values for $\textbf{u}$ which minimize that quadratic summation expression, while still satisfying the update equality constraints and the inequality constraints (which returns a vector of length $m$). Here each $x(k)$ is the error between between the actual/predicted state and the desired reference state at $k$.

Depending on the functions $f(x,u)$ and $g(x,u)$ this problem can be hard to solve or even to prove that you found the global optimum and often require an initial guess. For example a generic solver for these kind of problem in Matlab can be solved with fmincon.

However if $f(x,u)$ and $g(x,u)$ are linear functions, then this problem becomes "a lot" easier to solve. A side note: if $g(x,u)$ is not linear, but is convex, then you can approximate it with linear inequalities. In Matlab the function quadprog can be used, since the const function is already in a quadratic form.

In chapter 2 of the paper it is shown how you can simplify the problem when this is the case. Namely all equality constraints and future states on the system can be removed. It also introduces $\bar Q$, which allows you to calculate the cost function from $N$ to infinity, assuming that you are no long constraint by the inequality constraints. Such an infinite horizon in the cost function helps ensuring stability of the system. Finding this $\bar Q$ involves solving the Riccati equation.

For example for a problem with no inequality constraints and only linear equality constraints, such that,

$$
x(k+1+i) = A\, x(k+i) + B\, u(k+i) \qquad \forall\ i = 0, \ldots, N , \tag{3}
$$

with $x(k) = x_0$, you can find all other $N$ states of the system as a function of $\textbf{u}$ with,

$$
\begin{bmatrix}
x(k+1) \\ x(k+2) \\ \vdots \\ x(k+N)
\end{bmatrix} = \underbrace{\begin{bmatrix}
A \\ A^2 \\ \vdots \\ A^N
\end{bmatrix}}_{\textbf{A}} x_0 + \underbrace{\begin{bmatrix}
B & 0 & \cdots & 0 \\
AB & B & \cdots & 0 \\
\vdots & \vdots & \ddots \\
A^{N-1}B & A^{N-2}B & \cdots & B
\end{bmatrix}}_{\textbf{B}} \textbf{u}. \tag{4}
$$

The cost function from equation $(1)$, denoted as $J(k)$ in your linked paper, can then also be written as,

$$
J(k) = \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right]^T\! \underbrace{\begin{bmatrix}
Q & 0 & \cdots & 0 \\
0 & Q & & \vdots \\
\vdots & & \ddots & 0 \\
0 & \cdots & 0 & Q
\end{bmatrix}}_{\textbf{Q}} \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right] + \textbf{u}^T\! \underbrace{\begin{bmatrix}
R & 0 & \cdots & 0 \\
0 & R & & \vdots \\
\vdots & & \ddots & 0 \\
0 & \cdots & 0 & R
\end{bmatrix}}_{\textbf{R}} \textbf{u}. \tag{5}
$$

After expanding the brackets, then equation $(5)$ can also be written as,

$$
J(k) = x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0 + x_0^T\textbf{A}^T\textbf{Q}\textbf{B}\textbf{u} + \textbf{u}^T\textbf{B}^T\textbf{Q}\textbf{A}x_0 + \textbf{u}^T \left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right] \textbf{u}. \tag{6}
$$

The first term of equation $(6)$ does not depend on $\textbf{u}$ and thus will be a constant. The second term is the transpose of the third term because $\textbf{Q}$ is symmetric, since their results are also a constant (a one by one matrix) they are equal and can be added together. Equation $(6)$ can therefore be simplified even further to,

$$
J(k) = \underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0}_c + 2\,\underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{B}}_{f}\,\textbf{u} + \textbf{u}^T\,\underbrace{\left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right]}_{H}\,\textbf{u}. \tag{7}
$$

Because it is assumed that $Q$ positive semi-definite and $R$ is positive definite, therefore the minimum of $J(k)$ can be found by taking its derivative with respect to $\textbf{u}$ and equating it to zero. This allows you to get the following solution for $\textbf{u}$, similar to equation $(2.6)$ of the paper,

$$
\textbf{u} = -H^{-1}f. \tag{8}
$$

So using the model from equation $(3)$ allows you to express the cost function only in terms of the future inputs. The solutions of problems with inequality constraints can be found in a similar way, but may require numerical algorithms however that is an entire field of study by itself.


MFCC - Significance of number of features I have been doing some readings on the computation of Mel-Frequency Cepstral Coefficients (MFCC) and further use of Vector Quantizers (VQ) for recognition purposes. I am however stumped by the method of computation for those MFCCs regarding the number of suggested frequency bands. For example, why is the recommended number of mel filterbank filters 20 -40? Why do we eventually end up with 13 coefficients and not 40 as will have been found from the 40 log energies we shall calculate?

Number of filter banks

One of the last steps in the MFCC's calculation is measuring the energy in the filter banks. We do that because want to reduce the dimensionality of our input vector (amplitude spectrum), as well as capture its envelope. Those triangular filters are spaced over the Mel scale:

This means that we have very good resolution in low frequencies. Exactly opposite is true for higher frequencies. We do that because MFCC's are suited for speech-related tasks and most of the information is located in lower frequencies (i.e. formants).

So how many filter banks do we actually want? For example HTK by default is using 20 filter banks. You might increase that number, especially if you are dealing with signals that contain a lot of closely spaced frequencies and you want to resolve them between each other. It's totally up to you - in the end what really matters is the classification performance.

Number of coefficients

After taking the logarithm of energies in each filter bank, the last step (or second last if you are doing some liftering) is to calculate the Mel Frequency Cepstral Coefficients. We do that by fitting the cosines to calculated log energies using the DCT. This captures the periodicity in the reduced spectrum. Figure below should help you in understanding that process. You can imagine that for number of coefficients equal to number of filter banks this would correspond to capturing the alternating energy between each of the filter banks.

So how many MFCC's do we want to calculate? HTK by default uses 12 and in most applications it is more than enough. In general we don't want too many coefficients because:

It's all about reducing the dimensionality of our feature space.
One DCT properties is that it de-correlates and keeps most of the information in first few coefficients.

Again it all depends on your application and you should adjust this number based on the recognition performance. From my experience the gain from increasing the number of MFCC's was negligible compared to introducing the $\Delta$ and $\Delta\Delta$ coefficients.


Phase of the FFT of a sine function I'm trying to wrap my mind about what happens in Matlab/Octave when you compute the fft of a sine function. Theoretically, the Fourier Transform of a sine is given by two delta functions multiplied by the imaginary unit (one placed at the frequency of the sine and another at minus the frequency). I can see that in the modulus of the FFT (not exactly a delta because of the frequency leakage due to the windowing), but when I check the phase I get a signal that is non-zero over all the support (from 0 to the sampling frequency). I'm playing with this toy example in Octave:

fs = 200,t=-10:(1/fs):10;
x = sin(2*pi*30*t);
f = linspace(0,1,length(x))*fs;
fourierTrans = fft(x);
figure,plot(f,abs(fourierTrans))
figure,plot(f,angle(fourierTrans))

And the picture I get for the phase is:

Why is it that I don't see two deltas in the picture for phase of amplitude pi/2 and -pi/2?

Cheers

The phase of a sinewave that is not exactly periodic relative to the FFT aperture length is discontinuous at the edges.  If you want a sane/visualizable phase (not the phase of a discontinuity) for non-periodic frequencies, try generating your sinewave with a known phase at the center of the FFT (sin(2*pi*foo*(t-N/2))), and use an fftshift to re-reference the reported phase result of the FFT to the center of the FFT aperture as well.  Then the oddness/evenness ratio of your test function, and thus the phase result will remain stable as frequency is swept between FFT bin center frequencies.

The phase of a non-periodic-in-aperture sinewave will contain the phase of the convolution of the signal with the transform of the FFT window (something like a Sinc, which ends up spanning the full width of the FFT), which will also have a phase related to that of the signal.  Even a computed FFT of an exactly periodic sinewave will also produce some values that aren't exactly zero in all the other bins due to quantization noise in the computer arithmetic; thus those essentially zero, but not quite zero, magnitudes will have a phase that may or may not be zero.  And the phase of an actual zero magnitude vector is indeterminate, not necessarily zero.


Moving Average of sinusoid If you calculate a simple moving average over a window length that is equal to the period of the sinusoid, you get a straight line (=0 because the wave is symmetric around the X-axis): the wave has a period of 40, equal to the length of the window of the simple moving average: so the filter weights are 1/40, 1/40,... 40 times. This is kind of obvious because a sinusoid is like a reflection of the same shape around its axis.

Which non-trivial moving average weights would you need to find to get a filter that returns the same output as its input (so besides SMA(1) )? So I want to find the filter weights of some type of average over some window length, that when multiplied with their respective values of the sinusoid, while sliding the window over the sinusoid, add to a sum equal to the sinusoid.
ps: is it possible to answer using as little signal processing / engineering terminology as possible please? Like you're talking to a mentally challenged person.. :)

The fundamental property of averaging filters (in swear words, a linear system) is that the output of a sine is a sine of the same frequency, albeit of zero amplitude. So:

Which non-trivial moving average weights would you need to find to get a filter that returns the same output as its input?

An averaging in phase with the sine: take a filter full of zeroes, except $K$ values, pairwise distant by a multiple of $40$, with coefficient values $1/K$. This way, you will average $K$ sine value of the same amplitude, yielding a value with the same amplitude.

Hence, you have a bunch of non trivial filters. However, they are not very interesting in practice, because they shall be tuned to the period of a sine sampled with particular sampling.


How to check if a signal is mainly composed of a few impulse peaks? I am working on a problem that produces different types signals, and I need to detect the peaks after low-pass filtering, then find_peaks from scipy. For example, the one below:

I am however faced with some signals that have an impulse shape, as in:

or,

Filtering such signals destroys the peaks' locations, and sometimes lead to merging two or more peaks into one. Thus, I am not able to extract the peaks correctly.
Is there a method to infer if the signal has impulses (a few to many spikes)? This would help me directly extract the peaks easily.
I highly appreciate any input.

You can use kurtosis as a measure of how 'peaky' your signal is. Or the flatness measure, which is the ratio of the geometric mean of the signal to its arithmetic mean. Any signal which is not relatively flat will have a flatness value near 0.
$$\text{Flatness} = \frac{\sqrt[N]{\prod_{n=0}^{N-1} x(n)}}{\frac{1}{N} \sum_{n=0}^{N-1} x(n)}$$


Output signal in digital communication In digital communications systems we say that the received signal can be expressed as

$$y(t)= x(t)+ n(t)\tag1$$

if x(t) is the transmitted signal and $n(t)$ is the noise.

I am wondering whether a better model would be

$$y(t)= x(t-t_0)+ n(t)\tag2$$

To account for delay caused by channel. Is my logic correct?

If so then, what are the assumptions made when using (1) and what are we accounting for if we use (2)? ie what type of errors are we accounting for (delay spread, or time acquistion error...)

I am looking forward to hearing your responses and discussion on this. thanks

Equation (1) in your question is a pure additive noise channel without any further distortion of the signal. Equation (2) takes a delay into account, but this is usually not relevant. A pure delay does not affect the quality of the communication link in terms of bit error rate. (I'm not talking about the annoying effect of delay in voice communication).

So the model in Equation (2) is often not relevant and can be replaced by Equation (1). A more general and useful model, which includes Eq. (2) as a special case, is

$$y(t)=(x*h)(t)+n(t)$$

where $h(t)$ is the impulse response of the channel, and $*$ denotes convolution. This model represents a linear time-invariant channel with additive noise. Eq. (2) is obtained as a special case with $h(t)=\delta(t-t_0)$.


Image blind separation I am working in image cryptography using an ICA algorithm like JADE. To evaluate my work, I am used the SDR and PSNR.

I would like to know what the acceptable values of the SDR and PSNR in image blind separation are.

The best way to compare your method to others is to run the algorithms on a common data set. Quntative measures such as PSNR vary significantly across images.

It is also important to consider subjective analysis. After all, your research is aiming for improvements in appearance.


Resample audio signal with a low factors Currently, we need to handle a piece of an audio signal. We want to resample it with different factor rate (p/q) starting with 1.0002 to 0.9998 with a decrement of 0.000005. But Matlab function cannot handle this because of the Integer Limit. Are there any ways to solve it?

You can resample to any, including floating point approximations of irrational rates, using reconstruction by Sinc kernel interpolation (realistically, by using a windowed Sinc interpolator).


What is conceptually happening in this modulation function? I know this is a programming question but i am trying to understand the logic of this. I have come across below function which is basically written in the groovy, which i think is BPSK modulation function.
private final int SAMPLES_PER_SYMBOL = 150
private final float NFREQ = 1/15

private float[] bytes2signal(byte[] buf) {
  float[] signal = new float[buf.length * 8 * SAMPLES_PER_SYMBOL * 2]   // 8 bits/byte, 2 floats/sample
  int p = 0
  for (int i = 0; i < buf.length; i++) {
    for (int j = 0; j < 8; j++) {         // Not sure what is the significant of this ?
      int bit = (buf[i] >> j) & 0x01      // Not sure what is the significant of this ?
      float f = bit == 1 ? -NFREQ : NFREQ
      for (int k = 0; k < SAMPLES_PER_SYMBOL; k++) {
        signal[p++] = (float)Math.cos(2 * Math.PI * f * k)
        signal[p++] = (float)Math.sin(2 * Math.PI * f * k)
      }
    }
  }
  return signal
}

My question is what is significant of this two lines?
 for (int j = 0; j < 8; j++) {         // Not sure what is the significant of this ?
    int bit = (buf[i] >> j) & 0x01     // Not sure what is the significant of this ?

Can someone please explain this to me ? What is logically happening with this buffer shifting by 8 bit and And-ing with 0x01. Below is the test function for iteration i=0
------------------------------------------------------------
i = 0
j = 0
1 >> 0 = 1
1 & 0x01 = 1
f = -0.06666666666666667
------------------------------
j = 1
1 >> 1 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 2
1 >> 2 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 3
1 >> 3 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 4
1 >> 4 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 5
1 >> 5 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 6
1 >> 6 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
j = 7
1 >> 7 = 0
0 & 0x01 = 0
f = 0.06666666666666667
------------------------------
------------------------------------------------------------

My question is what is significant of this two lines

The code receives a buffer of bytes but the modulation is based on a stream of bits. The outer loop loops over all bytes in the buffer. The
 for (int j = 0; j < 8; j++) {

This loops over all 8 bits in the current byte
int bit = (buf[i] >> j) & 0x01

This tests whether the j'th bit in the current byte is 1 on 0, where j = 0 produces the LSb and j = 7 produces the MSb.
If your current byte is 0b11001001, then the variable "bit" will be 1, 0, 0, 1, 0, 0, 1, 1 in each of the 8 loop iterations. Depending whether the current bit is 1 or 0, the modulated signal is generated in phase or out of phase.


PLL phase and frequency characteristics I'm currently working with BPSK carrier phase and frequency recovery in GNU Radio. I have been able to use Costas loop and Frequency locked loop (four carrier acquisition) which are based on the control_loop class (this class implements a PLL). Until now, I never needed to configure the control_loop class since the wrapper classes (Costas and FLL) do this by default.
As it can be seen below, the control loop (PLL) is based on the proportional-plus-integrator. The control loop constructor takes three variables
control_loop(loop_bw, max_freq, min_freq). max_ and min_freq variables are used to keep the frequency (in the integrator part of the loop) between max_freq and min_freq. I have looked into the Costas and FLL source code and saw how they are configured.

Let's start with Costas loop.

The loop is configured to have max_freq/min_freq of +1 and -1 respectively
The phase error input to the control loop i.e. (Φ_i[n] - Φ_o[n]) is bound between -1 and 1.
After the phase error has been filtered, the output phase i.e. Φ_o[n] is wrapped (between 2pi and -2pi). The frequency is forced into the max_freq and min_freq [1,-1] region.

With the FLL,

The loop is configured to have max_freq/min_freq of (2*PI/samples_per_symbol) and -(2*PI/samples_per_symbol) respectively.
Unlike the Costas loop, the phase error input to the control loop i.e. Φ_i[n] is not bounded.
Like the Costas, the output phase i.e. Φ_o[n] is wrapped (between 2pi and -2pi). The frequency is forced into the max_freq/min_freq of (2*PI/samples_per_symbol) and -(2*PI/samples_per_symbol) region.

I basically have three questions

Why is the input phase error of the Costas loop limited in the [1,-1] region?
How do you choose the frequency region [max_freq and min_freq] of the integrator part of the loop? why do we need to limit the frequency range?
Why is the output phase Φ_o[n] (in both Costas and FLL) wrapped between -2PI and 2PI?

I have added links to the source code of the Costas and the FLL blocks

Costas Loop: https://github.com/gnuradio/gnuradio/blob/master/gr-digital/lib/costas_loop_cc_impl.cc
FLL band-edge : https://github.com/gnuradio/gnuradio/blob/master/gr-digital/lib/fll_band_edge_cc_impl.cc

Regards,
M.

Before I address your questions, you should understand:

a. the integral branch of the loop filter maintains a average phase increment in units of radians/sample.  It is not a frequency value, though in the right context it can be converted to a frequency value, using your sample rate as a conversion factor.

b. the total output of the loop filter is an instantaneous phase increment which is used to advance the numerically controlled oscillator's phase.

So to your questions:

  Why is the input phase error of the Costas loop limited in the [1,-1] region?

It looks like an arbitrary limit to keep the instantaneous phase increment output of the loop filter to some sane number of radians/sample.  (Correcting more than $\pm\pi$ radians/sample is pointless.)

The average phase increment maintained by the integral branch is limited to $\pm 1.0$ radians/sample as well.  In a no-signal, noise processing situation, the integrator output can grow without bound, so the limiters are needed.  The $\pm 1.0$ radians/sample is just an arbitrary limit to keep the loop near sane operating values for when a signal does come in.  A Costas loop that is applying a correction of more than 1 radian/sample is probably not anywhere near converged anyway.

  How do you choose the frequency region [max_freq and min_freq] of the integrator part of the loop? why do we need to limit the frequency range?

In your application context, which you really haven't provided, you have to decide what are reasonable limits for the average phase increment coming out of the integrator.  You must provide limits, because the integrator is unstable and its output can grow without bound, which would make the PLL useless.

For a Costas loop, the spacing in radians of your constellation points and the number of samples/symbols will dictate what are sensible limits for the maximum & minimum increment of the NCO in radians for each sample processed.  Anything outside of the interval $(-\pi, \pi)$ radians/sample would be nonsensical due to phase wrapping.

For a PLL that is trying to lock on to a tone in a particular frequency range $[f_{low}, f_{high}]$ in Hz, you can use your sample rate, $F_s$, to convert those to limits in units of radians/sample, i.e. $\left[\dfrac{\pi f_{low}}{F_s/2},\dfrac{\pi f_{high}}{F_s/2}\right]$ radians/sample.

  Why is the output phase Φ_o[n] (in both Costas and FLL) wrapped between -2PI and 2PI?

Because computation of $e^{j\phi_o}$, used for setting the NCO output, is generally numerically better with smaller angle arguments in most libraries.  Phase is cyclic so it's OK to just wrap down to small magnitude angles.


Power spectral density of a PAM signal Say I have a PAM signal as $x(t) = \sum_n a_n h(t-nT)$ where $\{a_n\} \in \{\pm 1\}$ are equiprobable random binary bits and $h(t)$ is bandlimited to $[-1/(2T),1/(2T)]$. ($x(t)$ is random process and cyclostationary).

My first question is that can I compute the energy of $x(t)$ as:
$$E_x = E\left\{(x(t))^2\right\} = \int_{-\infty}^\infty \lvert X(f)\rvert^2 df$$
If yes, then my next question is where I get stuck in the solution:
\begin{align}
X(f) &= \int_{-\infty}^{\infty} \sum_n a_n h(t-nT) e^{-j2 \pi ft} dt\\
&= \sum_n a_n \int_{-\infty}^{\infty} h(t-nT)e^{-j2 \pi ft} dt\\
&= \sum_n a_n  H(f) e^{-j2 \pi fnT}\\
&= H(f)A(e^{-j2 \pi fT})
\end{align}
Therefore,
$$E_x = \int_{-\infty}^\infty \lvert H(f)\rvert^2\lvert A(e^{-j2 \pi fT})\rvert^2 df$$

I know the answer is $E_x = \frac{1}{T} E_h E_a$ but I cannot make the connection.

Since $x(t)$ is a random process, its Fourier transform generally doesn't exist. In its original form, the process $x(t)$ doesn't even have a power spectrum in the conventional sense because it is cyclostationary (and not stationary). It can be made stationary by adding a random phase $\theta$ that is uniformly distributed on $[0,T]$ and independent of $a_n$:

$$x'(t)=x(t+\theta)=\sum_na_nh(t+\theta-nT)\tag{1}$$

This process is stationary and its power spectrum can be computed, which is usually done by first determining its autocorrelation function, and then applying the Fourier transform to the autocorrelation function. The corresponding derivation can be found in most textbooks on digital communications:

$$S_{x'}(f)=\frac{1}{T}|H(f)|^2S_a(f)\tag{2}$$

where $H(f)$ is the Fourier transform of the pulse $h(t)$, and $S_a(f)$ is the power spectrum of the discrete-time random process $a_n$.

If we can assume that the symbols $a_n$ are uncorrelated, the power spectrum $S_a(f)$ is constant, $S_a(f)=P_a$ (where $P_a$ is the power of $a_n$), and $(2)$ becomes

$$S_{x'}(f)=\frac{1}{T}|H(f)|^2P_a\tag{3}$$

Finally, the power (not energy) of $x'(t)$ is the integral over its power spectrum (note that $h(t)$ is band-limited):

$$P_{x'}=\int_{-\infty}^{\infty}S_{x'}(f)df=\frac{1}{T}P_a\int_{-\infty}^{\infty}|H(f)|^2df=\frac{1}{T}P_aE_h\tag{4}$$

where $E_h$ is the energy of the impulse $h(t)$.


How to calculate the entropy of random normally distributed samples? I have a hardware device that outputs a normal distribution of readings. It's (fitted) parameters are mean = 35 units, std.dev. = 8 units and I read the values as 8 bit integers.

I'm trying to determine a theoretical rate of entropy generation per reading, so for example 3.4 bits /reading or the like. Is there some formula I can make use of?

PS.  As entropy is actually generated by the observer rather than by the process itself, I'm particularly interested in the significance of the bit depth of the individual readings. I.E. what happens to the entropy rate as I alter the sample depth? I feel that I need to elucidate this point with a (facile) example...

Imagine measuring men's foot sizes, and you're using whole inches.  They would therefore range from 1 inch to 12 inches. That's 12 values and thus contain 3.6 bits of entropy /sample if feet were uniformly distributed in size. Now if you  change to a metric rule ranging from 1mm to 305mm, you'll be acquiring entropy at the rate of 8.3 bits /sample. This demonstrates that whilst God determines the size of men's feet, the entropy generated from measuring them is actually determined by the type of rule (or analogue to digital converter sample depth).

This would all be very simple if my actual case was uniformly distributed, but unfortunately it's a normal distribution.

It's (fitted) parameters are mean = 35 units, std.dev. = 8 units and I read the values as 8 bit integers.

The choice of parameters hint at a normal distribution. Is this the case? You might want to obtain a histogram and perhaps substitute mean and standard deviation with modal mean and range. This will also help with an estimation of entropy.

  I'm trying to determine a theoretical rate of entropy generation per reading, so for example 3.4 bits /reading or the like.

A first step would be to estimate the Shannon Entropy. In this case, your alphabet is the range of 256 symbols that are possible to be formed with an 8-bit integer and you would be estimating the entropy of a source that generates these symbols with probabilities dictated by the histogram of a long enough sequence of symbols.

If these 256 symbols appeared with equal probabilities in your sequence then you would need all 8-bits to represent them. However, if a subset of symbols appeared with higher probabilities than others then, on average, it would be possible to represent this sequence with less bits per symbol. To determine Shannon's Entropy, obtain the histogram of a long enough sequence, divide each entry by the sum of all entries so that it is now a probability density function (PDF), i.e. it has a total sum of 1.0, and then calculate:

$$H(X) = -\sum_{k=0}^{N} x_i \log_2(x_i)$$

Where $X$ is your PDF (derived above) and $N$ is the length of the histogram.

  PS. As entropy is actually generated by the observer rather than by the process itself, I'm particularly interested in the significance of the bit depth of the individual readings.

Entropy is a characteristic of the source. This particular way of deriving entropy applies to a particular kind of source called ergodic and this is why I am also emphasising the "long enough sequence" (For more information please see this and specifically section 5.)

I am not sure what you mean by "significance of the bit depth". The idea is this: Suppose that you have an alphabet with 256 symbols. Obviously, you need 8-bits to represent them. But, if you were to discover that your source is actually generating numbers in the range $0..32$ for 98% of the time and a single 256 value for the rest 2% of the total length of a given sequence, then, for 98% of the time you would be wasting 3 bits, because $2^{5} = 32$ and you would need all 8 bits for only 2% of the time. Imagine these bits on a sequence:

00000011  00000100  00001110  00001001  00001100  10001101  00011001  00001110

If you could somehow squeeze those 3 bits out of your stream, you could save time of transmission (therefore fit more data through a channel) or storage space.

This is why figures such as 2.5 bits make sense. You can't define half a bit BUT over two symbols (i.e. 3 bits + 3 bits = 6 bits) you would be saving 1 bit (because 2.5 bits + 2.5 bits = 5bits).

For more information on encodings that allow you to effect exactly that "squeezing", please see variable length coding.

Hope this helps.

EDIT:
If you are absolutely sure that the source is normally distributed, then substitute $X$ for an expression of the normal PDF and the sum for an integral over the sufficiently probable figures


DFT conjugate of $X^*[k]$, how to prove its formula in terms of $x^*[n]$? Trying to prove that:

$$
X^{*}[k] = \sum^{N-1}_{n=0} x^{*}\left((N-n)\right)_N\ W_N^{nk}
$$

Where:
$$((x))_N \text{ = x  modulus N}$$

$$W_{N}^{nk} = e^{-j\ 2\pi / N}$$

So I start out with definition for DFT:

$$
X[k]=\sum^{N-1}_{n=0} x[n]\ W_{N}^{nk}
$$

Then I conjugate both sides of equation.

$$
X^{*}[k]=\left[\sum^{N-1}_{n=0} x[n]\ W_{N}^{nk}\right]^{*}
$$

$$
X^{*}[k]=\sum^{N-1}_{n=0} x^{*}[n]\ W_{N}^{-nk}
$$

then somehow the -n changes to N-n which i'm guessing is this:

$$((-n))_N = N - n$$

I don't really understand it... but, seems its equivalent if you do this inside of an exponent… to just randomly take a modulus of a negative signed number and replace it with its modulus equivalent?

$$
X^{*}[k]=\sum^{N-1}_{n=0} x^{*}[n]\ W_{N}^{(N-n)k}
$$

Then, I get stuck on the proof because somehow $x^{*}[n]$ changes to $x^*((N-n))$,

and $W_{N}^{(N-n)k}$ somehow changes to $W_{N}^{nk}$.

and thus, i'm unable to obtain final result of:

$$
X^{*}[k] = \sum^{N-1}_{n=0} x^{*}((N-n))_N\ W_N^{nk}
$$

still wondering what the rules are to apply to get there?

So you want to show $X^*[k]$ ?

Actually it follows very simply from the fundamental properties

$x[n] \longleftrightarrow X[k]$
$x[-n] \longleftrightarrow X[-k]$
$x[n]^* \longleftrightarrow X^*[-k]$

and by combining the last two you have:

$x[-n]^* \longleftrightarrow X^*[k]$

which shows that $X^*[k]$ is given by the forward DFT of the signal $x[-n]^*$.

Note that, due to periodicity of DFT sequences, the negative index can be repaced by $-n = N-n$, and hence $x[-n]^* = x[N-n]^*$ too.

Direct derivation is as follows:

$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{nk} $$

where $W_N= e^{-j \frac{2\pi}{N} }$.

Take conjugate of both sides:
$$ X^*[k] = (\sum_{n=0}^{N-1} x[n] W_N^{nk})^* $$

$$ X^*[k] = \sum_{n=0}^{N-1} x^*[n] W_N^{-nk} $$

replace $n$ with $-n$ so that forward DFT appaers:

$$ X^*[k] = \sum_{n=0}^{N-1} x^*[-n] W_N^{nk} $$

which again shows that $X^*[k]$ is given by DFT of $x[-n]^*$. Note again that periodicity of $x[n]$ is utilized to rearrange the limits.


Using 8 bit samples for speech comparison I am making a speech comparison app in javascript.

I found this dsp.SE post:

  Speech comparison algorithm for rating on similarities,

where I read it is best to use mfcc and dtw algorithms.

I found the Node.js implementation of mfcc which I'd like to use. However, it is only working with 8 bit audio samples. Will that affect the quality of speech comparison, and by how much? If the difference is big enough I will edit the algorithm to support 16 bit samples.

Will that affect the quality of speech comparison, and by how much?

That is impossible to tell without knowing what the Node.js thing does internally; I think it's a bit much too ask for us to search for what you meant. As a comment: signal processing in JavaScript sounds like a bad idea, performance-wise and development-wise; it's really not what JS was designed and optimized for, and there's a significant lack of libraries, let alone efficient ones. I always urge people to use the right tools for their job, and JS isn't that here, I think. </comment>

That being said:

From a pure signal point of view, you can model the effect of quantization as noise. The problem here being that this kind of noise is neither uncorrelated to the signal, nor necessarily white. To make things a bit harder  on the compensation side, its amplitude typically isn't Gaussian even. Oh, well, but here goes quantization noise power, a figure that is both very important to understand the maximum SNR you can get out of a digital system and doesn't say much as long as you don't say how well the signal processing can deal with this specific kind of noise.

Luckily, for "small enough" quantization steps, a little stochastic consideration¹ says that the assumption that quantization noise (QN henceforth) is additive is pretty justified.

Now, assuming your signal amplitude is really uniformly distributed, and your ADC being perfectly uniform, the amplitude $\text{SNR}_Q$ (signal-to-QN-ratio) for an M-bit ADC becomes

$$\begin{align*}
\text{SNR}_Q &= 2^M\\
\text{SNR}_Q \text{[dB]}&= 20 \log_{10}(2^M)\\
 &= 20  \log_{10}(2)\,M \\
&\approx 6 M\text,
\end{align*}
$$
which implies that for 8bit, your $\text{SNR}_{Q,8b} \approx 48\,\text{dB}$, and for 16bit $\text{SNR}_{Q,16b} \approx 96\,\text{dB}$.

Now, speech definitely isn't uniform in amplitude; it's a bit hard to justify this model without knowing what your recording looks like, but I'd rather say it's composed of sines; in that case, you get an additional $1.8\text{dB}$ noise for both cases.

Point here is that I doubt that a "real world, non-studio equipment, non-anechoic-chamber silence" speech recording will ever be any close to $48\,\text{dB}$, so probably, no, that's a fine choice, and by the way you ask this question, I kind of doubt it will be trivial to extend the algorithm (which probably uses a lot of elegant numerical math internally) to 16bit, anyway.

¹ the amount of noise power correlated to the signal necessarily being very limited, and an mutually independent set of sufficiently many i.i.d. realizations (Normal due to CLT) added up still being Normal if their moments weren't too different ...


How to interpret these filter coefficients? I am trying to make sense of the filter coefficients.

The author claims that they are running a 3rd order lowpass Butterworth filter and were willing to give their design. They are not the most communicative type and it is typically weeks until I am able to get any sort of answer from them.

The coefficients are as follows:

1.22293 1.394312 1.994313
-0.3339159 -0.6572268 -0.9943293

They also provided the transfer function:

$$
H_k(z) = \frac{b_{k0}+b_{k1}z^{-1}+b_{k2}z^{-2}}{1-a_{k1}z^{-1}-a_{k2}z^{-2}},\quad k = 1, 2, 3.
$$

Whey I plug in these coefficients into the FDA toolbox in MATLAB, I get the following frequency response (sampling rate was set to $20833\textrm{ Hz}$):

What is the correct way to interpret these numbers? I am somewhat out of my element when it comes to the Digital Filtering and could use some help.

Sorry, but that does not make any sense

A 3rd order BW lowpass can be described most simply with a single number: the cutoff off frequency
It's the not numerator a denominator since this would have 4 numbers and the numerator would have the shape of x*[1 3 3 1] and the denominator should start with 4.
It's not the poles and zeros since all the zeros are at -1 and one pole is conjugate complex
It sure as heck isn't state space either
The transfer function is a boiler plate text book copy and also mostly wrong for a third order BW (which has one biqaud and one single real pole section).

You are not missing anything obvious and having you guessing blindly doesn't serve any useful purpose. You have to get back to the client and tell them that you need proper specification or you can't get the work done.


Transition table / matrix for convolutional encoder K=3 understandings Im trying to understand convolutional encoder $K=3$ , and I almost understand it but still confused about what's called transition matrix which in other words it means

what input value ($0$ or $1$) would produce the next state, given the current state.

So if we have $K=3$ which means $3$ shift registers, there's $2^{K-1}$ possibilities which means there's $4$ states as Im showing down: (it's in other words state machine)

And I fully understand it, and Im stuck on how do I built the correspond transition matrix for this state machine?
I've seen on the internet something weird that the transition table for convolutional encoder $K=3$ is:
transitionMatrix = [0 1 1 55; 0 0 1 1; 55 0 55 1; 55 0 55 1]

it was written that value $55$ is value of STATE HISTORY .
I didn't understand this matrix and how this matrix describe the state machine that I attached the photo above for convolutional encoder $K=3$
Could anyone please help me out for this?
I already did the state machine and I understand it, I didn't understand how do I build the transition matrix. I'm posting here in order to understand it and not just copying it or reading it from a book without knowing how do I implement it.

The OP's question reveals many mis-understandings of the issues.
The convolutional code in question is a nonsystematic $[2,1]$ convolutional code with constraint length $K = 3$. The specification is incomplete but the complete specification (we also need to know what the generator polynomials are) can be deduced from the state transition diagram that the OP has provided. This particular convolutional code is used as the first example in almost all texts on coding, and so the code is sometimes referred to as the ad nauseam code. Let's pick apart what things mean.  There is 1 input stream of bits entering the encoder during each clock cycle and 2 streams of bits leaving the encoder (in parallel, on two separate wires) during the clock cycle. The encoder itself contains one shift register (not $3$ as the OP thinks) of length K-1 = 2 meaning that it is capable of storing two bits. The encoder also contains XOR gates for computing the output bits from the encoder input bit and the shift register contents. The bits stored in the shift register are the immediate past input bit and the input bit immediately previous to the immediate past input bit, that is, the input bit from two clock cycles ago. Thus, this single two-bit shift register can be in one of $2^{K-1} = 4$ states rather than three shift registers somehow managing to be one of only $4$ states which has puzzled the OP (and would puzzle anyone else too). At the end of the clock cycle, the current bit enters into the shift register, becoming the "immediately past input bit" for the next clock cycle while the bit that served as the  "immediately past input bit" for the current clock cycle moves over one place to become the "input bit two cycles ago" for the next clock cycle.  At the start, the shift register is initialized with two $0$s in it. That is, the encoder is in state $00$ at the beginning of the encoding process.
During any clock cycle, K = 3 bits  viz. the current input bit, the immediately past input bit and the input bit from two clock cycles ago, are available for computing the 2 output bits,.  The state transition diagram of the encoder is a directed graph whose vertices labeled with the state that the encoder is in, and the label of the directed edge from one state to another state denotes "current input/output bits". For example, if the current state is $00$ and the input is a $0$, the encoder outputs $00$ and remains in state $00$. The edge showing this transition is thus labeled $0/00$ and the edge leaves and enters state $00$. On the other hand, if the input is a $1$, the output is $11$ and the encoder moves from state $00$ to state $01$. This the directed edge from $00$ to $01$ is labelled $1/11$. We also deduce from this that whomsoever drew the logic map for the encoder had the input entering the encoder from the right and the shift register shifting leftwards. If the input bits had been entering from the left and the shift register shifting rightwards, the next state would have been labelled $10$; you pays your money and you makes your choice.
Note that there is no edge from state $00$ to state $11$ or vice versa; there is no way that this particular encoder shift register can make a transition from state $00$ to state $11$ or vice versa in one clock cycle. Note also the symmetries in the state transition diagram. The two edges leaving each state have complementary input labels (obviously, if you think about it a bit, pun intended) and complementary output labels too. The two edges entering each state also have complementary output labels (but the same input label; can you see why?).
So, how does the encoder actually work? Well, if the input data stream is
$$ \cdots d[n-4]\quad d[n-3] \quad d[n-2] \quad d[n-1] \quad d[n] \quad d[n+1] \quad \cdots$$
then at the $n$-th clock cycle, the encoder shift register contains $d[n-2]$ and $d[n-1]$ as indicated below via the box, and $d[n]$ is the current input to the encoder shift register.
$$ \cdots d[n-4]\quad d[n-3] \quad\bbox[yellow,5px,border:2px solid red] {d[n-2] \quad d[n-1]} \quad d[n] \quad d[n+1] \quad \cdots$$
During this clock cycle, the encoder produces two code bits
\begin{align}
c_i[n] &= d[n]\oplus d[n-1] \oplus d[n-2],\\
c_q[n] &= d[n]  \oplus d[n-2].
\end{align}
The shift register then shifts, and during the $(n+1)$-th clock cycle the register contents are as shown below and the input data bits is $d[n+1]$.
$$ \cdots d[n-4]\quad d[n-3] d[n-2] \quad \bbox[yellow,5px,border:2px solid red] {d[n-1] \quad d[n]} \quad d[n+1] \quad \cdots$$
The generator matrix for this code is
$${\bf G} = \big[ 1 + x + x^2 \quad 1 + x^2\big]$$
meaning that if we regard the input data stream as a series
$$d(x) = d_0 + d_1x + d_2x^2 + \cdots,$$ then the encoder output streams can  also be regarded as series
\begin{align} c_I(x) &= (1 + x + x^2)d(x)\\
&= d_0 + (d_1 + d_0)x + (d_2 + d_1 + d_0)x^2 + \cdots\\
& \qquad \quad \cdots + (d_n + d_{n-1} + d_{n-2})x^n + \cdots\\
c_Q(x) &= (1 + x^2)d(x)\\
&= d_0 + (d_1)x + (d_2 + d_0)x^2 + \cdots\\
& \qquad \quad \cdots + (d_n + d_{n-2})x^n + \cdots
\end{align}
or, in short,
$$d(x){\bf G} = \big[c_I(x) \quad c_Q(x)\big] = \big[(1+x+x^2)d(x) \quad (1+x^2)d(x)\big].$$
Dedicated DSPissers will replace $x$ by $z^{-1}$ and say "Hey, you are telling me that $c_I(z) = (1+z^{-1} + z^{-2})d(z)$ and $c_Q(z) = (1+z^{-2})d(z)$, and I know that multiplication in the $z$ domain corresponds to convolution in the time domain, and so why don't you name these codes as convolutional codes and save a lot of explanations on dsp.SE?" and that's exactly what it is: the encoder output bits are obtained by convolving the encoder input sequence with the generator polynomial sequences.

I've seen on the internet something weird that the transition table for convolutional encoder $K=3$
is: transitionMatrix = [0 1 1 55; 0 0 1 1; 55 0 55 1; 55 0 55 1]; it was written that value $55$ is value of STATE HISTORY .
I didn't understand this matrix and how this matrix describe the state machine that I attached the photo above for convolutional encoder $K=3$

I too have no idea what that gobbledygook means.
In a comment, the OP asks about implementing a transition matrix. There may be a language barrier here in communicating with each other. I don't know what exactly implementing a transition matrix means or why one would want to go about implementing one. I know about people implementing the encoder from the state transition diagram, whether in hardware with flip-flops and gates, or via FPGAs, or in MATLAB or Python or C++ etc. When the generator matrix is not known, it takes a little analysis of the state transition diagram to figure out what the generator polynomials are, but once that is done, implementing the encoder is quite straightforward. Another way of looking at it is that the encoder implementation implicitly implements the transition matrix in that the transitions that the encoder makes follows the rules laid out by the transition matrix.


Smoothed image produces contour-like ripples, yet it disappears in saved image? I smoothed a grayscale image using a size 10 gaussian filter and display it in Octave (imshow) which resulted in ripples which appears to be contours for image intensity.

Yet if I save the image (imwrite), I get a proper smoothed image.(The ripple image above was taken as a screenshot.)

My guess is that imshow, for some reason, has a lower color depth at which it displays greyscale images than your image viewer.

The contours are actually just the quantization steps of the lower number of shades displayable.


Using under-sampled data, what determinations can be made about the Nyquist rate (correct sampling rate)? Using under-sampled data, what determinations can be made about the Nyquist rate (correct sampling rate)? If my instrumentation limits me to sampling data at a fixed but insufficient rate what, if anything, can be determined from that data about the sampling rate that would be required to properly represent the signal? Are there signal analysis techniques that would provide any indicators?

This may be a trivial answer, but you can of course use the aliasing information to determine the correct sampling rate, given that you know that the signal is undersampled. In the simplest case, given that the signal appears to be at $k f_s,\ k < 0.5$ one can devise that the frequency of the signal is $(1-k)f_s$. In general it may be at $(N \pm k)f_s$ for some integer $N$.


filter design for ECG I would like to ask how to design a digital filter to remove $50$ Hz interference on an ECG signal. Do I need to use notch filter?
Thanks.

Yes, you will need to make a 50Hz notch filter.

I would suggest you to look at this post.

If you tell what type of program you use (Python, R, etc) I could try to provide you some code aswell.


Determining the Nyquist rate of signals added and multiplied together Suppose I have a signal in the frequency domain $x(\omega)$.

Suppose I know that its Nyquist rate is $\omega_0$ (and let $\omega_m$ be the its maximum frequency).

I need to get the Nyquist rate of

$$x(\omega)+x(\omega)e^{-j\omega}$$

The answers say it's still $\omega_0$.

Why is it not

$$\omega_m+\omega_me^{-j\omega}$$

and if we substitute $2\omega_m = \omega_0$

$$\frac{1}{2}\omega_0+\frac{1}{2}\omega_0e^{-j\omega}$$

Let $Y(\omega)=X(\omega)+X(\omega)e^{-j\omega}$ with $X(\omega)$satisfying $X(\omega)=0$ for $|\omega|>\omega_m$. Then, by its definition, $Y(\omega)$ also satisfies $Y(\omega)=0$ for $|\omega|>\omega_m$, i.e., the maximum frequency remains unchanged, and, consequently, also the Nyquist frequency remains unchanged.

Note that if $x(t)$ is the signal corresponding to the spectrum $X(\omega)$, then $X(\omega)e^{-j\omega}$ corresponds to the signal $x(t-1)$, which clearly has the same maximum frequency as $x(t)$.


16-APSK NDA Carrier Synchronization Problem I'm working on a DVB-S2 Symbol Timing + Carrier Synchronizer with a focus on 16-APSK | 32-APSK modulation Scheme.
I understand that the first section could be achieved using a Polyphase Clock Sync approach (i.e. as implemented in GNURadio), the second functions is usually done using a Costas loops
(i.e. in the simple case of M-PSK modulation scheme)
Moreover, The Costas Loop used for carrier synchronization are designed to only accommodate PSK modulations. As such, a Costas loop used on any of the APSK modulations, may try to lock on to any of the PSK modulations corresponding to any of the constellation rings, each of which is a PSK constellation in itself. For example, for 32-APSK, a Costas loop may try to lock on to the 4 points in the inner ring, the 12 points in the middle ring, or the 16 points in the outer ring. The result will be an overall failure to achieve any reasonable carrier synchronization lock (figure represent the output of a 4-order Costas Loop after the Polyphase Clock Sync Scheme)

I'm interested in understanding the NDA or NDA-DD signal processing scheme behind and if there's some implementation for GNURadio environment.

Symbol timing and carrier synchronization are two separate "tasks" in a demodulator.
There is lot of literrature on the internet to understand synchronization on a DVB-S2 signal.
But not sure that GNU Radio has the good ready block for it.
For symbol timing synchronization, there is Gardner timing error detector that suits well for QPSK modulation. An adapted version of it for 16APSK and 32APSK in a second order control loop could do the job.
For carrier synchronization, if you try to look at the 16APSK and 32APSK in polar representation (the constellation), the angle value modulo a constant gives you a phase error. Also in a second order control loop you can track the carrier this way. The modulo constant depends on the order of the modulation scheme you use. That is probably also why the one in GNU Radio may not work.
What I described before is called fine synchronization. Which means that it can be not enough to lock on large offsets.
You also need a coarse synchronization to help the fine one to lock. Depends on your requirements.
I hope this helps


Compressed Sensing vs Common Compression Approaches I have a very basic question regarding common compressing approaches. Often, I have heard something like "The problem in our common way to compress data is, that we have to acquire a huge number of samples from our signal, because of the Nyquist-Shannon theorem. But then, when we compress these samples, we throw most of the information away."

I do not really understand the last sentence, probably because of lack of knowledge. So, my question is: what do you throw away and why?

Here are some sources:

Michael Lustig, at 0:26 he says

  Why we need to acquire all the data in the first place, and then throw most of it away.

Mark Davenport, on Slide 4 at 0:45 he says

  The very first thing that we do after we acquire some data is we compress it and seemingly throw away a lot of information that we painstakingly acquired.

Compression works by reducing redundancy and there are a number of ways that this can be done.

Very briefly, consider the pair of functions $C = f(S_{original}, \Theta_{f}), S_{recovered} = g(C, \Theta_g)$, where $f,g$ are the compressor and decompressor functions respectively, $\Theta$ is the set of parameters specific to $f,g$ and $C, S_{original}, S_{recovered}$ are the compressed, original and decompressed sequences respectively. If the compression is loseless then the error between the original and recovered sequences should be zero ($S_{recovered}-S_{original}=0$).

However, we might accept some small error so that $S_{recovered}-S_{original}<=e$. In this case, the compression is lossy.

What $f$ does is to reduce the redundancy in $S_{original}$ so that it can be represented using less data so that the length of $C$ is less than the length of $S_{original}$ ($|C| < |S_{original}|$).

Run-length encoding achieves this by looking for sequences of repeated values and substituing them for a pair of "value, count". The redundancy here is most obvious.

Huffman encoding, does this by creating a carefully constructed "dictionary" of codewords based on knowledge about the source of a signal. In this, frequently occuring values are mapped to short codewords and infrequently occuring values are mapped to longer codewords. In this way, instead of representing all values with the same amount of data (for example 8 bits), frequent values are represented with less data (less than 8 bits) and on average lead to describing a long $S_{original}$ with a short $C$. More frequently appearing characters than others, again, imply redundancy.

Lempel-Ziv-Welch encoding reduces redundancy by looking for repeating sequences. Once a repeating sequence is found it is placed in memory and assigned a codeword. The next time the same sequence is found it is substituted with its codeword. In this way, long sequences (of more than one values) are substituted with just one value (an index into the dictionary).

Other methods rely on preserving just a part of the frequency spectrum of a signal or image therefore reducing the overall ammount of data that needs to be described.

The key thing to remember here is redundancy. Repeating values, slow varying signals (which again implies nearby samples having similar values) or even components with negligible contribution. If a signal contains redundancy then it is predictable and if it is predictable then it can be captured in a model. (As a side note, please note that this redundancy, or amount of information in a signal can be quantified by Entropy).

So, compressed sensing is in a way similar to lossy compression in the sense that we are trying to recover an $S_{original}$ from a (shorter) $C$ accepting some small error $e$ in the recovery with some "decompressor" $g$.

BUT, to do this we need to make certain assumptions about the dynamics of $S_{original}$. For example, it is varying slowly, therefore it is not necessary to sample the full spectrum of frequencies that could generate some potential $S_{original}$

Another way to think about this is to look at a picture or a drawing and cut (or black)out a rectangle from it. That is, throw away some data or, don't acquire a part of data from the beginning. If the picture is predictable enough then it is easy to reconstruct the missing part. WHY is it predictable? Because it is a repeating pattern and we assume that the missing bit behaves similarly. If the picture is less predictable then we can still attempt to reconstruct a rectangle of it from the surrounding information but perhaps it will not be as accurate.

If we were to take this to the extreme (as it is done in practice anyway), then we could say that we don't really need to acquire all of the pixels of the image. We can acquire some of the data and by exploiting the redundancy of the signal and the fact that it is not going to be too unpredictable (i.e. we can model it) reconstruct the missing bits...almost as closely to reality.

This is exactly what iterative reconstruction techniques (a family of $g$ "decompressors") are doing by progressively trying to minimise the error between an $S_{recovered}$ at some iteration $n$ and the subset of known values contained in the actually acquired $C$ and stopping this process once the error drops below an "acceptable" $e$.

For more information please see this link which is a relatively simple iterative reconstruction scheme that shows how to "update" the unknown values from the known values. For even more information please this link.

Hope this helps.


understanding a signal from its mathematical formula I can understand from unit step function in the end of the equation that $n=[3 \quad4]$, but I cannot understand how delta function affects the signal in order to solve it and put it on matlab.

Answer:

What is misleading is that you ought to find the values of $x_1[n]$ for all $n$.

For the first term, assuming convention $\cdot^0=1$, you get $\delta[n]+\delta[n-2]$. For the second one, all terms except for $k=2$ cancel each other, so you only gain another $\delta[n-2]$. Then, as you guessed, the last term (without the initial minus sign), is 1 for $n=3,4$, and zero elsewhere.


How Does the Non Local Means Filter Work? As I understand,

A local mean filter is the one in which we take a pixel and calculate the mean of color-values of the pixels around them (but in a certain range of area) and the center pixel is replaced by that mean value.
The filter always has a specific size. For instance, 3x3, 5x5, and so on.
The filter sweeps the whole image row by row to achieve a filtered image (spatial domain), or, the FT of the image is multiplied by the FT of the filter (frequency domain).

But, I could not understand the case of a Non-Local Mean (or, median, or whatever) filter.

How is the the value of the non-local mean calculated?
How does that value replaces a pixel?
How does that filter sweeps the image?

Check out the description here.

Basically, it tries to find patterns of pixels all over the image, and average over these self-similarities rather than just using the pixels close to the current pixel.

So to answer your questions:

How is the the value of the non-local mean calculated?

Find all the pixels that are similar to the current pixel (non-local).

Find the average value of these.

How does that value replaces a pixel?

Find the non-local average as above.

Replace the current pixel with this non-local average.

How does that filter sweeps the image?

Well, provided you are creating a new image from the old one so that the mean calculation is being done on the old image and the mean assignment is being done on the new image, it shouldn't really matter.


Sample Rate & Highest Frequency Would I be right in saying that if a signal was sampled every 0.2ms, when converting it to digital. The sampling rate would be 5000(Hz)?

How would I go about working out the highest frequency it would have using the sampling theorem?

If you sample every 0.2ms you collect 5000 samples per second and hence the sampling rate is 5000 Hz. According to Nyquist, this means the highest frequency the signal can have so we don't experience aliasing would be 2500 Hz.

Though that wasn't the question (it was about the sampling theorem explicitly), it's worth pointing out that that's only a theoretical limit. In practice, we need to go a bit higher (leave some room) for various reasons, including the fact that we can never observe an infinitely long signal, as correctly pointed out by the commenters.


Digital beamforming and phase shift stage insertion I have read some literatures about the different topologies of beamforming and it came to my mind a question : where is applied the phase shift for each channel depending on the type of topology/technology used?

For analog beamforming, it is obvious as there is only one signal, the modulated carrier signal which is splitted to the different channels and phase shifting occurs afterward. The phase shifting stage comprises a gain/attenuation stage and a phase shifter stage. If the modulated signal is built by a digital stage, there is one DAC used for all channels prior to mixing and phase shifting.

In digital beamforming, the topology is a bit different: one DAC is used for each channel. This enables more possibilities, among them the possibility to phase shift digitally each channel when defining the modulation. However two types of digital beamforming exist:

one with LO mixing if the DAC is not fast enough to produce a RF signal,
another using RF-DAC where no LO mixer is required.

As in the first situation LO signal is generally the same for all channels, no phase shift is applied on LO-IF signal. Is the phase shift on IF modulation signal is enough or is it needed to phase shift the signal after mixing (to obtain an equivalent behavior as in analog beamforming?).

Best regards

In a first analysis, it does not matter where the phase shifting is inserted, as long as the mixing stage is linear.

Indeed, shifting a temporal signal is equivalent to multiplying it with a complex number integration the phase shift $\theta$ :

Unshifted signal : $s(t)$
$\theta$ phase shifted signal : $s'(t) = s(t).e^{j\theta}$
$f_0$ carrier mixed unshifted signal : $s_c(t) = s(t).c(t) = s(t).e^{2j\pi f_0.t}$
$f_0$ carrier mixed $\theta$ shifted signal : $s'_c(t) =s(t).e^{2j\pi f_0.t}.e^{j\theta}$ = $s(t).c'(t)$

where $c'(t) = e^{j(2\pi f_0.t+\theta)}$ is the shifted carrier.

In the end phase shifting the signal to be modulated is equivalent to phase shifting after carrier modulation.


How to handle edge effects with EMD (HHT)? In the application of the Huang Hilbert Transform (HHT) (or aka EMD, Empirical Mode Decomposition), there is the notable challenge of handling the edges (boundaries) of the signals; as noted in the answer link to its mention in other question.

What are some practical techniques that can be applied to signals in which there is 'noise' padding already present and when there is no redundant signal information at the boundaries? (eg how would a mirroring strategy be applied?)

One possibility is to mirror the extrema near both edges with respect to the edge time values. I think this is one of the best ways to reduce edge effects. This paper explains better the way to do it.
https://www.sciencedirect.com/science/article/abs/pii/S0888327007002701
Hope this is helpful to you.


How to differentiate a time domain signal in the complex transfer function? I have an input-output data set where the input is current and the output velocity. I am interested in the transfer function from current to acceleration though. So suppose:
$H(s) = \frac{I(s)}{V(s)}$
where $V$ is velocity and $I$ is current.
I thought that the solution would be to
$H'(S) = \frac{I(s)}{sV(s)} = \frac{I(s)}{A(s)}$ where A is the laplace transform of acceleration resulting in the transfer function I'm after. Now to MATLAB: I use 'tfestimate' to obtain the first transfer function $Txy$. Now I would suppose that substituting $s$ by $j \omega$ and implementing that to obtain the magnitude as $|\frac{Txy}{j\omega}|$ would yield the correct results. But I only get the correct results by using $|Txy \cdot j\omega|$ and comparing this to the differentiated signal of $v$. What error am I making? Is it actually the complete transfer function I want to multiply with $j\omega$ to achieve the differentiation instead of the denominator? Any help would be much appreciated!

The transfer function is commonly defined as output divided by input, and not so often the other way around (as you did in your question). Also Matlab uses that more common definition. So the transfer function $H(s)$ between current $I(s)$ (input) and velocity $V(s)$ (output) would be defined by

$$H(s)=\frac{V(s)}{I(s)}\tag{1}$$

Since acceleration $A(s)$ is the derivative of velocity, we have

$$A(s)=sV(s)\tag{2}$$

Consequently, the transfer function between current and acceleration is

$$G(s)=\frac{A(s)}{I(s)}=\frac{sV(s)}{I(s)}=sH(s)\tag{3}$$


How does multi-rate processing work? In the case with a 1-D input like audio, how can a system with one input do many different kinds of processing at different sampling rates and still, at the end of the processing, produce one waveform? For example an ADC/DAC reading and writing at 48kHz, but has an FFT running at 8kHz or another set of processes running at 16 kHz? Is it a matter of decimating and interpolating? How would all of these processes line up?

You are on the right track.  You can use decimation and interpolation (or resampling, which is a combination of both decimation and interpolation) to change the sample rate.  For instance, if you start with 48 kHz samples you can decimate by 6 to get down to 8 kHz, do some processing, then interpolate by 2 to get to 16 kHz, do some processing, then interpolate by 3 to get back to 48 kHz.

"Lining up the processes" is not in and of itself difficult.  You can decimate and interpolate whenever you want.  Doing it in a way that makes sense and is useful, on the other hand, requires knowing more about what you start out with and what you want to accomplish.