Numpy fft2 frequency shift I don't understand how to make frequency shift in fft2 or higher dimensions.  Could anyone explain it, please?

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

%matplotlib inline

numb_steps = 100
x = np.linspace(-10, 10, numb_steps)
y = np.linspace(-10, 10, numb_steps)

X,Y = np.meshgrid(x,y)

R = 1

Z = np.zeros((100, 100))
# circle
for i in range(numb_steps):
    for j in range(numb_steps):
        if x[i]**2+y[j]**2 > R:
            Z[i,j] = 0
        else:
            Z[i,j] = 1

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y,Z)

I only can get those plots

fig = plt.figure()
ax = fig.gca(projection='3d')

Z_fft = np.fft.fft2(Z)
FreqCompRows = np.fft.fftfreq(Z.shape[0],d=2)
FreqCompCols = np.fft.fftfreq(Z.shape[1],d=2)
FreqCompRows = np.fft.fftshift(FreqCompRows)
FreqCompCols = np.fft.fftshift(FreqCompCols)

S,D = np.meshgrid(FreqCompRows, FreqCompCols)

ax.plot_surface(S, D, np.abs(Z_fft))

plt.imshow(np.abs(Z_fft))

plt.xticks(range(len(x_ax)), np.round(x_ax))
plt.yticks(range(len(y_ax)), np.round(y_ax))

P.S. Yes, the ticks overlap is my problem as well

You can solve this using scipy.fftpack (sfft) instead of np.fft, because the sfft implementation can be directly used on 2-dimensionnal arrays so you don't have to do it in a convoluted way (ba dum tsss).

In your code header, add :

import scipy.fftpack as sfft

Then, to calculate the fft and shift the spectrum, use :

Z_fft = sfft.fft2(Z)
Z_shift = sfft.fftshift(Z_fft)

The obtained spectrum is then nicely arranged for image display :

plt.figure(4)
plt.imshow(np.abs(Z_shift))

Also, the way you are constructing the circle seems overly complicated, you can take advantage of python's syntax using boolean syntax :

X,Y = np.meshgrid(x,y)
Z = X**2+Y**2
R = 1
Z[Z > R] = 0;Z[Z > 0] = 1;

What I did here is using the meshgrid you created to calculate the circle equation result, in the form of Z. Then I test Z for all values above R, put them to zero, and the other values to 1.

My code is below :

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.fftpack as sfft # you have to import this package

plt.close('all')

numb_steps = 100
x = np.linspace(-10, 10, numb_steps)
y = np.linspace(-10, 10, numb_steps)

#constructing the disk in a simpler manner
X,Y = np.meshgrid(x,y)
Z = X**2+Y**2
R = 1
Z[Z > R] = 0;Z[Z > 0] = 1;

# Here, use sfft routines instead of np.fft
Z_fft = sfft.fft2(Z)
Z_shift = sfft.fftshift(Z_fft)
FreqCompRows = np.fft.fftfreq(Z.shape[0],d=2)
FreqCompCols = np.fft.fftfreq(Z.shape[1],d=2)
FreqCompRows = np.fft.fftshift(FreqCompRows)
FreqCompCols = np.fft.fftshift(FreqCompCols)

S,D = np.meshgrid(FreqCompRows, FreqCompCols)

# plotting your disk signal image
fig = plt.figure(1)
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y,Z)

# plotting "3D" fft non-shifted
fig = plt.figure(2)
ax = fig.gca(projection='3d')
ax.plot_surface(S, D, np.abs(Z_fft))

# plotting images of FFT's
plt.figure(3)
plt.imshow(np.abs(Z_fft))

plt.figure(4)
plt.imshow(np.abs(Z_shift))


What kind of filter is used here? I'm very new to the DSP world, and honestly I'm not even sure if this is the right place for this question.
Anyway, I found this piece of code on github, it's a kind of filter (I think), and I would like to know if this can be identified as a commonly used or well known concept. Does it have a name I can use to google? Where do I start to learn more about this, like how does one come up with this code.
The (relevant) code:
This is calculated once:
// Set the loop bandwidth to bw Hz.
w = 6.28f * 20 * bw * bsize / fsamp;
w0 = 1.0 - exp (-w);
w = 6.28f * bw * ratio / fsamp;
w1 = w * 1.6;
w2 = w * bsize / 1.6;

z1 = 0.0
z2 = 0.0
z3 = 0.0

This is calculated every loop:
// Run loop filter and set resample ratio.
z1 += _w0 * (w1 * err - z1);
z2 += _w0 * (z1 - z2);
z3 += _w2 * z2;

rcorr = 1 - z2 - z3;

In these snippets err is the input of the filter and rcorr is the output. bsize, ratio and fsamp are some constants.
For completeness, this code is from here on github.

It's part of a feedback loop -- presumably setting the sample rate in some real-time audio thingie, since I see keywords in the code like 'resample' and 'jack' and 'PLAY'.
// Run loop filter and set resample ratio.
So, this tells you what it does!
Elsewhere in the code, err is computed as the timing error between the input stream and the output stream.  So, this is a feedback loop.
_z1 += _w0 * (_w1 * err - _z1);
_z2 += _w0 * (_z1 - _z2);

Low-pass filter the error with a pair of first-order low-pass filters.  (_z1 is err filtered, _z2 is _z1 filtered).  Note that for this to work, _z1 and _z2 need to be persistent.  In C/C++, this means that they should be at least at file scope -- and in my humble opinion, this means they should have longer, descriptive names, even to the point of declareing a struct called "filter_state" -- it never hurts to hand your reader a clue-by-four.
I'm not going to dredge through all the loop analysis, but the cutoff frequencies of these filters should be significantly higher than the target loop bandwidth -- and the loop bandwidth is determined, in part, by the gain of the timing error detection process.  So you can't just slap this into your code willy-nilly.
_z3 += _w2 * _z2;
This is an integrator.
_rcorr = 1 - _z2 - _z3;

if (_rcorr > 1.05) _rcorr = 1.05;
if (_rcorr < 0.95) _rcorr = 0.95;

This does the final calculation of the command variable to the resampler.  Note, however, that we have a controller with an integrator, and we're limiting the output of the integrator, but we're not limiting the state of the integrator.  If you implement this code, search on "integrator anti-windup" for why this is a mistake, and why _z3 should be limited in magnitude, somewhere between $\pm 0$ and $\pm 0.1$, depending on your desired recovery characteristics from a large excursion.
_resamp.set_rratio (_rcorr);
This is actually setting the resample ratio, which in turn affects the value of err in the next iteration of the loop.


What is the advantage of Wiener filter for noise reduction of a 1D signal? I have a 1D signal, acquired by an accelerator sensor that measure the vibration of a structure. What is the advantage of using a Wiener filter for noise reduction compared to the other (more classic) filtering methods?
I have a signal + noise $u[n]= s[n] +w[n]$ measurement and a noise only $w[n]$ measurement, taken at rest.

The Wiener filter considers statistical behaviours of the noise and the signal, and thus, (theoretically) achieves optimum separation of them for a class of signals and systems, which is not the case for more classical approaches.
Wiener filter frequency response is such that, at those frequencies where noise power is dominant (a.k.a. low SNR), the gain is reduced, and the output is suppressed; causing noise-reduction. Whereas for frequencies at which the signal power is dominant (a.k.a. high SNR), then the gain is closer to one, and output is closer to the input.
Thus, to determine the gain of the Wiener filter, at a particular frequency, you must know the Power Spectral Densities (PSDs) of the noise and the desired signal.
Typically, you won't know these in advance very well; but they can be estimated from available data. The result is a departure from the optimum (ideal) performance.


CFO estimation and compensation I am reading about CFO estimation in OFDM system, assume for example that we estimate it based on CP as mentioned HERE and the estimated CFO for the received signal $r(t)$ is $β$. My question, how now can we compensate that CFO on the time-domain signal $r(t)$ ?

As with any signal: multiply with an appropriate $e^{-j2\pi \beta n}$.
(That, by the way, directly results from the last line of the answer you've linked to.)


About FM signal's proving My text book (Communication Systems analysis and design P.334, 335) have a FM signal is

  $s_{FM}(t) = A cos[2\pi[f_c+ks(t)]t+\phi]$ where k is a constant and $\phi$ is the arbitrary phase angles

then the textbook suggest

  while in FM, instantaneous frequency minus carrier frequency is a linear function of s(t)

So, I try to prove the $s_{FM}(t)$ is a FM signal

The instantaneous frequency is $f(t)= \frac {d\psi}{dt}=2\pi f_c+2 \pi s'(t) t +2 \pi s(t)$

and the instantaneous carrier ferquency is $f_{carrier} (t)=\frac{d\psi_{carrier}}{dt}=\frac {d}{dt} [2 \pi f_c t + \theta]=2 \pi f_c$. Assume the carrier signal is $Asin(2 \pi f_c t + \theta)$

So, $f(t)-f_{carrier}(t)=2 \pi s'(t) t +2 \pi s(t)$.

And my problem is why the text book suggest that signal is a FM signal? Clearly, I cannot prove the signal is linear.

As you have correctly derived, the difference of carrier and instantaneous frequency is

$$
g(s(t)) = \omega(t) - \omega_\mathrm{carrier} = 2\pi k s(t) + 2\pi k s(t)'t
$$

I will omit $t$ for convenience. If $g(s)$ is a linear function of $s$ it must fulfill the following condition:

$$
g(c(a + b)) = cg(a) + cg(b)
$$
So here:
$$
g(c(a+b)) = 2\pi k c(a + b) + 2\pi k \left[\frac{d}{dt}(c(a + b))\right]t\\
= 2\pi k c a + 2\pi k c b + 2\pi k c a' t + 2\pi k c b' t\\
= c(2\pi k a + 2\pi k a' t) + c(2\pi k b + 2\pi k b' t)\\
= cg(a) + cg(b)
$$
Thus, $g(s(t))$ is indeed a linear function of $s(t)$.


Detecting level in a Pulse Train Assume I have a square wave pulse train, where each high lasts for H and each low lasts for L seconds. I sample this waveform and have the sequence of numbers.

Question is what should be the sampling rate if I just want to know what level I am in(High or Low)?

I know by Nyquist I can find the sampling rate for full-reconstruction of the pulse train, but I don't need the full reconstruction. I just want to be able to:

1. Don't miss any High or Low

2. Detect all Lows and Highs

The minimum level time (i.e. minimum of L and H) must be greater than the sampling period. This way, there is no possibility for a pulse to be missed.

So, $T_s < \rm{min}(L,H)$.


Why level of noise can be magnified twice through each numerical differentiation? I was reading a paper and saw this is mentioned there, but I cannot figure out how this can analytically be proven?

Simply put: take an independent identically distributed  Gaussian noise (one observation in blue, left) and its gradient  (in green, right). The amplitude is roughly multiplied by the norm of the gradient filter (here $h=[1\,-1]$), which is $\sqrt{2}$, which you can see from the picture.

Their average energy  $E$ thus differ by a factor of $|h\|^2=\sqrt{2}^2 = 2$. As said by Olli, a similar question was answered in 2nd order edge detectors more susceptible to noise?. The part "can be magnified twice" is important, since this does not happen for all noises or all difference filters.

More generally, a random signal $x(t)$  passing  through a time-invariant filter with frequency response $H(f)$ produces a signal $y(t)$ with  PSD:

$$ S_y(f) = S_x(f) |H(f)|^2\,. $$

Since the Fourier domain is energy-preserving, you get that

$$\|H\|^2\ =|h\|^2=2\,.$$

This is a reason why people have developed noise-robust differentiation, often combined with smoothing (eg Savitzky-Golay filters). See examples at
Differentiation.

Note that in the continuous case, you should take care of non integrable processes, see for instance Variance of white Gaussian noise.


Re-recording audio vs. convoluting old audio data - differences? I want to transform a studio audio dataset into a smartphone audio dataset (used for a speech recognition application). Would there be any difference in the resulting audio files using one of the following two aproaches?

Replay all audio files and record with smartphone.
Convolute files with smartphone microphone impulse response, assuming the impulse response was determined using the same setup as
    in "1.".

Would there be any difference in the resulting audio files using one of the following two aproaches?

There would be a difference - probably negligible - which would represent the non-linearity or time-dependence in the system.


What is meant by product of two images? I would like to compute the product of two images. I tried with the code below. But the result is not like I expected. Can you tell me where is the problem in my code?

https://ibb.co/Njc8w3b

for (int i = 0; i < ImageWidth; i++)
{
    for (int j = 0; j < ImageHeight; j++)
    {
        PixelValueFirstImage = Giris1.GetPixel(i, j);
        PixelValueSecondImage = Giris2.GetPixel(i, j);
        Greyscale1 = Convert.ToInt16(PixelValueFirstImage.R * 0.3 + PixelValueFirstImage.G * 0.58 + PixelValueFirstImage.B * 0.12); //Converting RGB to grayscale
        Greyscale1 = Convert.ToDouble((Greyscale1 - 0) / (255 - 0)); // scaling between 0-1
        Greyscale = Convert.ToInt16(PixelValueSecondImage.R * 0.3 + PixelValueSecondImage.G * 0.58 + PixelValueSecondImageB * 0.12); //Converting RGB to grayscale
        Greyscale = Convert.ToDouble((Greyscale - 0) / (255 - 0)); // scaling between 0-1
        product = Greyscale * Greyscale1; // result of product
        product =(255*product)  ; //scaling between 0-255
        OutputImage.SetPixel(i, j, Color.FromArgb((int)product, (int)cproduct, (int)product);
    }
}

```

Multiplication of two signals is formally called as modulation; one of them is modulating the other.

In image processing most typically this will be called masking, when one of the images is a special selective mask for pixels.

However, when neither of the images is a mask, then their multiplication will blend (mix) them into each other according to their bright and dark regions...


What is differences between structured image and textured image I have been working on image processing. But there is one thing that has made me confused. What is difference between structured image and textured image. In my idea, both are the same. But recently, I read a paper that claimed that they are different. Somebody let me know what diffrerences they have?

Imho there is no "mathematical" description of either terms.

Both (texture and stucture) can be viewed as a repeating pattern. In the case of texture, the pattern might get a more random perturbation (think at an image of sand or grass).

In the case of a structured image you might get less random and more systematical perturbation (Think of a roof picture with tiles) Sometimes more info can be extracted from structured images. Imagine that you can actually get some 3D info by examining the way the pattern varies along one direction or another.

Of course I could be completely wrong and the "structure image" has a concrete definition for the author that you quote.

Edit:

I found two examples where 3D info is extracted from textural information. Chapter 9.5 "Shape from Texture" from "Introductory Techniques for 3-D Computer Vision" by Emanuele Trucco, Alessandro Verri (1998 ed) and Chapter 12.4 3D Cues Available in 2D Images from "Computer Vision" by L. Shapiro and G. Stockman (2000 ed)

There, you have the idea of deterministic texture vs stochastic texture. Also read Jedi Commymullah answer for some other insights.

Images are free:

http://pixabay.com/en/roof-roofing-house-roof-186391/
http://pixabay.com/en/grass-grassy-stalks-green-275986/


Scipy butter filter - remove DC offset w/ high sample rate I am having a hard time figuring out how to employ a high pass filter to remove the DC offset of my data signal with the "scipy butter" function because my sample rate is quite high.  The crux of my issue seems to be that my sample rate Fs is very high and my cut off frequency is very low (ideally just the DC offset component, but could be 1 to 10Hz with a sample rate of 125ksps).
Here is my example code:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# inital parameters
Fs = 125000         # sample rate (if changed to 125sps looks better)
n = 9               # filter order
Fc = 10             # cut off freq (Hz)
Nyq = Fs/2.         # Nyquist freq (1/2 Fs)
Fcc = float(Fc)/float(Nyq)  # normalized cut off freq

# find butterworth xfer functions
sos = signal.butter(n, Fcc, btype='hp', output='sos', analog=False)
# find freq/mag
w, h = signal.sosfreqz(sos, worN=1500, fs=Fs)
#convert to log scale
h_db = 20*np.log10(np.abs(h))

# plot filter response
fig1, ax1 = plt.subplots(1, 1)
ax1.plot(0.5*Fs*(w/np.pi), h_db)
ax1.set_title('db log')
plt.show()

I think this is working like it is supposed to because dropping the sample rate down to 125 instead of 125,000 yields a response that looks reasonable.  I think I'm missing something fundamental here with the Fc vs Fs.
Here is the response plot:

Here is a screenshot of the data I'm trying to apply this filter to:

As you can see there is a DC offset I'd like to filter out of my data file with a high pass filter.  Right now I'm just trying to get the filter response to behave so I'm not worried about processing the data yet - I just added the data file screenshot for reference.  I don't need anything fancy.  So my questions are the following:

Why do I have ringing and such an aggressive response?
How does the low cut off frequency and high sample rate affect the filter response?
Why did Scipy go from the b,a approach in the butter() filter to the "sos"?  I"m curious why they switch to this new sos appraoch.  My understanding is that it is more accurate and/or prevents errors in the b,a approach.  I'm unsure if this is just bolting second order transfer functions together though cascading simpler filters or what.
Also, you can notice (if you print(w) to see my freq axis) that my frequencies seem to run from 0 to 2*pi.  Is that normal or should it be from 0 to pi?
If question #4 is is "yes this should be 0 to 2*pi" than should line #23 in script:

actually be
ax1.plot(0.5*Fs*(w/(2*np.pi)), h_db)

so that my frequency range in the final result/plot runs from 0 to 1/2*Fs?
Thanks for any help!

Why do I have ringing and such an aggressive response?

Because you designed a really aggressive filter, which is borderline unstable. All the poles are less than 1/1000 away from the unit circle.

How does the low cut off frequency and high sample rate affect the filter response?

The lower the cutoff, the higher the order, and the higher the sample rate, the more it will ring (in terms of samples).

Why did Scipy go from the b,a approach in the butter() filter to the "sos"?

Because in transfer function representation the filter would actually BE unstable, i.e. the numerical precision of 64-bit floating point would NOT be enough to run this filter as a high order polynomial.

I"m curious why they switch to this new sos appraoch. My understanding is that it is more accurate and/or prevents errors in the b,a approach.

NO, NO, NO. Transfer function representation is the worst and should be avoided in favor of poles, zeros or second order section. Once you are in transfer function, you often can't get back to poles and zeros. It requires finding the roots of a high order polynomial which is a numerically ill defined problem.

Also, you can notice (if you print(w) to see my freq axis) that my frequencies seem to run from 0 to 2*pi. Is that normal or should it be from 0 to pi?

If you use a function, it's always a good idea to read the documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfreqz.html#scipy.signal.sosfreqz
Since you pass in the sample rate as an argument, the function already returns the frequency axis. By your attempt to normalize it again, you get complete non-sense. The dead give away here is the scale factor $10^9$ on your frequency axis. Your current axis tops out at 1.25 GHz.

If question #4 is is "yes this should be 0 to 2*pi" than should line #23 in script:

It should simply be ax1.plot(w, h_db)
The question you didn't ask:

What I should be doing ?

What @Ben said: use a first order filter AND either subtract the mean from the signal or seed the filter state with the signal mean.


Counting human movement of exercise using ROI or pre-set passing line with OpenCV I'm currently doing my final year Degree in computing studies.

The project I'm working on is a video of a webcam or computer webcam to have 2 pre-defined ROI (Region of Interest) on top and bottom of the video. It is to detect a person when he does a simple push up whether he reaches the area of both top and bottom ROI by detecting the person's head when he enters the ROI. I totally have very less experience in programming and really need guide and help.

During my research I've come across a journal/research at: Pattern Recognition Letters 33 (2012) 438–445. I've downloaded and gone thru it but can't seem to find its code or truly understand it. I got access with my University account at sciencedirect.com.

2nd idea I had was after viewing a video in Youtube: Traffic Counting based on OpenCV... I was very interested by that video. I was thinking if I could set 2 lines like he did and detect simple human motion up and down doing a simple accurate pushup exercise I would be so much happy but I can't seem to get contact or ask for any information.

I really need advice and help because I am really lost.

This might help partially -

You can, possibly track the person's head from his hair color(usually black), so can be distinguished from the background. This tracking just needs a simple RGB Filter like ( r<val2 && g<val4 && b<val3). Then, you can calibrate the program to respond when the centroid of that hair-region crosses your desired lines/boundaries.


How to convert text into wave I am trying to transfer text data via sound wave between two android phones. I am a newbie and would like to know how to proceed with it. I am stuck at the very first step on how to modulate text data as FSK or PCM raw data ? As I shall be implementing it on Android any Java pointers will be of great help

Text is just a string of bytes, and each byte is just a string of bits (typically 7 or 8), so the question boils down to: how to I send and receive a string of 1s and 0s using audio ?

You might consider using a simple AFSK scheme, where e.g. a 1 is represented by 2 cycles of a 2400 Hz sine wave and 0 is represented by 1 cycle of a 1200 Hz sine wave. So each bit period is the same length (8.14 ms) and the frequency for that bit period determines whether it's a 1 or a 0. Typically you would use a software UART to frame each character and handle the problem of synchronization between transmitter and receiver. This would allow you to send 1 ASCII character in around 10 bit periods. The bit rate is 1200 bps so you would be able to send around 120 characters of uncompressed text per second.


Spectral Inversion in Zero IF Receivers Is spectral inversion possible in zero IF receivers? If yes, then can somebody explain with an example how that happens and what can be done to correct it.
The following problem can be taken as reference:
A 40Mhz baseband signal is upconverted to fc of 1 GHz and transmitted.

Short Answer
In order to do a complete Zero-IF transmit and receive without inversion, the operations must be done as follows (or otherwise have the signs equally flipped on Q for both transmit and receive):
Zero-IF (Direct Conversion) Transmit:
$$I_{RF}(t) = I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)$$
Zero-IF Receive:
$$I(t) = \text{LPF}\{I_{RF}(t)\cos(\omega_{LO}t)\}$$
$$Q(t) = \text{LPF}\{-I_{RF}(t)\sin(\omega_{LO}t)\}$$
Yes spectral inversion can easily occur by having a sign change on Q on one side and not the other, showing that the way to correct it is to simply change the sign of the Q output.
Where:
$I_{RF}(t)$: Transmitter RF signal
$I_s(t)$: Baseband waveform I input at transmitter
$Q_s(t)$: Baseband waveform Q input at transmitter
$\cos(\omega_{LO}t)$: I LO
$\sin(\omega_{LO}t)$: Q LO
$\text{LPF}\{\}$: represents a low pass filtering operation.
Detailed Explanations
Consider the plots below that demonstrate exactly how this occurs:
For context first consider the typical Zero-IF Receiver architecture where we have a real signal that gets multiplied with a complex Local Oscilator, as implemented with two multipliers, one for the real or "in-phase" component as "I" and the other for the imaginary or "quadrature" component as "Q".  In the diagram the LO is shown to be fed into a quadature splitter (indicated by 90°) which produces two sinusoidal outputs in quadrature. The two outputs form a single complex source as real and imaginary components. One is designated to be "cosine" while the other is "sine" with "cosine" typically representing the real component. This is consistent with the Euler relationships listed further below.

Below shows the spectrums during this downconversion process. Note the Euler identities in the creation of the complex LO:
$$e^{+j\omega t} = \cos(\omega t) + j\sin(\omega t)$$
$$e^{-j\omega t} = \cos(\omega t) - j\sin(\omega t)$$
Multiplication in the time domain is convolution in the frequency domain, so the result of multiplying the real input with the complex LO (which is a single sided impulse in frequency) is to just shift the spectrum accordingly. So when we multiply by $e^{-j\omega t}$ the spectrum slides to the left as a result (and then low pass filter to get the resulting baseband signal), and when we multiply by $e^{+j\omega t}$, the spectrum slides to the right. Notice how the first case results in an inverted spectrum.

To reverse this inversion, we could simply change the sign of the sine (Q) output of the LO resulting in a negative frequency LO instead of a positive frequency. Alternatively and equivalently change the sign of the Q output of the receiver. To see this simply, consider a single baseband positive tone which in the time domain can be represented as a phasor rotating counter-clockwise on a complex IQ plane. If you negate Q (the imaginary component of the phasor), it would then rotate clock-wise instead, representing a single baseband negative tone of the same frequency.

The OP had additional questions specific to the transmitter side (direct-conversion transmitters for which I add the details below:
To explain the graphic we see how the transmitter structure can be viewed as the real component of a full complex multiplication; starting with the baseband Tx waveform on the top spectrum, the complex LO as $e^{+j\omega t}$ specifically to shift the spectrum to the right (for no inversion), resulting in what would be a complex RF signal on the third spectrum (which would require a real and imaginary datapath). By taking the real part of this, we get the conjugate symmetric real spectrum as on the bottom spectrum that can be fed to a single antenna.

To realize this, consider the complex baseband as $I_s+jQ_s$ multiplied with the complex LO as $e^{+j\omega_{LO}t}= \cos(\omega_{LO}t)+j\sin(\omega_{LO}t)$. The result of a full complex multiplication (requiring 4 real multipliers) would be:
$$(I_s(t)+jQ_s(t))\bigl(\cos(\omega_{LO}t)+j\sin(\omega_{LO}t)\bigr) $$
$$ = \bigl(I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)\bigr)+j\bigl(I_s(t)\sin(\omega_{LO}t)+Q_s(t)\cos(\omega_{LO}t)\bigr)$$
With the real part representing the transmitted RF signal specifically as:
$$I_{RF}(t) = I_s(t)\cos(\omega_{LO}t) - Q_s(t)\sin(\omega_{LO}t)$$
Notice the minus sign on the Q LO input as required to not have an inverted spectrum. If one were to make the mistake in using $+\sin(\omega_{LO}t)$ for Q then a spectral inversion would result! This could very well be a source of confusion since $e^{+j\omega t} = \cos(\omega t)+j\sin(\omega t)$, but one needs to review carefully the operation above when taking the real of the full complex multiplication to see why the negative sign is required to actually realize the positive complex LO!
Similarly with the Rx, in order to not introduce and inversion, we must multiply the received real RF signal by $e^{-j\omega_{LO}t}$ as shown earlier in the first spectrum plot on the left hand side. But in this case the RF signal is completely real, so we can show as in the product done out for the transmitter but in this case use $I_{RF}$ as the signal multiplied with the complex LO as $e^{-j\omega_{LO}t}= \cos(\omega_{LO}t)-j\sin(\omega_{LO}t)$ to get the complex baseband output (prior to low pass filtering). This results in:
$$(I_{RF}(t))\bigl(\cos(\omega_{LO}t)-j\sin(\omega_{LO}t)\bigr) $$
$$ = I_{RF}(t)\cos(\omega_{LO}t) -jI_{RF}(t)\sin(\omega_{LO}t)$$
With the final zero-IF receiver output as:
$$I(t) = \text{LPF}\{I_{RF}(t)\cos(\omega_{LO}t)\}$$
$$Q(t) = \text{LPF}\{-I_{RF}(t)\sin(\omega_{LO}t)\}$$


Poles and Zeroes - Z Domain 2nd order Butterworth transfer function:

$$H(z) = \frac{ 0.0201z^2 + 0.0402z + 0.0201 } { z^2 - 1.5609z + 0.6413 }$$

I'm a bit uncertain about the factorization at the denominator. Here's what I have so far:

\begin{eqnarray*}
H(z) &=& \frac{ 0.0201(z^2 +2z + 1)  }{ z^2 - 1.5609z + 0.6413 }\\
 &=& \frac{ 0.0201(z + 1)(z + 1)  }{ z(z - 1.5609) + 0.6413 }
\end{eqnarray*}
Zeros : $z + 1 = 0$, zero at $z = -1$.

How can I find the poles?

use the quadratic formula.  you know:

$$ a \, x^2 \ + \ b \, x \ + \ c \ = \ 0 $$

$\implies$

$$ x = \begin{cases}
 \frac{-b \, \pm \sqrt{b^2 - 4 a c}}{2\,a}  \quad \text{if } b^2 \ge 4 a c \\
  \\
 \frac{-b \, \pm j \sqrt{4 a c -b^2}}{2\,a}  \quad \text{if } b^2 < 4 a c\\
\end{cases}$$


Creating a 4th order Butterworth filter R I have sampled my data at 2000 Hz and want to generate a fourth order low-pass zero-lag filter with a cut-off frequency of 150 Hz.
I have put the code below that I have created but I have the following questions. Given that I have used the filtfilt function do I need to change the number to “2” so that it is a second-order and should I also double the cut-off frequency i.e change “W” to 0.3 from 0.15.
‘’’rawdata = data$Moment unfiltered
bf = butter(4,0.15,type = c("low"), plane = c("z"))
filtereddata = filtfilt(bf, data$Moment unfiltered)’’’
Thank you in advance!

do I need to change the number to “2” so that it is a second-order

Yes

and should I also double the cut-off frequency i.e change “W” to 0.3 from 0.15.

No
Strictly speaking, doubling up the filter will make the gain at the cutoff frequency -6dB instead of -3dB. So you have created a Linkwitz-Riley and not a Butterworth filter.
If you want a gain of exactly -3dB at 150 Hz you would need to increase the frequency slightly. You can work out the math by calculating the -1.5dB point of a Butterworth filter, but trial and error is probably quicker. 185Hz should get you fairly close.
In any case, the shape will still NOT be that of a Butterworth. Cascading two Butterworth filters does not result in another Butterworth filter. The zeros are in the right location but the poles are simply duplicated and are not Butterworth-distributed any more.


From time domain to frequency domain How is it possible to go from time domain amplitude values in floating point to frequencies?

I know the DFT is used, but what exactly is happening for it to be able to extract frequencies from amplitude values?

Is amplitude between -1 and 1 mapped from pressure levels? How are different frequencies detected in a block of amplitude values? What if two files have the exact same amplitude values but completely different frequencies?

The Mathematical structure of all Transforms can be understood by simply first seeing what a Fourier Series Representation is. In Continuous Time Domain, Fourier established that all periodic waveforms can be expressed as sum of Complex exponentials of some magnitude(which are the coefficients in a Fourier Series). The same can actually be applied for an aperiodic waveform assuming a periodic extension of it. So what the underlying Transform Relationships give us are nothing but the Analysis and Synthesis equations for what frequencies are present in a signal and what are their magnitudes. Rather, you can view these magnitudes as the Projections of the signal on the complex exponentials as to trying to figure out how much of a particular frequency is present in a signal. (You might want to remember the Analysis and Synthesis equations to see what I've said here).

Now for a certain set of functions(L2 functions or rather Limited Energy functions $\int_{-\infty}^{\infty}|x^{2}(t)|dt$ being finite and a set of conditions called Dirchlet's conditions being satisfied) have a special and easily identifiable Analysis equation called the Fourier Transform which comes straight from Fourier Representation of Series coefficients. Now in real world, we rarely encounter functions that violate these properties. If a function were to blow up to infinity in Energy we can never use it anyway.

Caution must be exercised because the Fourier representations involve a lot of intricacies which are rather ignored in daily life problem solving. One might learn a lot by looking at Existence, Measurability, Riemann Integrability etc. of functions and their energies.

Once you've tackled the problem of continuous time(both in terms of Series and Transforms) Discrete is just extrapolations right from Nyquist Shannon's sampling theorem. Now what exactly is a digital signal? It is just the sampled version of it. How quickly sampled is what determines the frequency content that you get  in the Sampled version. If you look at the formula of a DFT it is nothing but the Analysis equation of Fourier Series represented in $0,1...N-1$ points which are called "Bins". It is nothing but the Discrete version of the Fourier Series(think of Integral in the Analysis equation being replaced by Summation and the $\exp(-j\omega t)$ term in the transform being written as $exp(-2*pi*j*k/N)$ where k represents the number of points that you have from the sampled signal and $2*pi/N$ is nothing but the analogous term in the series which is the Fundamental Angular Frequency. I hope this explains much about how exactly DFT is able to get the Frequencies from a Time domain wave.

And for the last part, for 2 files to have exactly the same amplitude values at exactly the same points of time, it is impossible for them to differ in their frequency content because the Math behind it sees just the values. Maybe I'm misinterpreting your question or your phrasing is not proper.


What should the time-shift be when implementing a continuous wavelet transform on a computer? I'm currently researching implementation methods of the Continuous Wavelet Transform(CWT). On paper, the CWT produces infinitely many outputs on a finite signal since the scaling and shifting parameters are continuous.
When discretizing the CWT for a computer application, it is common to limit the scaling based on the bandwidth of the signal and the number of desired voices per octave.
However, I'm confused as to how the time-shifting should be limited. Since we're operating on a discrete input signal, we could technically just time-shift each wavelet by 1 sample to obtain the highest resolution. Yet, this seems wasteful for low-frequency convolutions.
It seems like the time-shift should be a function of scale, such that the time-shift is smaller for high-frequencies and larger for low-frequencies.
How is the CWT time shift typically determined when implementing the algorithm on a computer?

It's indeed wasteful for some signals, but unfortunately not much to do there except downsample afterwards. A fully rectangular representation enables reconstruction and analysis operations that aren't otherwise possible - further, doing 1 shift per sample can be much faster and take less (compute-) memory than e.g. skipping 1, via FFT convolution.
I provide a naive implementation with visuals here, also comparing against Python library implementations. Further, just today I've pushed out a ssqueezepy release, which has the best open-source CWT implementation I know of (I've yet to see if it beats MATLAB's) - can inspect its code.
If you insist on discretizing per-scale, I can only suggest a trial & error approach, as results will vary based on data, wavelet, and wavelet parameters used. You can then plot as a single 2D image by duplicating adjacent time-axis values for higher scales (but that's much effort and guesswork for little or no gain in speed).


Observability for Kalman Filtering? I wanted to know how observability of a stochastic state space system affects the performance of a Kalman Filter. Do we check for the usual observability matrix involving $\mathbf{C}$ (observation matrix) and $\mathbf{A}$ (state transition matrix) or there is a newer notion of stochastic observability? In case of non-observability, does the error between true and estimated states go out of bound ?

The only quick reference I can find is this one, which states:


Calculating a delayed response of a system Suppose I have a system:

$$y[n] = x[n] + x[n-1] + x[0]$$

To see if this system is time-variant, we can compare $H\{x[n-1]\}$ to $y[n-1]$ and check that they are different.

I think the delayed (by 1 unit) response to $x[n]$ is:

$$y[n-1] = x[n-1] + x[n-2] + x[-1]$$

Nevertheless, I'm having an issue understanding why. How does the $x[0]$ term (from $y[n]$) become $x[-1]$ while calculating $y[n-1]$?

Mathematically, to calculate $y[n-1]$, I thought it would be reasonable to replace $n$ by $(n-1)$ in the equation, which would result in $y[n-1] = x[n-1] + x[n-2] + x[0]$. I mean, if $y[n]$ response contains a constant tern $x[0]$, I expect this value to remain in $y[n]$ independent of the $n$ index.
Why that's not the case ?

I also have a problem with a different system:
$$y[n] = x[2n]$$

I don't understand why this system presents a delayed response:
$$y[n-1] = x[2n - 1]$$

(instead of $y[n-1]=x[2n-2]$) and why its response to the shifted (by 1 unit) input is $H\{x[2n - 2]\}$ instead of $H\{x[2n -1]\}$.

Looks like I have doubts about how to calculate the delayed response. Can you point out to problems in my examples to help clarify it?

The answer @Juancho liked to is indeed very useful. To me it looks like you don't understand time invariance fully, and that answer sould clarify things.

But, concerning your questions of getting the delayed response of the systems:

With the first system, $y[n] = x[n] + x[n-1] + x[0]$, you are right. The constant term will not change, since it is constantly $x[0]$, and thus not dependent of $n$. The correct delayed response is:

$$y[n-1] = x[n-1] + x[n-2] + x[0]$$

You can calculate that by substituting $n$ by $m-1$ (or $n-1$, it's the same which letter you use, but it's simpler for me to use a different letter while substituting). Now, whit $n=m-1$, you have:

$$y[n] = y[m-1] = x[(m-1)] + x[(m-1)-1] + x[0] = x[m-1] + x[m-2] + x[0]$$

When you get the final result (when you can't simplify any more), you can change the variable name from $m$ to $n$ again and get the first formula above.

With the second system, you employ the same logic to get the delayed response. Starting from:

$$y[n] = x[2n]$$

To get a delayed response you must substitute $n$ with $m-1$ and you get:

$$n = m-1$$
$$y[n] = y[m-1] = x[2 \times (m-1)] = x[2m - 2]$$

and now you can change the variable name again to get:

$$y[n-1] = x[2n - 2]$$.

These are explanations on how to calculate a delayed response (by 1). Looking at a response of a system delayed by 1 actually means that at time $n=5$, you will be looking at a response that you would normally get for $n=4$. The system is delayed, i.e. it behaves as the original, just with some time-delay.

To put it very simply, if I mark the delayed response by $y'[n]$, and the original by $y[n]$, this would hold:

$$y'[1] = y[0]$$
$$y'[2] = y[1]$$
$$y'[100] = y[99]$$

and more generally:

$$y'[n] = y[n-1]$$

For understanding time-invariance in systems, consult the linked question since that's precisely its topic.

To comment on the incorrect answers you proposed as delayed responses: they are actually responses to a delayed input.

What that means is that you will be feeding your original system with a delayed input. If your original input signal was $x[n]$, to get a response to a delayed input, you actually want to "feed" your system to a signal $x'[n] = x[n-1]$. Let's mark systems response to a delayed input with $y''[n]$ and calculate that for your first example:

$$ y''[n] = x'[n] + x'[n-1] + x'[0]$$

The system didn't change: It's still calculating the same operation: summing up the input at time $n$ with the input at time $n-1$ with the value of the input at time $0$. But, the input changed!. Substituting $x'[n] = x[n-1]$, you get:

$$ y''[n] = x[(n)-1] + x[(n-1)-1] + x[(0)-1] = x[n-1] + x[n-2] + x[-1]$$

To summarize, in order for a system to be time-invariant, the response of the system delayed by $k$ should be the same as the response of the system to the input delayed by $k$.

In this example, we were just examining delays by $1$. What proves that the system is not time-invariant in you case is the fact that $y'[n] \neq y''[n]$ because you have:

$$y'[n] = y[n-1] = x[n-1] + x[n-2] + x[0]$$
$$y''[n] = x'[n] + x'[n-1] + x'[0] = x[n-1] + x[n-2] + x[-1]$$


wavelength from peaks I have a simple 1sec signal of 1Hz added with 10Hz. I can detect all peaks (I use Python Scipy find_peaks) and I'm trying to deduce the 10Hz signal from those peaks. I can do it but it's not clean since the peaks spacing are not exactly the wavelength. Is there a formula that would adjust the wavelength according to the peaks height difference ? I'm trying to design a filter that removes the highest frequencies at each pass through the signal until there is none left.

and here's the Python code (without signal import or plotting):

_signal_build = np.zeros(len(_signal))
_sig_max, _ = find_peaks(_signal)
_sig_min, _ = find_peaks(-_signal)
for _current in range(0, len(_sig_max) -1):
      _wav_len = int(_sig_max[_current +1] - _sig_max[_current])
      _sine_start = _sig_max[_current] - round(_wav_len / 4)
      _sine_amp = 0.25*_signal[_sig_max[_current]] + 0.25*_signal[_sig_max[_current+1]] - 0.5*_signal[_sig_min[_current]]
      _sine = np.sin(np.linspace(np.pi, -np.pi + 2*np.pi/_wav_len, _wav_len)) * _sine_amp
      _signal_build[_sine_start:_sine_start + _wav_len] = _sine
      print(f"{_sig_max[_current]}-{_sig_max[_current+1]}, _wav_len: {_wav_len}, _sine_start: {_sine_start}, _sine_amp: {_sine_amp}")
_signal_filt = _signal - _signal_build

Thanks, cheers :)

Guillaume

In your comments you stated the purpose of this is to design a higher resolution transform than the FT offers. A consideration for the case of multiple frequencies that are not closely spaced is to use the adjacent bins in the DFT to interpolate a more precise frequency value. The simplest way to do this is to zero pad the time domain signal prior to computing the DFT, which will provide more samples of the Discrete-Time-Fourier-Transform (DTFT), revealing the frequency location of a localized single tone with much higher precision (the max value of the interpolated DFT will be the frequency location of a single tone).

Some precautions with this approach: Frequency resolution is theoretically the inverse of the time duration of the signal; for a non-windowed (or optionally stated using a rectangular window) each DFT bin in frequency represents the integrated energy under a Sinc shaped frequency response centered on that bin, with the first null of the main lobe of the Sinc function spaced one bin away. Thus we immediately see how this will not work to discern two closely spaced bins that may both be within the same main lobe (zero padding interpolates the frequency response but does NOT increase frequency resolution). For the case of one single frequency tone with zero padding the DFT result will be this Sinc function showing the additive correlation that would result to frequencies at any other locations under this curve, the peak of this will indeed be sufficiently close to the exact frequency location of the tone (within the precision of the number of samples we have provided with zero padding). However for two closely spaced frequencies, the interpolated result will appear as one Sinc function but will be the weighted and summed contribution of the two individual complex Sinc functions that each tone would give. Since the Sinc function shows relatively strong correlations even for tones spaced further away (the peak of the sidelobes of the Sinc only go down at rate $1/f$), it also shows how multiple tones can influence each other in the result depending on their relative locations. To combat this latter point, windowing should be used, which serves to increase the main lobe (increases the minimum distance two tones can be to discern each of them) while significantly decreasing sidelobes (so removes sensitivity to multiple tones further away).

Another precaution is that real signals appear as two tones in the DFT as demonstrated by Euler's relationship for a real cosine:

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

Where here each exponential frequency ($e^{j\omega t}$) Is what actually maps to an individual tone in the frequency domain.

So for the case of higher or lower frequency real tones or low number of samples these two tones that accompany a real cosine or sine function can also be too close to influence each other to not provide an accurate result.

This is demonstrated by the simple case of the DFT of 5 samples for a frequency that is exactly at 1.5 bins for both a cosine and single exponential frequency where the blue stem plot is the result for the DFT and the interpolated line plot is the samples of the DTFT from zero-padding:

And repeated for the case of a real tone using a cosine of the same frequency:

Here we see the influence of the two closely spaced tones (in terms of number of bins) since one would be located at 1.5 bins and the other at 3.5 bins so there respective correlations influence the actual location of the peak (the frequency locations are both 2 bins apart as well as 3 bins apart due to the cyclical nature of the DFT).  This example is the case of both a low real frequency in terms of number of bins and low number of samples. Even increasing the number of samples won't decrease the influence of the two tones as the two frequency locations will still be 3 bins apart regardless of number of time domain samples, so even windowing would not help in this case.

With that precaution hopefully better understood, here is a practical example of the approach with the case of 50 time domain samples and two real frequency tones at f1 = 12.3 bins and f2 = 20.5 bins, for the case of no windowing and windowing with a Kaiser window using $\beta=8$ (kaiser(50,8)). Zero padding was out to 10,000 samples, so providing for a resolution of the max value of 0.005 bins. With that the actual max values bin locations in the DTFT for both cases was:

No window: 12.255, 20.485

Kaiser window:12.300, 20.500

Also @CedronDawg has worked out the formulas for precisely computing the frequency of single tones which perhaps can be expanded to your case. See his blog post here: https://www.dsprelated.com/showarticle/1284.php


BPSK and QPSK demodulation computational complexity I wanted to ask about the real operations of BPSK and QPSK demodulation. How can we calculate it?
I mean for example if I have a received symbol +0.6, then I need to demodulate it to be 0 or 1. what's the computational complexity required for that. same for QPSK.

To demodulate BPSK and QPSK as described by the OP is to simply take the sign of the waveform (decision), which has trivial computational complexity. The real complexity lies in everything that occurs prior to decision in a practical receiver: There is typically an equalizer to remove multi-path distortion, a matched filter to optimize the signal to noise ratio at the decision locations (in the presence of white noise), an automatic gain control (AGC) loop to remove long term amplitude variation, a timing recovery loop to remove phase and frequency offsets in the timing clock (for the proper sample locations of the final decisions), and a carrier recovery loop to remove residual phase and frequency offset of the carrier frequency due to differences in the Tx and Rx clocks and Doppler due to possible motion between transmitter and receiver.
For instance, the block diagram below shows a practical receiver carrier recovery loop that also does the final decisions for QPSK; where for QPSK each decision would be a simple threshold detector on I and Q (notice the relative computational complexity of the decision block versus the rest of the loop, and this is just for carrier recovery and does not include the AGC and timing recovery operations described above that would be typical).


Multi-channel audio upsampling interpolation I have a four-channel audio signal from a microphone tetrahedral array. I wish to upsample it from 48 kHz to 240 kHz.

Is there a preferred interpolation method for audio? Does cubic interpolation (or any other) have any advantages over linear for the specific case of audio?

Assuming I am using cubic interpolation, do I interpolate each channel separately or is there any benefit in using a bicubic interpolation over all four channels?

Does cubic interpolation (or any other) have any advantages over linear for the specific case of audio?

You'd use neither for audio. The reason is simple: The signal models you typically assume for audio signals are very "Fourier-y", to say, they assume that sound is composed of weighted harmonic oscillations, and bandlimited in its nature.

Neither linear interpolation nor cubic interpolation respect that.

Instead, you'd use a resampler with a anti-imaging / anti-aliasing filter that is a good low-pass filter.

Let's take a step back:

When we have a signal that is discrete in time, i.e. has been sampled at a regular lattice of time instants, its spectrum is periodic – it repeats every $f_s$ (sampling freq.).

Now, of course, we rarely look at it this way, because we know that our sampling can only represent a bandwidth of $f_s/2$, we typically only draw the spectrum from 0 to $f_s/2$, for example:

S(f)
^
|---
|   \
|    \   ---
|     --/   \
|            \------\
+----------------------'---> f
0                    f_s/2

Now, the reality of it is that in fact, we know that for real-valued signals, the spectrum is symmetrical to $f=0$:

                     S(f)
                          ^
                       ---|---
                      /   |   \
               ---   /    |    \   ---
              /   \--     |     --/   \
      /------/            |            \------\
---'----------------------+----------------------'--->
-f_s2/2                   0                    f_s/2

But, due to the periodic nature of the spectrum of something that got multiplied with a "sampling instance impulse train", that thing repeats to both sides infinitely, but we only typically "see" the 1. Nyquist zone (marked by :)

       :                    S(f)                     :
       :                      ^                      :
       :                   ---|---                   :                   -------
  …    :                  /   |   \                  :                  /       \                      …
       :           ---   /    |    \   ---           :           ---   /         \   ---
       :          /   \--     |     --/   \          :          /   \--           --/   \
       :  /------/            |            \------\  :  /------/                         \------\
-------'----------------------+----------------------'---------------------------------------------'-->
     -f_s/2                   0                    f_s/2                                          f_s

When we increase the sample rate, we "just" increase the observational width. Just a random example:

                            S(f)
                              ^
                           ---|---                                       :------
  …                       /   |   \                                     /:      \                      …
                   ---   /    |    \   ---                       ---   / :       \   ---
                  /   \--     |     --/   \                     /   \--  :        --/   \
          /------/            |            \------\     /------/         :               \------\
-------'----------------------+----------------------'---------------------------------------------'-->
     -f_s/2                   0                    f_s/2              new f_s/2                   f_s

Try that! Take an audio file, let the tool of your liking show you its spectrum. Then, just insert a $0$ after every sample, save as a new audio file (python works very well for such experiments), and display its spectrum. You'll see the original audio (positive half of the) spectrum on the left side, and its mirror image on the right!

Now, to get rid of these images, you'd just low-pass filter to your original Nyquist bandwidth.

And that's really all a resampler does: change the sampling rate, and make sure repetitions and foldovers (aliases) don't appear in the output signal.

If you're upsampling by an integer factor $N$ (say, 48 kHz -> 192 kHz), then you just insert $N-1$ zeros after every input sample and then low-pass filter; it's really that simple.

In the ideal case, that filter would be a rectangle: Let through the original bandwidth unaltered, suppress everything not from there. A filter with a rectangular spectral shape has (infinite!) sinc shape in time domain, so that's what sinc interpolation is (and why it's pretty much as perfect as it gets).

Since that sinc is infinitely long, and your signal isn't, well, that's not really realizable. You can have a truncated sinc interpolation, however.

As a matter of fact, even that would be overkill: your original audio has low-pass characteristics, anyway! (simply because of the anti-aliasing filters that you invariably need before sampling the analog audio source; not to mention that high frequencies are inaudible, anyways.)

So, you'd simply go with a "good enough" low pass filter after inserting these zeros. That keeps the computational effort at bay, and also might be even better than the truncation of the sinc.

Now, what if your problem is decidedly not an integer interpolation? For example, 240000 / 44800 is definitely not an integer. So, what to do?

In this relatively benign case, I'd go for a rational resampler: First, we go up by an integer factor $N$, so that the resulting sampling rate is a multiple of the target sampling rate. We'd do the low-pass filtering as explained above, limiting the resulting signal to its original 44.8 kHz/2 bandwidth, and then apply a downsampling by $M$, i.e. anti-aliasing filtering it to the target 240 kHz/2 bandwidth, and then throwing out $M-1$ of $M$ samples.

It's really that easy!

In fact, we can simplify further: since the anti-imaging filter cuts off at 22.4 kHz, and the anti-aliasing filter only after 120 kHz, the latter is redundant, and can be eliminated, so that the overall structure of a rational resampler becomes:

Upsampling -> core filter -> downsampling

(in fact, we can even apply multirate processing and flip the order, greatly reducing effort, but that'd lead too far here.)

So, what are your rates here? For 44800 Hz in, 240000 Hz out, the least common multiple is 3360000 Hz = 3360 kHz, that's up by a factor of 75, low pass filter, and then down by 14. So, you'd need a 1/75 band lowpass filter. It's easy to design one using python or octave!


How to Model the Dispersion Which Is Introduced by an Electro Absorption Modulator (EML)? I need to model the dispersion which is introduced by an EML in Matlab. The parameters are:

Lambda: 1310nm

PAM4 power:
Pmax = 1.09 mW
Pmin = 0.13 mW

Bandwidth = 56Gbps

Dispersion is a frequency dependent phase (Instead of Linear phase).

Just transform your data into Fourier Transform and apply the phase per frequency.

Of course the data must be transformed back into time domain as @Deve noted.


Can PSNR > 48 used to validate the signals between two audio files? I have written some audio decoding code. I decoded an mp3 file and save as wav file.

In order to validate the code is working correctly, I compare the signal between the mp3 file and the decoded audio in wav file by using PSNR matlab (coding shown below).
The results obtained are shown in the table.

Signal1 = audio.mp3;
Signal2= audio.wav;
[R C]=size(Signal1);
err = sum((Signal1 -Signal2).^2)/(R*C);
MSE=sqrt(err);
MAXVAL=255;
PSNR = 20*log10(MAXVAL/MSE);

How can I comment on the result?

According to the post here, it says ‘If the reconstructed audio signal is exactly same as original signal then MSE =0. And if Max pixel value is 255 (8-bit representation), then the value of PSNR = 20*log(255) = 48dB.’

Since the PSNR of the songs are more than 48, can I conclude that the decoded audio signal in the wav file is valid?

If the reconstructed audio signal is exactly same as original signal then MSE =0. And if Max pixel value is 255 (8-bit representation), then the value of PSNR = 20*log(255) = 48dB

That statement is plain wrong. If the two signals are identical, then MSE = 0, then division by zero is undefinied and you don't get a PSNR at all, but "infinity".

  Since the PSNR of the songs are more than 48, can I conclude that the decoded audio signal in the wav file is valid?

PSNR is really not a very useful metric; you can get these values with a codec that simply broken, or one that sounds great.

Let's throw a bit of math in here: let's say we map our 8 bit to - to +1. Our original audio signal $S$ is well-scaled and its amplitude is pretty normally distributed¹, centered around 0, and has a variance $\sigma^2=\frac14$, so that clipping very rarely happens, i.e. $S\sim\mathcal N(0;\frac14)$.

Then, what's the PSNR of a codec that simply always produces "0"?

\begin{align}
\text{PSNR} &= 20\log_{10}\left(\frac{\text{Maxval}}{\sqrt{\text{MSE}}}\right)\\
&=20\log_{10}\left(\frac{\text{1}}{\sqrt{\text{MSE}}}\right)\\
&=10\log_{10}\left(\frac{\text{1}}{\sqrt{\text{MSE}}}\right)^2\\
&=10\log_{10}\left(\frac{\text{1}}{\text{MSE}}\right)\\
\end{align}

So, what's our MSE? It is

\begin{align}
\mathbb E\left((0-S)^2\right)&= \mathbb E\left( S^2\right)\\
&S \text{ is zero-mean}\\
&=\text{Var}\{S\}=\sigma^2\\
&=\frac14
\end{align}

Therefore,

\begin{align}
\text{PSNR} &= 10\log_{10}(4) \\
&\approx 6
\end{align}

So, what this would prove is that your codec is better than just producing the mean value of the information-theoretically worst possible continuous audio source – but that's about it.

¹ that assumption about the amplitude distribution is slightly wrong, but wrong in favor of your codec.


How is it possible to Filter out a person's voice out of 100 of other voices? So, I've just learned that the human voice is not a single sine wave its a bunch of unlimited Sine waves each having a different frequencies,

According to Wikipedia,

  The voice consists of sound made by a human being using the vocal
  folds for talking, singing, laughing, crying, screaming, etc. Its
  frequency ranges from about 60 to 7000 Hz.

So if the human voice is a composite signals, it contains various frequencies ranges between 60Hz to 7KHz.

Suppose if there is a group of people singing the same song all together, each person has its own ranges of voice frequencies,

For example,

If a person A has the following frequencies,
100Hz , 250HZ 6KHz, 10Hz, 87Hz, 52Hz, 2KHz.......

and a Person B has the following,
217Hz1, 11Hz, 12Hz, 2323Hz, 839Hz, 4KHz, 100Hz, 10Hz.....

there must be so many frequencies which are similar in both the person A & B, like in above example the frequencies 100Hz and 10Hz are common between two persons.

I was watching a TV Show name "Fringe" where they filter out the particular Man's voice from an audio file while there were other people voice present there too.

So how does exactly they filter someone's voice out of the voice of 100s of people if there are so many frequencies common among all of them does it have to do something with the amplitudes of person's frequencies ?

If the signal is recorded using just one microphone, you can use methods such as spectral subtraction. This method is more suitable for "constant" noise, like the noise from a fan or an idle engine. Other methods rely on statistics and perceptual models of speech. If the signal is recorded with several microphones, you can use blind source separation for separating the (speech) signals. As it stands today, you won't get perfect results. The typical end-result is always a trade off between "noise" and clarity of the speech signal of interest. More "noise" suppression --> more degradation of the signal of interest.


Harmonic distortion of sine signal. plotting resulting spectrum I'm trying to learn about how the spectrum of a sine signal changes when fed into a waveshaper and i've learnt a great deal from you guys here on stackexchange. Specifically how a transfer function such as $y = x + x^7$ adds the fundamental frequency of $x$ and then odd overtones up to the $7^{th}$ overtone.

Since i'm no mathematician and just now learning some dsp for an effect plugin i'm developing i need to be able to plot the spectrum of $y$ to learn more about the behavior how ploynomials affect this spectrum.

I could use Fourier tranform in matlab on $y$, but it feels like an unnecessary step.

I would like to be able to plot it like in this stackexchang question: Add odd/even harmonics to signal?

Thanks in advance.

Here is a way that you can approach this mathematically.  It's not the only way.  Another way to approach this is with Tchebyshev polynomials and I will leave that to someone else.

This approach requires understanding two ingredients.

an identity that is a consequence of Euler's formula:

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

the binomial theorem:

$$\begin{align}
 (a+b)^N &= \sum_{n=0}^N {N \choose n} \, a^{N-n} \, b^n \\
         &= \sum_{n=0}^N {N \choose n} \, a^n \, b^{N-n} \\
\end{align}$$
or
$$\begin{align}
 (a+b)^N &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^{N-n} \, b^n \\
         &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^n \, b^{N-n} \\
\end{align}$$

Now you might see where this is going.  So using the cosine identity (so i don't have to mess around with powers of $j$ in the sine identity):

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, (e^{j\omega t})^{k-n} \, (e^{-j\omega t})^n \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-n)\omega t} \, e^{-jn\omega t} \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\end{align}$$

It appears that we will have to deal with the two cases of even $k$ vs. odd $k$ differently.

For even $k \ge 2$:

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \frac{k!}{\tfrac{k}{2}! \, \tfrac{k}{2}!} + \sum_{n=\tfrac{k}{2}+1}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\
\\
 &= \frac{1}{2^k} \left(\frac{k!}{\left(\tfrac{k}{2}!\right)^2} +  \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{2 k!}{2^k} \sum_{n=1}^{\tfrac{k}{2}} \frac{1}{\left(\tfrac{k}{2}+n\right)! \, \left(\tfrac{k}{2}-n\right)!} \, \cos(2n\omega t) \\
\end{align}$$

Note that only DC and even-numbered harmonics come outa the even powers of $\cos^k(\omega t)$.

For odd $k \ge 1$:

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=\tfrac{k+1}{2}}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\
\\
 &= \frac{2 k!}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{1}{\left(\tfrac{k+1}{2}+n\right)! \, \left(\tfrac{k-1}{2}-n\right)!} \, \cos\big((2n+1)\omega t\big) \\
\end{align}$$

Note that only odd-numbered harmonics come outa the odd powers of $\cos^k(\omega t)$.

So now what if we pass this sinusoid:

$$ x(t) = A \, \cos(\omega t) $$

through this polynomial mapping function:

$$\begin{align}
 y(t) &= a_0 + a_1 x(t) + a_2 x^2(t) + a_3 x^3(t) + ... + a_{K-1} x^{K-1}(t) + a_K x^K(t) \\
      &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\
\end{align} $$

?

I believe we have to split the even-order terms from the odd-order terms in the summation

$$\begin{align}
 y(t) &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( x(t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( x(t) \big)^{2k+1} \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( A \cos(\omega t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( A \cos(\omega t) \big)^{2k+1} \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \cos^{2k}(\omega t) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \cos^{2k+1}(\omega t) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \left( \frac{1}{2^{2k}} \frac{(2k)!}{(k!)^2} + \frac{2 (2k)!}{2^{2k}} \sum_{n=1}^{k} \frac{1}{\left(k+n\right)! \, \left(k-n\right)!} \, \cos(2n\omega t) \right) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \left( \frac{2 (2k+1)!}{2^{2k+1}} \sum_{n=0}^{k} \frac{1}{(k+1+n)! \, (k-n)!} \, \cos\big((2n+1)\omega t\big) \right) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{n=1}^{k} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{n=0}^{k} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{n=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{n=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\
\end{align} $$

$\lfloor u \rfloor$ is the floor() function, which is the largest integer no larger than the argument $u$.

$$ \lfloor u \rfloor \le u  < \lfloor u \rfloor + 1 $$
or
$$  u-1 < \lfloor u \rfloor \le u $$

Also note that if $K$ is even
$$ \left\lfloor \tfrac{K-1}{2} \right\rfloor < \left\lfloor \tfrac{K}{2} \right\rfloor = \tfrac{K}{2} \quad \text{and} \quad \left\lfloor \tfrac{K-1}{2} \right\rfloor < \tfrac{K-1}{2} $$

If $K$ is odd
$$ \left\lfloor \tfrac{K}{2} \right\rfloor = \left\lfloor \tfrac{K-1}{2} \right\rfloor = \tfrac{K-1}{2} \quad \text{and} \quad \left\lfloor \tfrac{K}{2} \right\rfloor < \tfrac{K}{2} $$

$$\begin{align}
 y(t) &= \sum_{k=0}^{K} a_k \big( A \cos(\omega t) \big)^k \\
\\
     &= \sum_{n=0}^{K} b_n \cos(n \omega t) \\
\end{align}$$

where, apparently the DC component of the Fourier cosine series is:

$$ b_0 =  a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{(k!)^2} \left(\frac{A}{2}\right)^{2k}  $$

The coefficient for the even harmonic terms:

$$ b_{2n} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{\left(k+n\right)! \left(k-n\right)!} \left(\frac{A}{2}\right)^{2k}  $$

And the coefficient for the odd harmonic terms:

$$ b_{2n+1} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{(2k+1)!}{(k+1+n)!(k-n)!} \left(\frac{A}{2}\right)^{2k+1}  $$


Confusion Regarding Bi Linear Transform I was reading my book where the z-transform of a signal is derived to be ${1-e^{-2bT}z^{-1}}$ . Then it goes on to say that by applying the bilinear transform we can get

$$\frac{2(1+bT+(bT-1)z^{-1})}{T(1+z^{-1})}$$

I know a little bit about bilinear transform but this example is using it in a different context. Can someone please help me understand how the bilinear transform is being used here.

The bi linear transform is the transform from the Laplace Transform Domain to the Z Transform.

The Laplace Transform Domain is a regular plane.
This transform transforms vertical lines in the Laplace domain into circles in the Z Domain.
Hence the Fourier Vertical Line in Laplace Domain (The Y Vertical Lines) is transformed into the unit circle in the Z Transform.

I'm a little inaccurate here, since the whole left side of the plane is transformed into the inner part of the Z Plane and the right hand of the Laplace transform it transformed to the outer part of the Z Plane (Out of the unit circle).

This transform is used to derive Digital Filters from their Analog Laplace domain representation.


FM demodulation with arctan I am trying to demodulate a FM signal with a RTL-SDR and python using the arctan method. From here I know that
$$
\frac{d}{dt}\arctan(\frac{q(t)}{i(t)}) = \frac{i(t)\frac{d}{dt}q(t)-q(t)\frac{d}{dt}i(t)}{[i(t)]^2}
$$

and so far I have successfully been able to demodulate using the RHS and the following python code

def discrim(x):
    X=np.real(x) # X is the real part of the received signal
    Y=np.imag(x) # Y is the imaginary part of the received signal
    b=np.array([1, -1]) # filter coefficients for discrete derivative
    a=np.array([1, 0]) # filter coefficients for discrete derivative
    derY=signal.lfilter(b,a,Y) # derivative of Y,
    derX=signal.lfilter(b,a,X) # " X,
    disdata=(X*derY-Y*derX)/(X**2+Y**2)
    return disdata

However I would also like to get the same results using the LHS and naively assumed that python code as follows would do that

def discrim_bad(x):
    X=np.real(x) # X is the real part of the received signal
    Y=np.imag(x) # Y is the imaginary part of the received signal
    b=np.array([1, -1]) # filter coefficients for discrete derivative
    a=np.array([1, 0]) # filter coefficients for discrete derivative
    x = np.arctan2(Y,X)
    der=signal.lfilter(b,a,x) # derivative
    return der

Unfortunately this is not the case. I assume it has something to do with the discretization and most likely the discrete derivative approximations but I am not quite able to figure it out.

The $\arctan2()$ function has a jump at $\pm\pi$ which will cause problems when trying to compute a derivative.

Since you're using the first difference discrete derivative approximation to compute $d\phi/dt$, there's a simple way to avoid the jumps.

First note that your signal samples can be written polar form:

$$s[n] = i[n]+jq[n] = r[n]\cdot e^{\phi[n]}$$

so you can use multiplication by a complex conjugate to perfom angle subtraction before taking the $\arctan2()$:

$$s[n]\cdot s^*[n-1] = (i[n] +jq[n])(i[n-1]-jq[n-1]) = r[n]e^{\phi[n]}r[n-1]e^{-\phi[n-1]}=r[n] r[n-1] e^{\phi[n]-\phi[n-1]}$$

Thus, to get the discrete derivative you want, just take:

$$\mathrm{Arg}[(i[n]+jq[n])(i[n-1]-jq[n-1])]$$


Z-transform of an FIR filter QUESTION

Compute the Z-transform of $y[n] = x[n] + 2x[n-1]$. and find the poles and zeros.

I just bombed an interview where I couldn't do this (because I have no grounding in fundamentals and have worked all from cookbooks). I'm not too bummed because I would have been working at something I wasn't prepared for. You get one guess who would ask this kind of question in a first phone interview.

ANSWER (30 minutes post interview):

$x[n-1] = z^{-1}$

$1 + 2z^{-1}$.

Why?

Because the z-transform of an FIR is the coefficients of the impulse response multiplied by the delays.

If you know (or better: understand) the basic $\mathcal{Z}$-transform relation

$$x[n-k]\Longleftrightarrow z^{-k}X(z)$$

then from the given time domain equation you immediately get

$$Y(z)=X(z)+2X(z)z^{-1}=X(z)(1+2z^{-1})$$

Since the transfer function $H(z)$ is the ratio of $Y(z)$ and $X(z)$, you get

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


Non invertibility of system $y[n]=x[n]-x[n-1]$ using transform method? For the system to be invertible, we should have different outputs for different inputs.

In terms of constant functions say,

$$X_1[n]=3 \quad \forall n \in \mathbb{Z}$$ and
$$X_2[n]=4 \quad \forall n \in \mathbb{Z}$$

It is perfectly clear that the above system gives 0 for both the inputs hence can be concluded as non invertible. But how can this be explained in terms of transform method (Bilateral z transform or unilateral z transform )?

The impulse response of your system is $h[n] = \delta[n] - \delta[n-1]$ and its Z-transform is accordingly

$$H(z)=1-z^{-1}$$

Now, you can evaluate the frequency response of this system:

$$
|H(z)|_{z=\exp(j\omega)}=|1-\exp(j\omega)|
$$

you see, that $H(z)|_{z=\exp(j0)}=0$ (i.e. it vanishes for $\omega=0$). Now, in order for a filter to be invertible, its Z-Transform on the unit circle (i.e. when setting $z=\exp(j\omega)$ must not vanish. If it does not vanish, the inverse is given by $\frac{1}{H(z)}$. If it vanishes, the inverse does not exist.


Effect of Sampling frequency on RMS estimation I have a simple, but interesting question for you.
I have to compute the RMS (root mean square or standard deviation) value of a time signal which represents in my case a velocity mesasurement over time, of course sampled an stored in a computer. In your opinion, does the sampling frequency have a strong impact on the RMS value or it is just the TOTAL number of samples that matters the most?
Can you please, if possible, refer to some books where this topic is dealt with?

Thanks a lot, I'm looking forward to reading your answers.

Bye

E.

If the sampling rate meets the nyquist criteria (sampling rate more than twice the highest frequency in the signal) then it won't be a factor in the accuracy of the RMS value.

Since you are measuring velocities the nyquist criteria could be a problem if you don't know the frequencies involved.  The frequencies would be related to the acceleration - higher acceleration means higher frequencies.  I don't know how to tell you to determine the needed sampling rate give the expected acceleration.

If you can't determine the highest frequency (or are limited in your sampling frequency,) then you might run your sampling at the highest available rate, but filter the analog signal with a lowpass with a cutoff frequency lower than half of the sampling rate.  This will avoid the problems that arise when sampling a fast signal with a too slow sampling rate, at the expense of throwing away information if the sampling rate really is too slow.

The number of samples only comes into it when trying to calculate the RMS value of a repeating signal, line a sine wave.  In that case you have to make sure to get a least one full cycle of the lowest frequency.  More could help make up for inaccuracies in the measurement or reduce some of the noise.  In this case, the longer sampling period averages out any measurement errors and noise.

Summary: Make sure that your sampling is fast enough, that the samples are evenly spaced, and that you (obviously) have all the samples for the period of interest.


Why is the name "Phase vocoder"? I've read that phase vocoder is the STFT, but done in windows (i.e. using windowing).

Why is it called "phase vocoder"?

A vocoder is a type of synthesizer originally developed for speech (or the human voice, thus the "vo" prefix), even before the use of digital processors.  The Bell Labs Voder was demonstrated in 1939, and VoCoders were reportedly used for encrypted voice communication during WW2 (according to Wikipedia).

DFTs can be used to analyze which audio spectrum peaks are to be used to "code" the speech, and short STFT windows are used to get a fine enough time resolution to isolate typical phonic atoms over reasonably stationary time periods.

But the frequency bin centers of a short FFT are not spaced with sufficient frequency resolution to capture pitch, overtones, and modulations thereof, accurately enough to sound human.

But, by taking into account the FFT phase information, the difference in the phase of FFT result bins between adjacent offset windows can be used to adjust (or interpolate) the apparent or actual frequency of the voice codec synthesizers to more accurately model narrow-band frequency peaks than by just using the magnitude alone of the STFT FFT results.  Thus the "phase" prefix to the analysis-resynthesis vocoder process.

Note that if you slide two Windows apart, the amount of phase change of stationary spectrum should be linear with frequency.  If analysis shows a little more phase change, the frequency is higher, a little less phase change, lower.  And this phase difference will be measurable even for differences in frequency much less than the Fs/N bin spacing of each FFT.  Thus there is information gain in using succesive frames and phase, rather than single frames and magnitude only, for the frequency estimation of narrow-band spectral peaks.  Phase can either be stretched or compressed across frames during FFT resynthesis, or the frequency of higher resolution oscillators can be adjusted to match the estimated frequency during synthesis.

With increases in DSP performance, phase vocoders are now used for music as well as speech, but the "vo" prefix remains.


Zero Phase Filter: Determining Initial Conditions for Forward Backward Filtering Is anybody familiar with Gustafson's algorithm for minimizing transients in forward backward filtering [1]? I'm trying to implement it and my first guess was to check Matlab's filtfilt.m, since they are referencing the paper. In the Matlab function also a linear equation system is solved in order to find initial conditions zi that minimize the startup transients, but the relationship between reference and code is not obvious to me. The only lines of code regarding the minimization are (nfilt is the length of the coefficient vectors):

zi = ( eye(nfilt-1) - [-a(2:nfilt), [eye(nfilt-2); zeros(1,nfilt-2)]] ) \...
 ( b(2:nfilt) - b(1)*a(2:nfilt) );

Can anybody point me in the right direction on how those lines relate to the algorithm described in Gustafson's article?

[1]  Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992.

For anyone who is interested, i coincidentally found a paper describing the method implemented in matlab's filtfilt.m. A link to the paper is attached. At least to my understanding matlab's filtfilt.m doesn't implement the Gustafson algorithm.

Sadovsky, P.; Bartusek, K: Optimisation of the Transient Response of a Digital Filter, Radioengineering Vol. 9, No. 2, 2000


Getting 16 bit image statistics in Photoshop The histogram in Photoshop for a 16-bit image downsamples to 8 bits, showing a tonal range of 0 to 255. I appreciate to show tonal ranges for a 16 bit image a horizontal zoom function would be necessary.

However the mean, std deviation,... statistics also get downsampled to 8 bits with a corresponding loss of accuracy.

Is there a way of seeing the true mean, std deviation etc of a 16 bit image?

A lot of tools do this by default (even matlab's image processing tools for the most part). I'd load the image up in matlab and use the hist, std, mean, etc function by treating the image as a vector. Tools like EnvI also handle 16 bit images properly.

The matlab/octave code would be roughly as follows (assuming 16 bit greyscale png):

img=double(imread('image.png'));
mean(img(:))
std(img(:))
hist(img(:),number of bins)


Modelling transmitter and receiver impairments in communication system A paper I read studied the effect of transmitter impairments and receiver impairments by modeling both as additive white Gaussian noise. The noise for transmitter impairmenets are added before the antenna before transmission. The noise modelling receiver impairments is added as an additional white Gaussian noise to the thermal white Gaussian noise. A conclusion was made as follows

  The main reason for this is that the noise added at the transmitter (modeling Tx EVM) experiences the same frequency- and time-selective fading as the desired signal. In contrast, the noise added at the receiver (modeling Rx EVM) does not experience frequency-selective fading (only time-selective fading).

Does the quoted conclusion above make sense?

It makes sense to me. If the ideal transmitted signal is $s(t)$, the signal that is actually transmitted is $s(t)$ plus noise and distortion. This is the signal that suffers attenuation and fading, and to which an extra heaping of noise and distortion is added in the receiver. Modeling transmitter imperfections is not something that is included in most introductory textbooks, but it's very important to implement a system with good performance. You may see for example "RF System Design of Transceivers for Wireless Communications" by Qizheng Gu for more.

Having said that, modeling the transmitter imperfections as just Gaussian noise is a bit optimistic, but it may be a step in the right direction.


How do I check for linearity for the following piecewise-defined system? The problem at hand:

Where I'm currently stuck:

I'm not entirely sure about how to move on from this point, I'm trying to find the superposition of the responses of the two individual signals so I can later check if that's equal to the response of their (the two signals') superpositions. The intervals in each response involve different signals and thus I'm unable to find the superposition of responses.

First, you are not trying to find the linearity of a signal, but the linearity of a system defined by the transform above between input and output signals.
If you suspect your system is not linear (the different expressions depending on the sign of $x(t)$ is a big red flag), then you only need a counterexample rather than a general proof of linearity.
In this case, if you scale your input by a negative number, the output is not scaled accordingly. That is: $T{\alpha x} \neq \alpha T{x}$ if $\alpha$ is negative.


Are there any measurements or test for whether an audio signal is stationary? I need to divide an audio signal into frames in which the signal can be assumed to be stationary. How can I measure whether an audio signal is stationary in a frame?

I'm a DSP student and perhaps this is not the most complete answer. I have to answer the same question in a partial work and this is what I got so far:

As said by André pure stationary signal is a signal in which mean and variance does not change in time. In nature (real life, man interpreted music, etc.) this is very hard to find.

So you may:

Visualize the signal in an spectrogram. If the signal is stationary the signal will appear as horizontal bars of same intensity in the spectrogram:

Even in this signal, made from a filtered white noise and a tone the media and the variance changes a little in time.

So, split the segments of your signal under test to find the whether the media and variance doesn't change 'much'. How 'much' depends on you:

% test for stationary B

xnbf = buffer(signal,100);

mxn = mean(xnbf);

sxn = std(xnbf);

figure; errorbar(mxn, sxn, 'o'); title('Signal A')

Also, check this link, which includes a comprehensive list of methods to test whether a process is stationary or not:
https://quant.stackexchange.com/questions/2372/how-to-check-if-a-timeseries-is-stationary


DOA for bistatic radar Consider a circular antenna array (UCA) with $N$ antenna elements as part of the CW bistatic radar system.
For each element we receive a signal of length $L$: $x_n[k]$, $k=1...L$, $n = 1...N$.
Оne of the antenna elements is directed towards the transmitter.
Thus, we can find cross ambiguity functions (CAF) as cross-correlation functions between the reference channel and the remaining $N-1$ channels.
Further, in each CAF, we detect targets and determine the ranges to them.
My question is how to define directions on a target.
Standard explanations include beamforming. The way I understood this method: the samples from each of the channels, except for the reference for same time point, are used to form the directional diagram, i.e. from the $N \times L$ of the signal samples, $L$ diagrams are formed. Each directional diagram includes lobes directed to all targets as well as side lobes.
How to find out which direction corresponds to each target with its own range and how to exclude side lobes if their amplitude can be higher than the amplitude of the lobe directed at a target with a weak reflected signal?

You use complex symbols and at the same time for me it seems very strange, the concept you describe. Its hard to guess if you are a professor and actually want to know exactly what you are asking, or a student and want to know general concept about finding targets with antenna array.
If you want general answer, then your whole aproach is far from what is used today. You cant effectively detect direction from signal's strength alone. And having a working antenna array you are much better off using phase shifters that are supposed to be in that array, that will perform a convolution for you in hardware almost for free.
https://en.wikipedia.org/wiki/Convolution
https://en.wikipedia.org/wiki/Phase_shift_module
In simple terms, you change signal delay until signal added up from two antenas becomes most powerful. This value of delay will tell you the direction of the target. Do it for every pair and you can get the direction very precisely. This is what beamforming usually is.
If you ment exactly what you are asking, then it is very complex. You can collect the signal strength for each antena, then try to make a simulation, simulate what strength would each antena show if signal would be in a particular location. You use your array data, such as each antenna's direction and position to make this simulation. And you add each particular antena's radiation pattern. In simulation you are likely better off using monte-carlo method, calculate the collection of signal strengths for a few points. Select the best fitting point, calculate a few points nearby that point, with more compactly placed points. Proper antenas wont work at all for this, regardless of solution, even with the one i offered because they have symmetrical radiation pattern and will together only show how far the target's direction is from the array's center direction. Best antenas for this tasks are likely half-covered, to provide very unsymetrical radiation pattern, that will help with detecting the direction. I dont think anything like this exist today, and I've never seen anything like this in history. It could be beneficial only before phase shifter were popularised. I agree that it can be beneficial for very low tech radar, if you are actually making this, i would like to read it.


Filtering artefacts and filtering short signals I have a signal from EEG sensors and I try to denoise it from AC frequencies. For that reason, I estimated PSD of my signal and found that 50 Hz and 100 Hz are likely to represent noise. I constructed Butterworth filter of order four and got much clearer signal, but at the start ([0:150] segment) there is even more distortion. Why is it so? If it helps, I use lfilter from scipy.signal.

Besides, in the future I want to break a signal into smaller pieces (say, length of 100). I have tried denoising them already and it seems like such filters do not work on short segments. What can I do with this?

All real-time filtering (as opposed to post processing) wirh FIR and IIR filters will have start up transitions based on the state of the filter at start up. For optimum rejection of AC noise , Instead of a Butterworth Filter consider using an 2nd order IIR notch filter with the notch set at your AC frequency (such as 50 Hz). A design for this is further explained on this site and is simple and effective for this purpose, and happens to be demonstrated with a notch at 50 Hz. (Transfer function of second order notch filter).
To have shorter segments, partition into segments at the output of the filter so that the filter’s memory can be maximized to meet the filtering requirement as designed (a filter’s transition bandwidth between passing a frequency of interest and rejection is dependent on the span of its memory—- for tight rejection long memory is required.)
If the application is for very short durations at completely different instances in time, that can make no use of the prior immediate time domain signal, then such traditional FIR and IIR filtering techniques will not be suitable beyond what they are capable of for the given time span. That said, the IIR nulling approach with its rejection optimized at the interference frequency may be sufficient and simple.


Why do we need circular convolution? I am undergoing a course on DSP and we were recently introduced to the amazing concept of DFT.

But I have some apprehensions against it. DFT is supposed to be the methodology of fourier transform, as implemented in hardware or software. Now, if I want to convolve two signals' DFT, I could simply linearly convolve them. There would be no wraparound error, because all that I have to do to implement linear convolution is tell my computer that the array now flips and then multiplies and slides - i.e., linearly convolves.

Circular convolution was developed to avoid wraparound error; however, the computer does not know that DTFT is a periodic function, and, thus, shifting the DFT array will render me with a circular shift.

I know my opinion is rather radical, but I am sure that this is because I have a very abstract level acquaintance to this concept. Can someone enlighten me more, herein?

I think your confusion comes from the fact that when we talk about "circular convolution" we are not doing a convolution operation in the frequency domain.  Rather we are multiplying the two DFT arrays together, element by element.

When we switch back to the time domain via an inverse DFT, we have effectively done a circular convolution in the time domain.


Different definitions of unit step signal While I was learning signal theory, I have come across different definitions for the unit step signal. For example.
$u(t)=\begin{cases}
1 & t\geq0\\
0 & t<0\\
\end{cases}$
$u(t)=\begin{cases}
1 & t>0\\
0 & t<0\\
\frac{1}{2} & t=0\\
\end{cases}$
$u(t)=\begin{cases}
1 & t>0\\
0 & t\leq0\\
\end{cases}$
However, in the context of LTI systems the outputs are not different.
i.e LTI systems with these inputs produce same outputs. How can one justify this?

The Wikipedia article on Heaviside has a section that discusses the 3 different conventions on $H(0)$:

https://en.wikipedia.org/wiki/Heaviside_step_function

Essentially, the step function is used in probability theory as well but is called something different.

I believe that the "symmetric" $H(0)=\frac{1}{2}$ version is probably most suitable in Linear Systems Theory. In discussions on the Gibbs phenomena, the truncated Fourier series interpolates discontinuities at the mid-value which in this case is $\frac{1}{2}$.

Also, if you look at the article, there are a number of limits of continuous functions like Erf, that approach the step as some parameter approaches some value. If you were to set up a physical experiment, you would never be able to realize an ideal step function, you would need infinite bandwidth, so in reality, the best you could achieve would be something like a one of those continuous functions.  $H(0)=1/2$ is a vestige of physical realizability.


Why is the sound field intensity due to $K$ point sources given by $ I(p,\omega) = \sum_{k=1}^K \sigma_k^2(\omega) \delta(p - p_k)$? I am trying to understand the following piece of text. I am not used to dealing with sound intensity and power so I'm not familiar with the derivation of the formula $(*)$ below.

Statement: 1.

We assume there exist $K$ point sources in the far-field indexed by the letter $k$. Each propagates in the direction of the unit vector $p_k$. Then the signal coming from direction $p ∈ S$ is given by $x(p, ω, t) = \tilde x(p, ω)e^{jωt}$,where $\tilde x(p, ω)$ is the emitted sound signal by a source located at $p$ and frequency $ω$. The intensity of the sound field is then
$$
I(p,\omega) = \mathbb{E}[|x(p, ω, t)|^2] = \sum_{k=1}^K \sigma_k^2(\omega) \delta(p - p_k), \tag{$*$}
$$
where $\sigma_k^2(\omega)$ is the power of the $k$-the source and $\delta(p)$ is the Dirac delta function on the unit circle.

Question

How is $(*)$ derived? The formula for intensity on Wikipedia is
$$
I(r) = P/A(r),
$$
where $P$ is the power and $A(r)$ is the surface area of a sphere of radius $r$. The formula in $(*)$ doesn't seem to make any use of a sphere of radius $r$..and it also features Dirac deltas which aren't in the Wikipedia formula.

So what is the relationship between the Wikipedia formula and the formula $(*)$?

I believe the confusion comes from the fact that equation $(*)$ is referring to a received intensity (in units of W/m^2), whereas the second equation is referring to a radiated intensity.

The assumption behind equation $(*)$ is that the received power as a function of angle is the surface integral of the incident intensity over some aperture area $A \subset S$. Because the surface area element $dA$ has units of area, the result of the integral will be in units of power (W).

The Dirac delta distribution $\delta(p - p_k)$ expresses that the power density is coming from exactly one direction $p_k$ and nowhere else. In this context, the Dirac delta distribution possesses the following property:
$$
\int_A \delta(p - p_k) dA' = \begin{cases} 0 & p_k \notin A \\ 1 & p_k \in A \end{cases},
$$
which is analogous to a "density" where all the "mass" is concentrated at one point $p_k$.

To calculate the total received power $P$, you integrate the received intensity over the receive aperture surface $A$:
$$
P = \int_A I(p,\omega) dA'
$$

Each term of the integral will look something like this:
$$
\int_A \sigma_k^2(\omega) \delta(p - p_k) dA' = \begin{cases} 0 & p_k \notin A \\ \sigma_k^2(\omega) & p_k \in A \end{cases},
$$
so in other words, you are just adding up the power coming from each source that your receive aperture $A$ captures, and each source contributes power $\sigma^2_k(\omega)$.

Implicit in this formulation is that all the sources are incoherent, i.e., the waves do not have a fixed phase relationship and thus the average power of the superposition is the sum of the contributions from the individual sources.


Power Spectral Density of a Filter I need to calculate the output power spectral density of the following digital filter

My calculations are as follows:
$y\left(n\right)\:=\:x\left(n-1\right)+d\left[x\left(n-1\right)+x\left(n\right)\right]$
$Y\left(z\right)\:=\:X\left(z\right)z^{-1}+d\left[X\left(z\right)z^{-1}+X\left(z\right)\right]$
$Y\left(z\right)\:=\:X\left(z\right)\left[\left(1+d\right)z^{-1}+d\right]$
$\frac{Y\left(z\right)}{X\left(z\right)}=\:H\left(z\right)=\left(1+d\right)z^{-1}+d$
Did I get this equation right? would somebody help me please because this diagram is kind of confusing for me.

A digital filter doesn't have power spectral density, a signal has. I guess what you want is the transfer function of this filter, i.e., $H(z)=Y(z)/X(z)$.
Let's add two intermediate signal $u(n)$ and $v(n)$ in the diagram.

We have:
$$
y(n) = v(n-1) + du(n) \tag{1}
$$
$$
u(n)=v(n-1)+x(n)\tag{2}
$$
$$
v(n) = x(n) -du(n) \tag{3}
$$
Take Z-transform and we get
$$
Y(z) = V(z)z^{-1} + dU(z) \tag{4}
$$
$$
U(z) = V(z)z^{-1}+X(z)\tag{5}
$$
$$
V(z)=X(z)-dU(z)\tag{6}
$$
Furthermore, we can derive that
$$
V(z)=X(z)-d\Big[V(z)z^{-1} + X(z)\Big]
$$
$$
\frac{V(z)}{X(z)} = \frac{1-d}{1+dz^{-1}}
$$
$$
Y(z)=V(z)z^{-1} + d\Big[V(z)z^{-1}+X(z) \Big] =(1+d)V(z)z^{-1} + dX(z)
$$
Thus,
$$
H(z) = \frac{Y(z)}{X(z)} = (1+d) \frac{V(z)}{X(z)} z^{-1} + d = (1+d) \frac{(1-d)z^{-1}}{1+dz^{-1}} +d
$$