Differential Pulse Code Modulation (DPCM) In the following figure we see an encoder and a decoder of a Delta Modulation System, which is a simple DPCM System with a 1-bit quantizer:

We assume that signal $x(n)$ is known, that $\hat{x}(-1)=0$ and that there is not noise in the telecommunication channel.

I want to define the predicted signal $\tilde{x}(n)$, the difference signal $d(n)$, the quantized difference signal $\hat{d}(n)$, the encoded signal $c(n)$ and the reconstructed signal $\hat{x}(n)$

I have written the following equations:
$$d(n) = x(n) - \tilde{x}(n) \; \; \; \; (1) $$
$$\tilde{x}(n) = 0,8 \cdot \hat{x}(n-1) \; \; (2)$$
$$\hat{x}(n) = \tilde{x}(n) + \hat{d}(n) \; \; \;\;(3)$$
In addition, duo to the feedback, the error $e(n)$ for the reconstruction (between the initial signal on the transmitter and the reconstructed signal on the receiver) is equal to the error of the quantizer. So, I also get:
$$\hat{d}(n) = d(n) + e(n) \; \; \; \; (4) $$
$$\hat{x}(n) = x(n) + e(n) \; \; \; \; (5) $$

After solving the system of the equations, it seems to me that the number of the uknowns exceeds the number of the linear independent equations. If it so, I need one more equation since not all the equations above are linear independent.

So my question is what am I missing? Thanks in advance!

Eq. $(5)$ in your question is really just a consequence of Eqs $(1)$, $(3)$, and $(4)$, so you have four independent equations. But that's all you need because there's a simple one-to-one correspondence between $\hat{d}[n]$ and $c[n]$ (given by the encoder). Note that the quantization error $e[n]$ is usually modeled as a zero-mean white random process with a uniform PDF.

You generally cannot compute an analytical solution of all the sequences in the coder due to the quantizer, which is of course a non-linear device. However, you can compute the sequences sample by sample by simply simulating the DPCM system as shown in the diagram in your question.

Assuming that the prediction filter has an initial state of zero (i.e., $\tilde{x}[0]=0$), and that the quantization step $\Delta=1$ for simplicity, we can easily compute all sequences given the input sequence. With an input sequence of x=[2,3,5,7] we obtain

$\begin{array}{cccrcr}
n & x[n] & \tilde{x}[n] & d[n] & \hat{d}[n] & \hat{x}[n] \\
0 & 2 & 0.00 & 2.00 & 1 & 1.00 \\
1 & 3 & 0.80 & 2.20 & 1 & 1.80 \\
2 & 5 & 1.44 & 3.56 & 1 & 2.44 \\
3 & 7 & 1.95 & 5.05 & 1 & 2.95
\end{array}$


what does the output of ifft in ofdm look like? i have a problem in ofdm modulation .... there is a stream of bits that is converted from serial to parallel and is applied to an IFFT block. based on the definition of IFFT the output should be an complex valued number. how does it represent the real and imaginary parts in the digital domain( bits 0,1) ?

You are correct that the output of the IFFT will normally have both a real and imaginary component.  That is fine and normal.  When the signal is modulated from baseband to the carrier frequency the real part is typically modulated by a sine wave and the imaginary part by a cosine (or vice versa).  This allows both parts to coexist without interfering with each other, since sines and cosines have zero cross-correlation.

Regarding how the numbers are represented digitally, they are represented the exact same way that real numbers are represented.  The only difference is that there is an extra set of numbers for the imaginary part of the signal.


Sinusoid amplitude loss in power spectrum estimation using Pwelch method I am trying to estimate the amplitude or power spectrum of a sinusoidal signal using Pwelch method. The following is the matlab code I use.

rng default
fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + 0*randn(size(t));
[pxx,f] = pwelch(x,fs,0,fs,fs,'onesided'); plot(f,10*log10(pxx))

And it gives the following result

By default, a Hamming window of length $f_s$ is used in the pwelch function. Since the frequency resolution is 1 Hz, the frequency of the sinusoid lies exactly in one of the discrete frequency bins of the FFT and no scalloping loss is involved.
The question I have is that why the amplitude or power of the peak is -4.357dB instead of -3dB. I tried other windows with the same length. The rectangular window gives a peak of exactly -3 dB and the Hanning window, -4.78 dB. The power loss maybe related to the type of windows, I guess. But I cannot figure out why.

Can anyone explain this?

Thank you in advance.

I think your difficulty lies in having an overlap size of 0.

I'd recommend having an overlap that gives a Constant Overlap Add result when the windows are laid on top of each other.

For a rectangular window that is 0.
Hann it is 50% blocksize
Hamming seems to have many, but try 75%

Overlap values for different Windows

And don't forget that the PSD 'bin' represents all the content within that delta-f.  For example, if you had a signal at 100 Hz and 100.25 Hz, the 1 Hz delta-f wouldn't resolve them to 2 spikes .

RMS from a bin is sqrt(Area of bin).


When Two Sampled Sinusoidal Are Orthogonal? If two analog sinusoidal are orthogonal with duration T, their minimum frequency difference should be 1/2T in case of no phase offset between them.

And if there is phase offset, the difference is 1/T.

And in N sampled domain the orthogonality depends on sampling frequency.

The condition is $f_{s}=Nf_{k}/k$ obtained from $f_{k}=kf_{s}/N$.

($f_{s}$:sampling frequency, $f_{k}$: frequency difference, $k$: integer bigger than 0, $N$:# of sample)

Ex.)$N=2048, k=1, f_{k}=15kHz$ then $f_{s}=30.72MHz$

I like to know how $f_{k}=kf_{s}/N$ is derived

Thanks

If you have two sinusoidal length $N$ sequences

$$x_1[n]=\sin\left(n\omega_1+\phi_1\right),\qquad n\in[0,N-1]\\
x_2[n]=\sin\left(n\omega_2+\phi_2\right),\qquad n\in[0,N-1]\tag{1}$$

then for orthogonality we require

$$\sum_{n=0}^{N-1}x_1[n]x_2[n]=0\tag{2}$$

For the given sequences $(1)$ we obtain

$$\begin{align}\sum_{n=0}^{N-1}x_1[n]x_2[n]=&\sum_{n=0}^{N-1}\sin\left(n\omega_1+\phi_1\right)\sin\left(n\omega_2+\phi_2\right)\\=&\frac12\sum_{n=0}^{N-1}\cos\left[n(\omega_1-\omega_2)+\phi_1-\phi_2\right]-\\&\frac12\sum_{n=0}^{N-1}\cos\left[n(\omega_1+\omega_2)+\phi_1+\phi_2\right]\tag{3}\end{align}$$

where I've used

$$\sin(\alpha)\sin(\beta)=\frac12\big[\cos(\alpha-\beta)-\cos(\alpha+\beta)\big]\tag{4}$$

The two terms on the right-hand-side of $(3)$ have the form

$$\sum_{n=0}^{N-1}\cos(n\omega+\phi)=\cos\left(\frac{N-1}{2}\omega+\phi\right)\frac{\sin(\omega N/2)}{\sin(\omega/2)}\tag{5}$$

For large $\omega$ (i.e., $\omega$ close to $\pi$), the magnitude of the term $(5)$ becomes small. Assuming that $\omega_1+\omega_2$ is large (but not greater than Nyquist), people tend to neglect the last term on the right-hand side of $(3)$.

From $(5)$, the first term on the right-hand side of $(3)$ is zero if

$$\frac{\sin((\omega_1-\omega_2) N/2)}{\sin((\omega_1-\omega_2)/2)}=0\tag{6}$$

is satisfied, i.e., for

$$\omega_1-\omega_2=\Delta\omega=\frac{2\pi k}{N},\quad k\in\mathbb{Z},\quad k\neq 0\tag{7}$$

With $\Delta \omega=2\pi\Delta f/f_s$ we obtain for the frequency difference

$$\Delta f = \frac{kf_s}{N},\quad k\in\mathbb{Z},\quad k\neq 0\tag{8}$$

Note that choosing the frequency difference according to $(8)$, the two sinusoids will be approximately orthogonal.
If we don't neglect the second term on the right-hand side of $(3)$ we can make them exactly orthogonal. However, in that case we cannot anymore freely choose one of the two frequencies, but the frequencies need to satisfy

$$f_1=\frac{(k_2-k_1)f_s}{2N},\quad f_2=\frac{(k_1+k_2)f_s}{2N},\\\quad k_1,k_2\in\mathbb{Z},\quad k_2>k_1\tag{9}$$


Algorithms for detection of fingertips or nails Can you give me some advices on what could be the best algorithm(s) to use for detection of fingertips/nails in image. First thing that crossed my mind was Viola - Jones. After rethinking I concluded that maybe it would be possible to use just Hugh transformation after applying edge detection. But I'd like some more advices. Also, as this is going to be a student's project with a purpose of learning, I'am not allowed to use OpenCV or similar frameworks. Below is typical image that will be processed. (note it's not top-down view). There is no need for thumb detection.

typical image http://www.deviantpics.com/images/BwgPX.jpg

I would consider using neural network or SVM to fit the model. The difficulty with this approach is that you must collect a lot of data - both positive and negative examples. But you can generate a lot of artificial data(by scaling or rotating images you have already collected). After you collect data you can use "moving window" of few sizes to detect nails on images not present in training set.
I don't know how much time you want to spend doing the project. But implementing efficient learning algorithm yourself is difficult. However there is already library for SVM which I would use for this purpose.


How to understand Kalman gain intuitively? The Kalman filter algorithm works as follows

  Initialize $ \hat{\textbf{x}}_{0|0}$ and $\textbf{P}_{0|0}$.

  At each iteration $k=1,\dots,n$

  Predict

  Predicted (a priori) state estimate    $$ \hat{\textbf{x}}_{k|k-1} =
 \textbf{F}_{k}\hat{\textbf{x}}_{k-1|k-1} + \textbf{B}_{k}
 \textbf{u}_{k} $$
    Predicted (a priori) estimate covariance $$
 \textbf{P}_{k|k-1} = \textbf{F}_{k} \textbf{P}_{k-1|k-1}
 \textbf{F}_{k}^{\text{T}} + \textbf{Q}_{k}$$
  Update

  Innovation or measurement residual     $$ \tilde{\textbf{y}}_k =
 \textbf{z}_k - \textbf{H}_k\hat{\textbf{x}}_{k|k-1}$$ Innovation (or
  residual) covariance   $$\textbf{S}_k = \textbf{H}_k
 \textbf{P}_{k|k-1} \textbf{H}_k^\text{T} + \textbf{R}_k$$ Optimal
  Kalman gain      $$\textbf{K}_k =
 \textbf{P}_{k|k-1}\textbf{H}_k^\text{T}\textbf{S}_k^{-1}$$
    Updated (a posteriori) state estimate      $$\hat{\textbf{x}}_{k|k} =
 \hat{\textbf{x}}_{k|k-1} + \textbf{K}_k\tilde{\textbf{y}}_k$$
    Updated (a posteriori) estimate covariance     $$\textbf{P}_{k|k} = (I -
 \textbf{K}_k \textbf{H}_k) \textbf{P}_{k|k-1}$$

The Kalman gain $K_k$ represents the relative importance of the error $\tilde{\textbf{y}}_k$ with respect to the prior estimate $\hat{\textbf{x}}_{k|k-1}$.

I wonder how to understand the formula for the Kalman gain $K_k$ intuitively? Consider the case when the states and outputs being scalar, why is the gain bigger, when

$\textbf{P}_{k|k-1}$ is bigger
$\textbf{H}_k$ is bigger
$\textbf{S}_k$ is smaller?

Thanks and regards!

I found a good way of thinking intuitively of Kalman Gain $K$. If you write $K$ this way

$\displaystyle \quad\ \bf{K_k} = \bf{P_k^-\, H_k^{\rm T} (H_k P_k^-\, H_k^{\rm T} + R_k)^{-1}}
  = \bf{\frac {P_k^-\, H_k^{\rm T}}{H_k P_k^-\, H_k^{\rm T} + R_k}}$

you will realize that the relative magnitudes of matrices ($R_k$) and ($P_k$) control a relation between the filter's use of predicted state estimate ($x_{k}⁻$) and measurement ($ỹ_k$).

$\displaystyle \quad\
  \lim\limits_{\bf{R_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
    = \bf{H_k^{-1}}$

$\displaystyle \quad\
  \lim\limits_{\bf{P_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
   = \bf 0$

Substituting the first limit into the measurement update equation

$\displaystyle \quad\
  \bf{\hat x_k} = \bf{x_k^-} + \bf{K_k}(\bf{\tilde y_k}-\bf{H_k}\bf{x_k^-})$

suggests that when the magnitude of $R$ is small, meaning that the measurements are accurate, the state estimate depends mostly on the measurements.

When the state is known accurately, then $H P^⁻ H^T$ is small compared to $R$, and the filter mostly ignores the measurements relying instead on the prediction derived from the previous state ($x_k⁻$).


What is the right way to downsample using Fourier method? I want to know what is the right way to downsample a sampled signal using Fourier transform as the implementation in scipy.signal.resample confuses me.

Reading through the code it first converts the signal to frequency domain, then discards the middle half of the frequencies (i.e. second and third quarters), then inverse-transform back to time domain. Therefore it is trying to keep the lowest and highest frequencies of the signal.

Due to Nyquist-Shannon we know that we cannot reconstruct a frequency of $f$ without at least $2f+1$ samples of the signal. Then when downsampling, shouldn't we instead discard the second half of the frequencies?

I already read some related answers but no one seems to address this specific question:

Scipy resample, "fourier method" explanation

Python's $\tt resample$ vs $\tt resample\_poly$ vs $\tt decimate$

Note that the result of a DFT is inherently periodic, so the second half of the resulting output vector does not correspond to the highest frequencies, but to the negative frequencies. The highest frequency (Nyquist) is in the middle of the vector, after that you get the most negative frequency and the last element is the smallest negative frequency just below DC. That's why when taking out the middle part of the vector you actually remove the highest frequencies.


Suggestions on how to use FIR model in sentence retrieval task Estimating an impulse response from input-output measurements is called system identification. When the impulse response of the
linear system $h(n)$, is of finite duration (M samples) or can be
effectively modeled with a finite impulse response system.  System identification, that is, the modeling and identification of a system from knowledge of its input $x(n)$ and output signals $y(n)$ is known as non-blind whereas if the knowledge of the input and the  system both are unknown by the receiver end then it is termed as blind system identification. In wireless communications and in signal processing research papers I have seen that the channel impulse response is modeled as a moving average. In communications, there is an additive noise known as the measurement noise $e(n)$. Mathematically expressing,
$$y(n) = x(n)*h(n) + e(n)$$

In signal processing and estimation, the objective is to estimate the channel impulse response and recover the input from the noisy observations.

How can I use the above model where the input signal is a sentence. So, a sentence would be represented by assigning each unique word by a symbol. Let's say, I have the sentence "the quick brown fox jumped over the lazy dog" and I have the vocalbulary and mapping as

$the \mapsto 1$

$quick \mapsto 2$

$brown \mapsto 3$

$fox \mapsto 4$

$jumped \mapsto 5$

$over \mapsto 6$

$lazy \mapsto 7$

$dog \mapsto 8$

So, the sentence would be mapped to $[1,2,3,4,5,6,1,7,8]$.

This is the input signal or the input data.

I cannot understand what the role of the channel impulse response would be if

I want to estimate after decoding the input speech signal from noisy observations.
What would be the noise term and
would there be any "channel impulse response" ? What would be the best way to use the FIR model for this kind of data?

How can I use the above model where the input signal is a sentence.

makes no sense; your whole system is not linear (which is one of the requirements for your whole convolution approach); for example, if you you play back the words "fox" and "jumped" at the same time, you don't hear "lazy".

Not linear means your system can't be represented by an impulse response.

Your signal model is overly simplistic, too. You might map words to numbers as soon as you know the words, but that's only happening after speech recognition (and might be helpful in understanding logically what has been said), but not before.

So, your whole approach shows you haven't dealt with a lot of speech recognition before – which is OK, but it means that you can't just come up with some approach and assume it makes sense from the start.

I'd recommend reading a book on speech recognition. There's certainly a big wealth of literature out there.


Why am I not getting a flat phase for a Gaussian pulse when doing a Fourier transform in Python? I have been trying to obtain a spectrum and a spectral phase of a Gaussian pulse using the Fast Fourier Transform provided with numpy library in Python. Here are the results:

It is known that the spectral phase of a Fourier-limited Gaussian pulse should be flat (i.e. equal to some constant across the whole spectrum). However, after having Fourier-transformed the pulse, I am getting a strange spectral phase as a result. I find this phase using the function angle provided with numpy which should calculate the phase of a complex number. Is there any idea what I a missing here?
Here is the code that I use.

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import speed_of_light
from scipy.fftpack import fft, fftshift, ifft

# construct Electric field of a Gaussian pulse in time domain
def buildPulse(delays, wavelength = 700, fwhm_fs = 20, time_position = 0, phi0 = np.pi/2):
    frequency = speed_of_light/wavelength*10**9 # frequency in Hz
    period = 1/frequency*10**15 # period in fs
    sigma = fwhm_fs*(2*np.sqrt(2*np.log(2)))**(-1)
    envelope = np.exp(-(delays-time_position)**2/2/sigma**2)
    Ef = envelope*np.sin(2*np.pi*frequency/10**15*delays-phi0) # electric field
    return Ef, envelope, frequency, period, fwhm_fs

# fast Fourier Transform
def FFT(Ef,t0,tf,t_step,delays,fwhm_fs):
    # get frequency axis in Hz and wavelength axis in nm
    frequencies = np.linspace(10**(-25), (2*((tf-t0)/((tf-t0)/t_step+1)))**(-1), delays.size//2)*10**15
    clbr_const = (frequencies[3]-frequencies[2])*10**15 # freq axis calibration const
    print("clbr_const = ",clbr_const)
    wavelengths = speed_of_light*frequencies**(-1)*10**(9) # in nm
    wavelengths = wavelengths[::-1] # reversing wvelength axis to make it increasing
    # fast Fourier transform of Ef (electric field in time domain)
    spectrum = fft(Ef)
    # half of the spectrum-vector
    spectrum_half = spectrum[0:len(spectrum)//2]
    # normalization factor
    A = 2*(tf-t0)/spectrum.size/fwhm_fs
    return spectrum,spectrum_half,frequencies,wavelengths,A

# construct time-axis (delays in fs)
t0 = -200 # first time-point in fs
tf = 200 # last time-point in fs
t_step = 0.01  # time-step in fs
delays = np.linspace(t0, tf, int((tf-t0)/t_step)+1)

# construct the electric field of a Gaussian pulse
Ef, envelope, frequency, period, fwhm_fs = buildPulse(\
    delays, wavelength = 700, fwhm_fs = 20, time_position = 0, phi0 = np.pi/4)
print("frequency = ",frequency," Hz\nperiod = ",period,' fs')

# fast Fourier transform
spectrum,spectrum_half,frequencies,wavelengths,A = FFT(Ef,t0,tf,t_step,delays,fwhm_fs)
print("len(spectrum) = ",len(spectrum))

# amplitude spectrum
amp_s = A*np.abs(spectrum_half)
# real and imaginary parts of the spectrum
real_s = A*np.real(spectrum_half)
imag_s = np.imag(spectrum_half)
arg_s = np.angle(spectrum_half)

# plotting the results
plt.figure()
plt.plot(delays,Ef,delays,envelope)
plt.xlim([-100, 100])
plt.legend(["Electric field (Ef)","Envelope"],loc="upper left")
axs[0,0].set_xlabel('time (fs)')
plt.figure()
plt.plot(frequencies/10**12,amp_s)
plt.xlim([0, 1100])
plt.legend(["abs(fft(Ef))"],loc="upper left")
axs[1,2].set_xlabel('frequency (THz)')
plt.figure()
plt.plot(frequencies/10**12,arg_s)
plt.xlim([0, 1100])
plt.legend(["angle(fft(Ef))"])
plt.xlabel("frequency, THz")

The group delay will be flat assuming the filter is linear phase (which means the phase is linear vs frequency). Group delay is the derivative of phase with respect to frequency. Either unwrap your phase and detrend it or use the group delay function that is part of scipy.signal to evaluate directly. Any FIR filter (not just the Gaussian) will be linear phase if (and only if) the filter coefficients are symmetric or anti-symmetric.


Playing back an audio file at an arbitrary rate with upsampling and downsampling Background: My overall goal is to create an audio file player which can playback an audio file at an arbitrary rate. The rate can change continuously as the file plays. This means upsampling to slowdown the playback when the rate < 1, and downsampling to speedup the playback when the rate > 1. (I'm not concerned with maintaining pitch.)

Let's focus on the slowdown (upsample) part to make the question simpler.

I read that to upsample by a rational fraction L/M, I should first

Upsample by L
Apply an interpolation (lowpass) filter
Downsample by M

However, this doesn't work well for an online, live context, because L or M can be very large depending on the desired playback rate, requiring large buffers for upsampling/downsampling. I want to process the audio file in small chunks as it plays, so my buffers must be kept to a reasonable size.

I decided to try the following algorithm instead, for upsampling by an arbitrary floating point rate:

// pseudo-code upsample algorithm
sloppy_upsample(int source[], int destination[], float rate, int frames):
  for each frame i:
    destination[i] = source[round(i * rate)];

This is a sloppy "zero order hold" algorithm that probably doesn't get the sample rate exactly correct, but it actually sounds surprisingly good for my purposes. I also like this algorithm for its simplicity. But, of course, there's audible aliasing.

I know that I can apply a lowpass filter before upsampling to prevent aliasing, but I'm not sure what frequency to set the cutoff at. In similar situations such as this one, the cutoff could be set to less than L * 0.5.

Is it possible to select a good cutoff for my sloppy_upsample function above that prevents aliasing?
Is there another approach for upsampling/downsampling that works well for live playback?

Some solutions that are ordered by complexity, simplest first.

The simplest solution is the one you have already implemented- pick the nearest value.
The next simplest is linear interpolation.  You use the two nearest values and use a weighted average based on how far the point is from each.  For instance, if the point you want is .3 samples past $x_1$ and .7 samples before $x_2$, the value would be $y = x_1*.7 + x_2*.3$.  This is still pretty simple and should have much better sound performance than your current solution.
The next simplest is quadratic interpolation, or, more generally, polynomial interpolation.  Polynomial interpolation is also sometimes referred to as spline interpolation.
The most complicated solution is Lanczos resampling.  This is where you use a windowed sinc filter to low-pass filter the data.  You pick the particular set of taps to use from the filter based on the particular delay (i.e. at what point you are between the samples).


Calculate Max Number Of Harmonics for given Sample Rate? How do I calculate the maximum number of harmonics a complex wave can have based on the sample rate?

Say I Have a 20hz fundamental and a sampling rate of 44100 hz. My instinct is to take 44100/2 = 22050(nyquist) and the take the fundamental and multiply by 2 while the frequency is less than nyquist and the count how many times I had to multiply.

I would like to know if i'm going in the right direction and if i'm not... what is the correct way to do this?

EDIT:
I would like to add that I am creating a table based oscillator with multiple tables with different harmonic contents.I am trying to figure out how to calculate the harmonics each table should have if I have an arbitrary number of tables. So table 1 would have say 500 harms so a max freq of 44.1hz and table n would have x harms and have a max freq of y.

No, you don't just multiply by 2, you multiply by every integer.  You can have a 2nd harmonic, a 3rd harmonic, 4th, etc.  (Only multiplying by 2s would tell you the number of octaves, not the number of harmonics.)

So take Nyquist frequency, and divide by the fundamental frequency

  22050 Hz / 20 Hz = 1102.5

Then round down to the nearest integer.  So you can fit 1102 harmonics.  The 1102th harmonic is 1102*20 = 22040 Hz, while the 1103th harmonic is 1103*20 = 22060 Hz, which is too high.

(Note that some consider the fundamental as the 0th harmonic while others consider it the first harmonic.  I am counting it as the first, 1*f.)


Norm MPEG-1 Layer III (Mp3) Filter banks distortion So, I'm working with the filter banks used in Mp3, in Matlab. I am obtaining quite a huge distortion in the output and I don't really know whether I'm making something wrong.

I've reviewed the literature and I'd say I'm making it properly, so I don't really understand what I'm making wrong.

So I'm basing the work on a filter with 512 coefficients that can be found on page 9-13, here: https://www.mp3-tech.org/programmer/docs/jacaba_main.pdf (I think it comes from a standard so no big deal there, I guess).

I want to make a filter bank with 32 branches. I want to, based on that Low Pass filter prototype, calculate the 32 branches that will cover the entire bandwidth, so I make:

(note "h" is the variable with the filter coefficients)

branches=32;
hh=zeros(branches-1,length(h));
n=0:511;
for k=0:1:branches-1
   hh(k+1,:)=h.*cos(((2*k+1)*pi*(n-16))/(2*32));
end

That one for the generation of the analysis filters coefficients and the following one for the synthesis filters:

branches=32;
hh=zeros(branches-1,length(h));
n=0:511;
for k=0:1:branches-1
   hh(k+1,:)=h.*cos(((2*k+1)*pi*(n+16))/(2*32));
end

As per the way of applying the filters to the input signal and to reconstruct the output signal:

yy=zeros(branches,floor(length(wav_in)/32)+1)
for k=1:1:branches
  yy(k,:)=decimate(filter(hh(k,:),1,wav_in),branches);
end

This will provide a matrix with all the outputs of the shifted analysis filters (yy).

Now, in order to reconstruct the signal, with the synthesis filters:
(Note it's a different function so we will not have problems with the variables names)

yy = zeros(1,length(xd(1,:))*32)
for k=1:1:branches
  yy=yy+(filter(hh(k,:),1,upsample(xd(k,:),branches)));
end

The filter responses look good (it seems like that to me, at least):

If I use as an input a 1 kHz tone, I have:

The output is "After". As you can see the signal is very distorted. In general, I am following the approach in "How sound is processed in MP3 player?" by T. Dutoit and N. Moreau.

I think it is not working and honestly, I don't have a single clue on what I'm doing wrong. (Maybe how I'm applying the filter banks?).

I have tested your code.

There are problems with the decimation and upsample stages. When you run the analysis and synthesis filters at the undecimated rate, then you get the perfect reconstruction (plus a filter shift)

First of all in case of applications which are sensitive to group delay, you should better use the conv / filter function and then compansate for the group delay.

Also, the decimate function applies a lowpass anti-aliasing filter, which kills the signal content at the $k$-th channel considering the fact that its output is a bandpass signal rather than lowpass. Hence you should avoid the anti-aliased decimation but instead employ a bandpass decimation strategy.
A direct subsampling by 32 would provide a convincing result.

Same goes with interpolation part. After you expand your signal, you should select the proper band of image spectrum by a suitable bandpass filter, rather than a lowpass filter from a conventional interpolator. A direct expansion may convince as its output will already be bandpass filtered later by the polyphase filterbank.

That being said, your main problem is the document that you based your code. That's far too out of signal processing context. I suggest you read the MPEG audio standard book of BOSI. But it's up to you of course.

The following is a working corrected MATLAB / OCTAVEcode:

% S0 - Load the prototype lowpass filter impulse response:
% --------------------------------------------------------
load h2.mat;         % h[n] is the prototype lowpass filter... of length 512
L = length(h);

figure,subplot(2,1,1)
stem([0:L-1],h);title('The Prototype Lowpass Filter');
subplot(2,1,2)
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(h,4*L)))));
grid on;

% S1 - Create the 32 band filter bank by cosine modulation
% --------------------------------------------------------
numbands = 32;                  % number of bands
n=0:L-1;

hha=zeros(numbands,L);         % band of filters
for k=0:1:numbands-1
   hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end

figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
title('ANALYSIS FILTER BANK');
hold on
for k=1:numbands
    plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end

% S2 - Design the synthesis filterbank:
% -------------------------------------
numbands = 32;                     % number of banks
n=0:L-1;

hhs = zeros(numbands,L);           % band of filters
for k=0:1:numbands-1
   hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end

% S3 - Generate a test input signal:
% ----------------------------------
N = 2*1024;      % NOTE: According to the MPEG specification a determined LENGTH must be applied.
wav_in = sin(0.05*pi*[0:N-1]);        % pure sine tone

figure,subplot(2,1,1)
plot(wav_in);title('input signal')
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(wav_in,4*N)))));

% S4 - Apply a wav file to it,  get every output of the every sub-band:
% --------------------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands));   % decimated outputs..
yyu = zeros( numbands , N);

for k=1:1:numbands
    temp = conv( wav_in,hha(k,:));
    yyu(k,:) = temp(L/2+1:L/2+N);
    yyd(k,:) = yyu(k,1:numbands:end);
end

% S6 - Apply synthesis on the decimated signal
% --------------------------------------------
ys = zeros(1, N);
yi = zeros(1,N);

for k=1:1:numbands
    yi(1:numbands:end) = yyd(k,:);
    temp = conv(yi,hhs(k,:));
    ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;

% Display the result:
figure,subplot(2,1,1)
plot(ys);title('Synthesized Back');
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(ys,4*N)))));

With the outputs being:

1- The prototype filter:

2- 32 channel filter bank:

3- A pure tone input:

4- An reconstructed back output:


Determine the Unit Impulse Response of ARMA Filter for 8 First Indices I have this problem, I've checked for almost 3 hours in books and internet for the unit impulse response for n = 0 : 7 but nothing.

I've benn working on how to get h[t] but I end up on a dead end at:

[

From there, I don't know how to solve this problem, Someone could help me pls?

As you wrote correctly the function is given by:

$$ y \left[ n \right] = \frac{1}{2} y \left[ n - 1 \right] + \frac{1}{8} y \left[ n - 2 \right] + \frac{1}{2} x \left[ n - 2 \right] $$

Assuming initial conditions of 0 we can easily see that for unit impulse (Namely all zeros but $ x \left[ 0 \right] = 1 $:

$$ \begin{align*}
y \left[ 0 \right] & = 0 \\
y \left[ 1 \right] & = 0 \\
y \left[ 2 \right] & = \frac{1}{2} \\
y \left[ 3 \right] & = 0.5 * y \left[ 2 \right] = \frac{1}{4} \\
y \left[ 4 \right] & = 0.5 * y \left[ 3 \right] + 0.125 * y \left[ 2 \right] = \frac{3}{16} \\
y \left[ 5 \right] & = 0.5 * y \left[ 4 \right] + 0.125 * y \left[ 3 \right] =
\frac{1}{8} \\
y \left[ 6 \right] & = 0.5 * y \left[ 5 \right] + 0.125 * y \left[ 4 \right] = \frac{11}{128} \\
y \left[ 7 \right] & = 0.5 * y \left[ 6 \right] + 0.125 * y \left[ 5 \right] = \frac{15}{256} \\
\end{align*}
$$

Since it has IIR Component it will go on infinitely.

MATLAB Simulation:

The full code can be found in my StackExchange Signal Processing Q43918 GitHub Repository.


Aperture Jitter in Digital to Analog Converters I have noticed in the literature that there is an internally induced timing jitter inside Analog to digital converters called (aperture jitter). It is caused by the sample-and-hold circuit of the ADC.
So what about the digital to analog converter? does it contain a sample-and-hold circuit? does it induce aperture jitter?
thanks

Yes absolutely there is also the equivalent of aperture jitter results in a DAC. This is due to the effects of jitter in the sampling clock itself, and variation in the electronics of the time duration for the translation from the different digital levels to analog levels in the output.

In either case the jitter is not necessarily just due to a S/H circuit but the combined effects of that interface as well as the clock itself, and power supply noise can also translate to aperture jitter given the resulant changes in threshold levels.


Unexpected frequency components found in Periodogram I am trying to make sense of Digital Signal Processing with the help of R language. So, the issue came up when I tried to create a periodogram for a simple sine wave. The regular periodogram shows 2 frequency components as expected (expected frequency & aliased frequency). However, when I looked at the contents of the Power values for this periodogram, I realised that these apart from these 2 frequency components there are other components too, which have relatively very small values(order of e-30). Since, compared to the major 2 frequencies, all the other frequency components appear to be ~ 0, therefore the simple periodogram plot doesn't show them. However, if I plot the periodogram, after applying log function, then these frequency components(the insignificant ones) become visible.

Code :

t = 1:100            #Time instances
freq_t = 0:99/100    #Fundamental frequencies : t-1/n ,where n is the no. of time instances
y = sin(2*pi*0.06*t) #Sine wave with the time instances 't'
fft_y = fft(y)       #Performing FFT on the 'y' signal
spec_y = (abs(fft_y)^2)/100 #To obtain the modulus values of the complex values present in 'fft_y'
log_spec_y = log10(spec_y) #applying log at base 10 to the 'spec_y'

Note : The major frequency components (0.06 & 0.94) are marked in red

So, my issue is :

Why do these unwanted frequency components even exist? Shouldn't
they have their corresponding power values as 0?
How will these minor frequency components affect any sort of analysis that can be performed in the frequency domain?

Speculation for the first issue:
Now, for existence of those minor frequency components, the only reason that I can come up with is that the sine wave I made, is basically a sampled signal, hence discrete in nature. In order to approximate this sampled signal, apart from the major frequency components a lot of other minor frequency components should have been used, and that's the reason they show up like this. For instance, consider the Fourier series for square waves. Here, the fitted Fourier series looks like a summation of multiple sine/cosine signals with different frequencies. And the highest frequency of all those sine/cosine signal determines the frequency of the square wave.

As far as the second issue is concerned, I am pretty much clueless about that.
Hence can you help me resolve these 2 issues?

First of all, if you use a logarithmic graph, it is customary to use dB. In order to convert a FFT in decibels, we use 20*log10, not log10.

Second of all, the "unwanted" frequencies are 10^-30 times smaller than the wanted frequencies, i.e. they are insignificant. . The unwanted frequencies are most likely due to the limited precision of floating-point arithmetic and not due to the sampling process.

https://en.wikipedia.org/wiki/Floating-point_arithmetic

That should answer your second point too...


How to prove that frequency mixer is a nonlinear system? I know that frequency mixer $y(t)=x(t)cos(\omega_ct)$ is a nonlinear system since its output has frequency components that are not present in the input. To prove its nonlinearity, it shuldd be shown that the system either violates scaling or additivity. However, I could not do that. Can someone please help me?

If the system is nonlinear then if $y_1(t)$ is the response to the signal $x_1(t)$, and  $y_2(t)$ is the output given input signal $x_2(t)$ then the response to the signal

$$x(t)=a_1x_1(t)+a_2x_2(t)\tag{1}$$

with arbitrary constants $a_1$ and $a_2$ will generally not be equal to

$$y(t)=a_1y_1(t)+a_2y_2(t)\tag{2}$$

However, for the given system an input signal (1) always results in an output signal (2). The conclusion is that the system is linear.

The confusion arises because you think that linear systems cannot add new frequencies to the input signal. This is generally not true. Only linear and time-invariant (LTI) systems have that property. Linear but time-varying systems can produce new frequencies in the output signal, and a modulator is a classic example of a linear time-varying system.


Finding linear chirp parameters from a low quality spectrogram, possibly with image processing techniques? tl;dr: What is a fast way to get the equation of this line?

I have to extract the "line of best fit" from this data (shown below), and many others like it (so I need a robust method). It is a a spectrogram of a linear chirp. The x axis is frequency and the y-axis is time. The linear chirp came from a truck outputting 4-140Hz linearly over exactly 8 seconds, but attenuation in the ground and some other physical phenomena mean that we are seeing a signal for around 10 seconds.

I looked into using the Hough transform to find the gradient of the line, and then cross-correlating the data with a line of this gradient. But the Hough transform is too slow to perform en mass, which is what I need to do.

I thought about finding the maxima for frequency ranges (5-10, 10-15, etc.) and plotting a line of best fit through these, but in this particular example there are only two or three areas of strong intensity.

The other spectrograms look like this but the high intensity regions will move around to different points on the line, and some may have more noise.

Whats a good way of doing this algorithmically?

My current idea is:

Look for highest intesity regions.
Look at many pixels in their neighbourhood and estimate the "center of mass".
Draw line between them and find equation of line.

Is there a better way?

I think you can use standart method of fitting data to Straight Line. See for example explanation in 15.2 of "Numerical Recipes in C" http://apps.nrbook.com/c/index.html (page 661-662)

To prepare arrays of $x_i$, $y_i$, $\sigma_i$, you must for each point on your image (frequency-time):
$x_i$ := frequency
$y_i$ := time
$\sigma_i$ := 1/intensity (or 1/intensity$^2$)
size of arrays will be N*M (N, M - x and y size of spectrogram in points)
Calculate a and  b (y = a + b*x) by using formula (15.2.6) from http://apps.nrbook.com/c/index.html
If it will be some problem with outliers, use some robust fitting method. See for example 15.7 of  http://apps.nrbook.com/c/index.html. But I think there will not such problem with your data.

This method is faster then Hough transform

EDIT
It is very important to choose good expression for $\sigma_i$. May be better results will be with $\sigma_i$ := $1/\sqrt{intensity}$. May be will be better to use $\sigma_i$ := 1/(intensity-background) for (intensity > background) and exclude all points with (intensity < background). If there are some limitations on parameres of straight line (a,b), you can exclude some x-y points from $x_i$-$y_i$ arrays. But this modifications depend on actual data. It will be very interesting to see the your results.


Difference between demodulation and band-stop Filter I was wondering what is the difference between demodulating a signal that is carried with a frequency $F_{\rm carrier}$ and applying a band-stop on a band around $F_{\rm carrier}$ ?

Filtering is a linear operation on the signal, and as such can alter amplitude an phase of signal components, but not their frequencies.

Demodulation is quite a general concept, but normally includes non-linear operations such as shifting the signal in the frequency domain from the carrier to baseband (or an intermediate frequency) in order to extract the unmodulated signal.

So you are trying to compare very different operations.


Are all of the edge descriptors necessary to differentiate edges? My question is about edge descriptors produced from an edge detection method applied on an object image. As we know, an edge has four descriptors edge normal, edge direction, edge position, and edge strength.  These descriptors are used for differentiating between the detected edges.

Are all of this descriptors are necessary? When can I dismiss any of these descriptors?

I am thinking about assumptions based on the properties of the input image, e.g., if the input image has no noise, can I neglect the edge strength and assume that all resulting edge points are true edges?

You can calculate the edge direction from the edge normal by atan2(y,x)link.
The other descriptors are independent of each other and thus only the problem itself can tell if they are necessary or not.


Is Savitzky-Golay filter a special case of an FIR lowpass filter? The Savitzky-Golay filter is explained here. It looks to me that any Savitzky-Golay filter can equally be considered an FIR filter. Am I right about that?

Yes. The Savitzky-Golay filter is a weighted sum of a finite number of data samples. So it is indeed an FIR filter (albeit, not causal, because it uses samples from the past and the future). And it is indeed a low-pass filter (smoother) as seen from the frequency response of the filter on the same page.


Using MFCC to an ANN Speech Recognition System I'm developing an Artificial Neural Network based Speech Recognition System using MFCCs.

Suppose I have 260 input nodes in the ANN, and this number of nodes corresponds to the number of MFCCs that I will use. During feature extraction the number of total coefficients vary with respect to the duration of the sound file. This poses a problem if the ANN was trained just for 260 coefficients.

So most likely the system will fail if a different sound duration which yields lesser or greater number of coefficients is used to test the Neural Network. My question is how do I go about this problem? I have seen several papers in the net talking about Speech Recognition using ANN but I haven't seen something concerning this problem

You should process the speech in blocks of a few milliseconds at a time. Each block will be the same length. Then you can apply your ANN to those features.


“The Fourier transform cannot measure two phases at the same frequency.” Why not? I have read that the Fourier transform cannot distinguish components with the same frequency but different phase. For example, in Mathoverflow, or xrayphysics, where I got the title of my question from: "The Fourier transform cannot measure two phases at the same frequency."

Why is this true mathematically?

It's because the simultaneous presence of two sinusoidal signals with the same frequency and different phases is actualy equivalent to  a single sinusoidal at the same frequency, but, with a new phase and amplitude as follows:
Let the two sinusodial components be summed like this :
$$ x(t) = a \cos(\omega_0 t + \phi) + b \cos(\omega_0 t + \theta) $$
Then, by trigonometric manipulations it can be shown that :
$$ x(t) = A \cos(\omega_0 t + \Phi) $$
where
$$A = \sqrt{ a^2 + b^2 + 2 a b \cos(\theta-\phi) } $$ and
$$ \Phi = \tan^{-1}\left(\frac{ a \sin(\phi) + b\sin(\theta) } { a \cos(\phi) + b\cos(\theta) } \right) $$
hence you actually have a single sinusoidal (with a new phase and amplitude), and therefore nothing to distinguish indeed...


Compare two signals, find times when they are affected by the same processes I have two (noisy) signals $S_1(t)$ and $S_2(t)$. They contain information on three distinct events occurring in a physical system. Event $A$ affects only $S_1(t)$, event $B$ affects both $S_1(t)$ and $S_2(t)$, and event $C$ affects only $S_2(t)$. I want to identify the sequence of events in the system by analyzing the signals.

Is there some generic way to identify if two strings of data contain correlated events and find the timing of said events?

More specifically, my signals take the form
\begin{equation}
\textrm{d}S_i (t + \textrm{d}t) = S_i^{\textrm{true}}(t)\textrm{d}t + \textrm{d}W(t),
\end{equation}
where $S_i^{\textrm{true}}$ is the noiseless value of the signal calculated from the state of the system at time $t$, $\textrm{d}W(t)$ is a Wiener element (uncorrelated Gaussian noise), and the statement is to be read as an update formula. The l.h.s. is what is available to me. In case of an event the true values change sign in the following manner:
\begin{align*}
A &\Rightarrow\ S_1^{\textrm{true}} \to -S_1^{\textrm{true}}, \\
B &\Rightarrow\ S_1^{\textrm{true}} \to -S_1^{\textrm{true}} \textrm{ and }  S_2^{\textrm{true}} \to -S_2^{\textrm{true}},\\
C &\Rightarrow\ S_2^{\textrm{true}} \to -S_2^{\textrm{true}}.
\end{align*}
The way I am trying to do things right now: I am using a running average to get rid off some of the noise, then I find the second derivatives of the smoothed signals. If the true value changes, there will be a peak in the second derivative. I then determine whether there is a peak in both signals or just one at the given time. However, this does not work very well, as the derivatives pick up the noise, too.

Here is a plot of my raw and running window averaged signals:

I have no experience with signal processing. Is there a canonical procedure? If the question is too broad, where can I find recipes to deal with issues as such?

You can model your system as a hypothesis testing problem. For example, let $M_1$ and $M_2$ denote the two systems which are correlated. Let the observed signal in $M_1$ and $M_2$ be given as
\begin{align}
y_1(t) = s_1(t) + n_1(t) \\
y_2(t) = s_2(t) + n_2(t),
\end{align}
respectively, where $n_i(t) \sim \mathcal{N}(0, \sigma_i^2)$, for $i = \{1,2\}$, and $s_i(t)$, for $i = \{1,2\}$, is the original signal. Since the two systems are correlated, you can model them as
\begin{align}
\begin{bmatrix} y_1(t) \\ y_2(t) \end{bmatrix} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \begin{bmatrix} s_1(t) \\ s_2(t) \end{bmatrix} + \begin{bmatrix} n_1(t) \\ n_2(t) \end{bmatrix},
\end{align}
where $-1 < \rho < 1$ denotes the correlation constant between the two systems. Now, we can model the system as a hypothesis testing problem as follows:
\begin{align}
\mathcal{H}_0: \rho = 0, &\quad \mbox{uncorrelated model} \\
\mathcal{H}_1: \rho \ne 0, &\quad \mbox{correlated model}.
\end{align}
If you assume that $s_1(t)$ and $s_2(t)$ are known, then find the maximum likelihood estimate (MLE) of $\rho$ under $\mathcal{H}_1$. If you assume $s_1(t)$ and $s_2(t)$ are unknown, then use generalized likelihood ratio test, i.e., find the MLE of $s_i(t)$ under each hypothesis and replace $s_i(t)$ with its MLE.


convolution expression h(n,k) = h(n-k) As a new to DSP concepts  want to understand the meaning of impulse response and why we do $h(n-k)$ in convolution and $h(n+k)$ in correlation.

The correlation is used to find the time shift between signsls while convolution represents system response to predefined input.
Since the system response depends on previous input, rether then future input, the sign of the time distance is negative.


Possibility of temporal super resolution I've been reading some papers in the context of super resolution in imaging and video (for example works by Michal Irani [1], [2] and so forth), and was wondering if a similar idea can be applied with improving the temporal resolution of waveforms. In other words, if there are multiple recorders that are deterministically delayed from one another, will I be able to super-resolve in time a waveform? was this done?

Broadly similar techniques are used in a variety of applications when it's necessary to capture more samples of an event than is possible with available equipment.  The idea is really quite simple: an interpolation is just "stretching" but it does not provide any new data, which is why interpolation requires a low pass filter.  By contrast, if you add another recorder (AD converter) at a fixed delay (usually by setting up another record that is trigger by a phase shifted version of the master clock), then you are "subsampling" the analog data.

You can improve temporal resolution by using multiple ADs at phase shifts, capturing the same analog signal. You can improve spatial resolution using multiple transducers.  You can do both at the same time, although there are practical limits to how much certainty either measurement can have.

These techniques are used in probably all high end imaging: PET, airport scanners, CAT, medical ultrasound imaging, ground penetrating radar, etc.

Some of the discussion in those papers covers practical concerns that affect only video.


What does it mean that the normalized power of a received equals one? I have read in that lecture, page 2, here that,  the symbols are independent, and normalized to have power 1, which means $E(x x^H) = I$ .. what does that mean? could you please explain that for me?

In mundane words: the concept is very close to that of orthogonality. For two components of $x$, say  $x_i$ and  $x_j$, if they are different ($i \neq j)$, the average outcome of their product is zero (or $E(x_i x_j^H)  =0$). If   $i = j$, one assumes that the energy is normalized to one (or $E(x_i x_i^H)  =1$).

Resultingly, the product matrix is a unit matrix: one for the diagonal ($i = j$), and zero outside.


100 kHz sinewave contaminating a signal sampled at 48 kHz If sampling at a rate of 48 kHz, with no anti-aliasing filter on the input, and the input signal is contaminated by a 100 kHz sine, at what frequency would the sinewave signal appear in the output spectrum?

In your case, since 100kHz is more than twice 48 kHz, it will be folded twice in the sampled version, then there are two pairs of positions: +/-4 KHz and +/- 44Khz.


Graphical DT Convolution Attached is the question 1 part B, solution and my attempt. My attempt seems to be not correct. can someone comment?

To elaborate the point of my comment:

What you did there is a kind of circular convolution, where you thought of the impulse response $h(t)$ as periodic, i.e. the system triggered by $\delta(0)$ would oscillate until infinity. The periodicity is represented by your "wrapping around" the impulse response. The output you calculated is actually just one period of the output signal.

This is different from the linear convolution that the task formulation obviously expects to be applied here. With the linear convolution, the impulse response is thought to be the graph representing $h(n)$ padded by zeros from $-\infty$ to $0$ and from $3$ to $\infty$. Here, if the system is triggered by $\delta(0)$, you would get $h(n)$ and then all zeros to infinity.


OFDM Times of Arrival of multipath Assume i have received a known analog OFDM symbol (such the Short-Training-Field (STF) in WiFi systems), call it $s(t)$ with a known DFT per subcarrier, in the presence of multipath and absence of noise.

$$r(t) = \sum\limits_{i=1}^p \gamma_i s(t-\tau_i)$$
where $\gamma_i$ and $\tau_i$ are complex gains and the delays of the multipath, respectively. Without loss of generality, let $\tau_1 > \ldots > \tau_p$.

My question is simply, what is $\tau_1$? I could imagine what $\tau_2 \ldots \tau_p$ are. They are just relative delays with respect to $\tau_1$. But i don't seem to quite grasp how to quantify or understand what $\tau_1$ means. Is it absolute ? Is it with respect to something ?

Any suggestion is much appreciated.

Thanks

I don't fully understand your question, since you correctly state

  $\tau_i$ are delays of the multipath

therefore, $\tau_1$ is simply the delay incurred by the shortest path.

And no, $\tau_{i\ge2}$ are not relative delays; look at the formula: you just add up different independent paths with different delays; every path has its own delay, which is totally independent from $\tau_1$.


Bilateral Laplace transform and existence of Fourier transform I was reading from Athanosios Papoulis' "The Fourier integral and its applications." and they referenced the bilateral Laplace transform and Fourier Transform as:

$$F(p)=\int_{-\infty}^{\infty}e^{-pt}f(t)dt$$
$$F(\omega)=\int_{-\infty}^{\infty}e^{-j\omega t}f(t)dt$$

and stability indicates that the real part of $p$ must lie between $a$ an $b$.(For stability I'm guessing?)

As per the textbook, if we take the values of $p$ to be purely real and ignore the imaginary axis, then $F(\omega)$ doesn't even exist.

logically, I can't see how that could happen.
I was wondering if anyone could explain how that is so.
As well as why we limit the real part of $p$.(guessing it's similar to the unit circle in Z transform).

The bilateral Laplace transform converges in a vertical strip $a<\text{Re}\{p\}<b$, called the region of convergence (ROC). Compare this to the bilateral $\mathcal{Z}$-transform which converges in an annulus centered at the origin of the complex plane: $r_1<|z|<r_2$. For causal signals we have $b=\infty$ and $r_2=\infty$.

If the vertical strip $a<\text{Re}\{p\}<b$ does not include the imaginary axis, i.e., if $0<a<b$ or $a<b<0$, the bilateral Laplace transform does not converge for $p=j\omega$ because the imaginary axis is not inside the ROC. Consequently, the Fourier transform does not exist because the corresponding integral does not converge. For the $\mathcal{Z}$-transform the analogous case would be that the ROC does not include the unit circle, in which case the discrete time Fourier transform does not exist.


Implement the 3D Structure Tensor in C/C++? I am trying to figure out the details on how to implement the 3D structure tensor in C/C++ in an easy but efficient way and need some advice!

For a discrete function $ I(x_i,y_j,z_k)$ the 3D structure tensor is given by:
$$ S=\begin{pmatrix}
W \ast I_x^2 & W \ast (I_xI_y) & W \ast (I_xI_z)\\
W \ast (I_xI_y) & W \ast I_y^2 & W \ast I_yI_z \\
W \ast (I_xI_z) & W \ast I_yI_z & W \ast I_z^2 \\
\end{pmatrix}$$
where W is a smoothing kernel, $\ \ast $ denotes convolution and subscript denotes partial derivative with respect to.

The calculation of the structure tensor have two main steps:

calculate the partial derivatives of the function in a way that is
robust to noise over a window.
smooth products of the partial
derivatives over another larger window.

I start by looking at step 2.
I want to use a Gaussian as smoothing kernel. The normal distribution in 3 dimensions is separable:
$$ g(x,y,z) = g(x)g(y)g(z) $$
where
$$ g(x) = \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{1}{2}(x/\sigma)^2) $$
etc.

The Fourier transform of the normal distribution in 3 dimensions is given by:
$$ G(kx,ky,kz) = G(kx)G(ky)G(kz) $$
where
$$ G(kx) = \frac{\sigma}{\sqrt{2\pi}}exp(-\frac{1}{2}(kx\sigma)^2) $$
etc.

How do I implement the Gaussian smoothing?

The simplest way to implement the Gaussian smoothing would be to loop over a 3D Gaussian kernel for each point in I.

Since the Gaussian is separable however it should be more efficent to perform a convolution with a 1D Gaussian in the x direction followed by the y direction and z direction.

Another even more efficient approach (?) would be to do the convolution in the Fourier space (where it becomes a multiplication):
$$ g*a=\mathcal{F}^{-1}[GA] $$

Next I look at step 1.

The derivative of the Gaussian is given by:
$$ G_x(x) = -x/\sigma^2G(x)G(y)G(z)$$
etc., where subscript denotes partial derivative with respect to.

The Gaussian derivative can be used to estimate the partial derivatives of I in a way that is robust to noise:
$$ I_x=-x/\sigma^2G(x)*G(y)*G(z)*I $$

Here I am faced with the same decision:
implement it without using separability,
implement it using separability or implement it in Fourier space?

Any advice or comments?

There are many sides to this.

For typical kernel sizes, Fourier transform is not the most efficient method for convolution, if you are correctly utilizing the CPU for sliding window operations. That, of course, takes a good deal of software optimization. There are libraries to do that though (for example many of the deep learning libraries do 3D convolutions in CPU/GPU - I guess most of them do not use FFT). Given the specific symmetric nature of the Gaussian kernel, it is also very suitable for sliding operations. If the kernel size is very large, you could also try recursive convolutions.
The Fourier domain operations might be considered to be a better alternative if you have some luxury to pre-compute certain things. For example, if the image size is not changing and etc.
Probably the best answer would be to try and actually see it. I would prefer the separable convolutions. Vast amount of resources are also available for this.


Feature extraction for sound recognition and classification I am building an application that would "listen" to the microphone input, analyse it, and compare the analysis to a pre-analysed and pre-classified sound bank (small - maximum 20 sounds). It will then show the user what sound it was.

Now, I have a vague idea on how to implement this. I would like to choose a set of features that would best represent the sounds. The issue is that the sounds in the sound bank could be whatever the user recorded. From strong onsets and short sounds, to long onsets and long sounds.

The current features I'm thinking of are:

Spectral Centroid
Spectral Flux
Spectral Rolloff

What do you think? Would these be sufficient to properly classify the sound? Also, as these features output a single value for specific sound buffer, how would you go about handling the feature vector that represents the whole sound? I am using kNN for classification, and was wondering what's the best way to compare two feature vectors? would cross-correlation be a feasible technique?

Thanks a lot!

P.S
I have seen that a similar question was asked here, but it doesn't fully answer my issues.

Have a look @ librosa, a simple python library for audio analysis, implementing common features.  Here is a great introduction and example notebook.


How to add two signals? If x(t) is the curve shown in the image

how can I find (x(2t)+x(2t+2))?

I know how to find the individual terms but don't know how to add them!

And how can I express my answer in terms of the unit step function?

First express $x(t)$ in piece-wise form:

$$
x(t) =
\begin{cases}
0 ~&,  \text{ for}  & t  < -2 \\
0.5(t+2) &, \text{ for}  &-2 < t < 0 \\
-t+1 &, \text{ for}  & 0 < t < 1 \\
0 &, \text{ for }   &1 < t \\
\end{cases}
$$

Now, express $x_1(t)=x(2t)$ in a similar form:
$$
x_1(t) = x(2t) =
\begin{cases}
0 ~&,  \text{ for}  & 2t  < -2 \\
0.5(2t+2) &, \text{ for}  &-2 < 2t < 0 \\
-2t+1 &, \text{ for}  & 0 < 2t < 1 \\
0 &, \text{ for }   &1 < 2t \\
\end{cases}
$$

Recognize how the formula and intervals are adjusted for $x_1(t)$
$$
x_1(t) = x(2t) =
\begin{cases}
0 ~&,  \text{ for}  & t  < -1 \\
t+1 &, \text{ for}  &-1 < t < 0 \\
-2t+1 &, \text{ for}  & 0 < t < 0.5 \\
0 &, \text{ for }   &0.5 < t \\
\end{cases}
$$

Now, express $x_2(t)=x(2t+2)$ in a similar form:
$$
x_2(t) = x(2t+2) =
\begin{cases}
0 ~&,  \text{ for}  & 2t+2  < -2 \\
0.5(2t+2+2) &, \text{ for}  &-2 < 2t+2 < 0 \\
-(2t+2)+1 &, \text{ for}  & 0 < 2t+2 < 1 \\
0 &, \text{ for }   &1 < 2t+2 \\
\end{cases}
$$

Recognize how the formula and intervals are adjusted for $x_2(t)$ as well
$$
x_2(t) = x(2t+2) =
\begin{cases}
0 ~&,  \text{ for}  & t < -2 \\
t+2 &, \text{ for}  &-2 < t < -1 \\
-2t-1 &, \text{ for}  & -1 < t < -.50 \\
0 &, \text{ for }   &-0.5 < t \\
\end{cases}
$$

Now add those two piecewise defined signals $x_1(t)$ and $x_2(t)$, considering the intervals carefully:

$$ y(t) = x_1(t) + x_2(t) = x(2t) + x(2t + 2) =
\begin{cases}
0 ~&,  \text{ for}  & t < -2 \\
t+2 &, \text{ for}  &-2 < t < -1 \\
-t&, \text{ for}  & -1 < t < -0.5 \\
t+1 &, \text{ for }   &-0.5 < t < 0 \\
-2t+1 &, \text{ for}  & 0 < t < 0.5 \\
0 &, \text{ for }   &0.5 < t \\
\end{cases}
$$

Expressing this using unit steps (alone) is not possible, you need the ramp functions. Then you can provide a representation using scaled and shifted ramps.

A MATLAB plot of the resulting signal is provided below ;


FIR filters always linear? Are Finite Impulse Response filters always linear? I looked at https://en.wikipedia.org/wiki/Nonlinear_filter and its says that if a filter does not obey linearity, it is non-linear filter. Thats understood? But does does linearity has to do with the length of impulse response?

What am I missing here?

Yes, FIR filters are linear.  No, the length of the impulse response does not have anything to do with whether it is linear or not.  For instance, IIR filters are also linear.

Linearity is determined by whether the filter obeys the superposition principle, nothing else.  So, FIR and IIR filters are linear, median filters are not.

EDIT: Technically Dilip is correct.  Finite and infinite impulse response filters can be non-linear if their impulse responses incorporate non-linear operations, such as squaring.  I was going by the "industry" meaning of FIR and IIR, which are linear because they only use multiplication and adding, but in the mathematical sense FIR and IIR filters can be non-linear.


Power spectral density of sigma-delta modulator I am simulating a 3rd order SDM and I want to plot the PSD of the output.

As you may know, the output is a sequence of pulses ranging from -3 to +4.
The number of samples is $N$.

The way I approached the problem is the following:

Use fft on the output sequence
normalize (divide by $N$)
Square the above signal.
multiply by 2 ( because I need the one-sided band -- not sure about this ).

The plot should then be $10\log_{10}(signal)$ if I am not mistaken.

Problem is that I get a big deviation from the ideal curve of the expected noise.

Am I doing something wrong?

Everything you have written sounds correct so I assume you forgot to consider the BW for each FFT bin and are comparing to an expected spectral density that is per Hz.  Without windowing or zero padding the BW in each bin will be 1/T where T is the time length of your signal (so in other words 1 bin wide).  Therefore you want to reduce the total power measured in each bin by 10Log(1/T) to get the equivalent power per Hz spectral density.

If you are using windows you can refer to fred harris' classis paper on windowing for the Equivalent Noise BW of all common bins as well as windowing implementation loss for an accurate computation.

https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf


Is there a set formula for calculating FIR coefficients (for a LPF ir 96 taps)? I need to calculate a series of coefficients for several FIR filters having 96 taps is there a standard or set formula for calculating these, it is a low pass filter.

is there a standard or set formula for calculating these, it is a low pass filter.

No. There's very many approaches to filter design, and it very much depends on what you want to achieve with your filter.

Things like:

What do you need as a low pass filter? Is it important to have linear phase?
Is it going to be narrow, or wide? (i.e. is the stopband bigger or smaller than half of the input bandwidth)
Is it going to be implemented in hard- or software?
What about numerical precision? If software: fixed or floating point? SIMD'ability? Enforcing linear phase to be able to exploit coefficient symmetry?
What's more important:

Passband ripple?
strength of the strongest sidelobe?
Cumulative stopband power leakage?
Maximum steepness for given length?
Numerical stability (i.e. well-conditioned taps)?
...

Now, assuming you "just" want a FIR for usage with e.g. Python, you go ahead and use your favourite scientific/signal processing toolkit, and use one of the many window design methods. Often, things will get automatically chosen correctly for you. In essence, scipy has an extensive documentation on its filter design routines, and I urge you to look at that.


How can I describe a sum of sinusoids by an equivalent lowpass signal? I consider a sinusoidal signal with fixed amplitude $A_{\rm in}$ and frequency $f_c$.

This signal is sent into a nonlinear amplifier, whose AM/AM and AM/PM characteristics are known. Using these characteristics for the given amplitude, I can determine the output amplitude $A_{\rm out}$ and the output phase $\phi_{\rm out}$.

Now however, I consider a sum of two sinusoids with different amplitudes $A_{\rm in_1}$ and $A_{\rm in_2}$ and different frequencies $f_{\rm c_1}$ and $f_{\rm c_2}$. First, I want to describe this input signal by its equivalent lowpass signal. I want to do this in MATLAB. I know the hilbert() command gives the analytic signal. So I suppose when my first sinusoid is called $x_1(t)$, the equivalent lowpass signal can be described by $\mathcal H(x_1(t)) \cdot \exp(j2\pi f_{c_1} t)$. The same is valid for the second sinusoid. But  I don't know how to represent the sum by an equivalent lowpass signal. In a second step, I want to apply the AM/AM and AM/PM characteristics to the signal. But now, I have two different amplitudes.

I found this:
Consider a multiple carrier input:
$x(t) = \sum(A_i(t)*\cos[(\omega_c + \omega_i)* t + \theta_i (t)],i,N)$

This can be written as $x(t) = A(t)*\cos[\omega_c *t + \theta(t)]$

where
 $A(t) = [x_I^2(t) + x_q^2(t)]^{1/2}$

$\theta(t) = tan^{-1} [x_I(t) / x_q(t)]$
I really don't understand this equation given for the resulting phase. Why is it valid? I mean X_I and X_Q are not just real and imaginary parts of x(t), are they? I know by Euler's formula, that the cos corresponds to the real part of the exponential and the sin corresponds to the imaginary part. But (see below) X_I and X_q are both sums of cos or sin waves of different amplitudes and different phases. How can I just add them?

$X_I(t) = \sum(A_i(t) * \cos[\omega_i * t + \theta_i (t)],i,N)$

$X_q(t) = \sum(A_i(t) * \sin[\omega_i * t + \theta_i (t)],i,N)$

The sum of two sinosoids with different frequencies $f_{c1}, f_{c2}$ does not add up to a single sinosoid. Only, if both sinosoids have the same frequency, you can combine them into a single sinosoid using the equation you state at the end of your post. See also for example this post.


Phase rotation without Fourier Transform Is there any known math equation or other method to perform phase rotation on signal without decomposing it to frequency domain?
In frequency domain it is obvious. You just need to perform phase shift for each frequency.
But decomposing signal from time domain to freq domain, then perform phase shift on each frequency, then again composing it to time domain it is very computationally demanding for frequency-rich signal.
That's why I wonder if I can achieve it without decomposing signal to frequency domain?
To be clear I want to perform on the signal something like that:

But I don't know the frequencies in the signal. I just have time domain signal. Is that possible?

You can use the Hilbert transform, then multiply the real part with the sine of the angle you wish to transform, and the imaginary part wth the cosine. A quick code in Octave:
t=[0:0.01:4];
s=hilbert(1-2*mod(t,1));
a=real(s);
b=imag(s);
phi=[0:0.5:6];
test=a.*sin(phi)+b.*cos(phi);
plot(a+4,"",test');

I know there are Hilbert transforms made with IIRs, but their phase is not exactly linear, so a FIR is a better, but costlier choice. Its order will be a function of the lowest frequency of interest, which will count as the transition width. If your bandwidth is large, the order will be large.

A FIR is a filter, which means it acts in time domain (but, implicitly, the frequency domain will be affected). To give an example, the signal you're showing of 1 Hz, sampled at 32 Hz, filtered with a FIR with a Kaiser window with As=60 and 1 Hz transition width:

The output isn't sampled because I am brute-adding the two inputs, but that's how the result would look like. The green number is the order. It's a symmetrical FIR so you'll need that many delays and half as many multipliers. You could try a polyphase approach and that will reduce the number of delays, but not the multipliers.


Delaunay triangulation and Voronoi diagram in image processing I am reading about image processing and it was mentioned in a couple of articles "Delaunay triangulation" and "Voronoi diagram".

I read on what the two things mean and how they are acquired, but why are they important in image processing? After searching on the DSP stackexchange and overall the Internet, I did not seem to find an explanation about this, so if I have missed a thread, with an explanation, I am apologizing in advanced.

...why are they important in image processing?

Because they offer ways by which to summarise a large set of points.

The Voronoi Diagram (VD) and Delaunay Triangulation (DT) of a set of points are duals of each other. If you have one of them you can construct the other.

The VD splits space into a set of tiles, one for each point in the set. Each tile represents those locations in space that are closest to a particular point in the set.

Imagine for example, that you take a handful of marbles and scatter them across the floor in a room. Then, you divide the room in a fine grid of locations. Let's say every 10cm x 10cm. You select one of the locations on this grid and measure its distance to all the marbles across the floor. This location will be closest to one of the scattered marbles. So, you assign this location in space to that marble. Then you repeat this process for all the locations in the fine grid, progressively assigning them to closest marbles.

In the end, you will have the VD of the set of points (positions of marbles on the floor). To visualise it, you simply plot each point in space choosing a different colour for each cluster.

To now construct the DT from the VD, you simply connect every two points in the set (the marbles on the floor) that share a boundary, to produce the VD's dual.

The whole point about this process is that it provides a way to relate the points in our set through rules and introduce some structure where before there was none.

The rule that governs the VD (measure distance, assign membership to the shortest distance) is "translated" to another rule in the DT. The triangles that are formed by connecting VD points that share a boundary have the maximum smallest angle of all other possible ways of connecting the points in triangles. Even for a moderate size of a set of points there is a huge number of possible ways by which the points can be connected into triangles. If you did form all of them and measured the smallest angle in the triangles of each one of these ways, the DT would be the one where the smallest angle of its triangles would be the largest.

The most obvious example where Delaunay Triangulation is applied is the construction of mesh geometry from point clouds. Say for example that you have a lidar which gives you a point cloud of measurements to points in its environment. If you do a Delaunay Triangulation on that set of points, you get a mesh that you can then render.

If you don't have a lidar, perhaps you have a number of cameras that are taking pictures of the same object and you wish to reconstruct the 3D representation of the object from the 2D views of each camera.

In these cases, the connection between DT and the result is obvious: We have a set of points and we would like to find the surface they are likely to define.

But, there are other applications where the "summarisation" is more pronounced.

Say for example that you want to do object recognition that is independent from scale or resolution. You take a camera, shoot a picture and then obtain some set of descriptors off of the object.

The descriptors are simply a cloud of points with specific geometrical relationships between them. Whether you shoot the object up close or from (reasonably) far away, these relationships will be preserved.

To summarise these points, you might want to use DT. In fact, it is used in face recognition or similar pattern recognition and shape-based retrieval applications as a descriptor of a set of points. In those, the DT is not just a way of connecting the points together but because of its properties (such as the minimal angle mentioned before) it can also be used to extract features from the resulting DT.

The VD has had applications in image compression. The main idea there is that you are trying to represent a large area of homogeneous (or approximately homogeneous) pixels with the value of a single "dominant" pixel. So, instead of storing $M \times N$ pixels, you store a set of key point positions and their colour and reconstruction involves obtaining the VD of the set of points which is now "painted" with the pixel colours.

Here (or here) is an example of that.

The artistic aspect of VDs (or DTs) applied to images is a consistent fad, it comes back again and again, you basically reduce the amount of triangles drastically and it creates a sort of "stained glass effect" on the whole image.

Hope this helps.


Why Does the Kalman Filter Remove Only Gaussian Noise? What and where in the derivation of the Kalman filter is the assumption of Gaussian noise? Why and how does this assumption help?

A Gaussian process is completely specified by its mean and variance.  A Kalman Filter updates the process mean which is the state, and its variance.  These are the sufficient statistics.

The measurement noise is reduced but the process noise  is part of the recursive state history and is tracked.

Your heading question and subsequent paragraph aren't fully consistent.  One can use a linear Kalman Filter when the noise isn't Gaussian in many circumstances but it wouldn' be optimal.  The derivation assumes Gaussian (or Normal) noise and all deterministic inputs are known, as well as knowing initial values.


power spectrum of passband of ofdm is it right shape or not? I tried to get power spectrum of ofdm passband signal.
Below is the whole matlab code for that with BPSK+ofdm.
and I got the power spectrum as below.
I wonder that is that right form or not?

Nfft=64; % number of sub carriers
Ncp=Nfft/4; % number of cyclic prefix
Nuc=52; % number of used carriers

d=2*(randn(1,Nuc)>0.5)-1; % bpsk mapped data

X=zeros(1,Nfft);
X([2:Nuc/2+1,(end-Nuc/2)+1:end])=d; % only used sub carrier
x=ifft(X,Nfft); % time domain
sbb=[x(end-Ncp+1:end) x]; % to add cyclic prefix

%% passband

fc=4; % carrier frequency
sps=16; % samples per symbol
Tsym=1; % symbol duration time
Ts=Tsym/sps; % sampling time
fs=1/Ts;

Tmax=(Nfft+Ncp)*sps*Ts-Ts;
t=0:Ts:Tmax;

%return

% I-ch up-sampling
I=real(sbb);
I=repmat(I,sps,1);
I=I(:).';
sI=I.*cos(2*pi*fc*t);

% Q-ch up-sampling
Q=imag(sbb);
Q=repmat(Q,sps,1);
Q=Q(:).';
sQ=Q.*-sin(2*pi*fc*t);

s=sI+sQ; % transimt signal

%% plot power spectrum

f=(0:(numel(s)-1))/numel(s)*fs;

figure
plot(f,10*log10(abs(fft(s))))
%plot(f/1e6,10*log10(abs(fft(s))))
xlabel('frequency[Hz]');
ylim([-40 20]);
title('full span');
grid on

figure
plot(f,10*log10(abs(fft(s))))
%plot(f/1e6,10*log10(abs(fft(s))))
xlabel('frequency[Hz]');
ylim([-40 20]);
xlim([3 5])
title('main band');
grid on

Doesn't look right. a 64-FFT with only the center 52 carriers occupied should look pretty "rectangular".
I'd recommend not doing the upconversion. OFDM is designed to be an equivalent-baseband system, so you can definitely investigate whether your OFDM works directly in baseband. It's way more intuitive that way.


Transfer Function Block Diagram Confirmation Can someone confirm if this is the correct block diagram for the following transfer function?

The original equation provided was: y[n+1] = y[n] + 0.01x[n]

Which I rearranged into H(z) = Y(z)/X(z) = 0.01 / (z - 1) and then also rearranged into y(n) = y(n-1) + 0.01x(n-1)

And left me with the following block diagram,

That is a correct way to do it, but leaves you with (depending on what you're trying to do) an excess of delay blocks.  This also works:


How can frequency dispersion in a vibrating string be modelled using an LTI filter? According to JOS on stiff-string synthesis, stiff strings (like on a piano) introduce inharmonicity (i.e. the harmonics of the tone are not all in tune) due to different frequency components of the wave travelling at different speeds. This inharmonicity is called "dispersion".

According to another page on the same site about modelling this effect, he says it can be done by putting an LTI (linear, time-invariant) filter of some sort at one end of the delay line (in waveguide modelling, a delay line is used to model a physical string) before feeding the output back into the delay line.

Somehow, this LTI filter introduces inharmonicity, or frequency distortion. But by definition, an LTI filter does not introduce any new frequencies into the signal (and thus can be accurately represented by a frequency-response graph).

How can inharmonic tones be introduced by an LTI filter? Am I missing something obvious? Do I not understand my DSP basics well enough?

From your question, it sounds like you think the all-pass filter is an independent post-processing step that takes as its input the output of the classic Karplus-Strong model (delay line, LP filter, close the loop). If that would be the case, yes, it would be very surprising indeed that new harmonics would appear! But here, the allpass filter is part of the closed loop. You can intuitively think about what it does in these terms: it is a short delay, with a delay value which is frequency-dependent (phase shifter). This gets added to the delay provided by the main $N$ samples delay line, and causes the lower frequencies to "see" a longer loop. So the highest harmonics will have a period close to $N$, and the lower harmonics will have a lower period. The overall transfer function of the loop is no longer a comb with regular spacing between teeth, but a comb with a slowly increasing gap between teeth.

What you are forgetting also is that this whole system (delay line, all-pass filter combined with loss LP filter, close the loop) which is a LTI, is excited by a burst of white noise. So the input signal contains all frequencies, and there's nothing surprising or contradictory with the fact that altering the transfer function of a LTI can change its frequency response so that its resonant frequencies are shifted, or so that new bumps appear.

If you need to convince yourself... System A is a band-pass filter with a frequency of 1kHz. System B is a band-pass filter with a frequency of 1kHz ; in parallel with a band-pass filter with a frequency of 2kHz. Both are excited by white noise. In the second case you have "added" a LTI filter in the system; and a new frequency has appeared. There's no magic in this - you have just shaped the frequency response of the system differently to sculpt the input in new ways.


Optimal receiver for a binary system employing waveform $s=A\cos^{2}(\pi f_c t \pm \alpha)$ I want to design an optimal (minimizing error probability) receiver for system of two waveforms (with $0<\alpha<\frac{\pi}{2}$):
\begin{align}
s_1 &=A\cos^{2}(\pi f_c t + \alpha)\\
s_2 &=A\cos^{2}(\pi f_c t - \alpha)
\end{align}
It seems like a modified BPSK, as:
\begin{align}
s_1&=A\cdot\ \frac{1 + \cos(2\pi f_c t + 2\alpha)}{2}\\
s_2&=A\cdot\ \frac{1 + \cos(2\pi f_c t - 2\alpha)}{2}
\end{align}
But I don't know what to do next.

Awesome, your transformation into a sum of a constant and a cosine is wonderful!
Since both signals are cosines shifted from the origin by the same $\frac A2$, that doesn't chang their euclidean distance.
Thus, subtracting $\frac A2$ doesn't change your decider's error rate. Also, it's technically super easy – just a high pass filter that removes the constant, no matter what it was.
The you're left with
\begin{align}
\tilde s_1 &= \frac A2 \cos(2\pi f_c t + 2\alpha)\\
\tilde s_2 &= \frac A2 \cos(2\pi f_c t - 2\alpha)& \beta := 2\alpha\\[1.5em]
\tilde s_1 &= \frac A2 \cos(2\pi f_c t + \beta)\\
\tilde s_2 &= \frac A2 \cos(2\pi f_c t - \beta),&0<\beta<\pi
\end{align}
Now, IQ-mixing this down with $f_c$ will indeed give us some phase-modulation-type thing, with the two constellation points with $\pm \beta$ as phase – but that's not a BPSK, since that is not always a phase difference of 180°!
However, it's always symmetrical to the I-axis in your constellation diagram. So, you can just remove the real part of your baseband signal and end up with something where the sign contains all the information.
Note that this, geometrically, is a Maximum Likelihood decision only if your noise distribution (which haven't told us) is is symmetrical after the transformation – which is the case for the usual complex Gaussian noise we assume, but I'd be careful whether that's the right assumption here - such non-symmetric PSK constellations are (surprisingly) actually used in a couple of cases, most of which involve interesting residual carrier systems with nonlinearities that shape the noise in a non-symmetric way.


Shift signal in frequency domain I want to shift the signal in frequency domain without going back to the time domain.

Is it enough to shift the samples of the spectrum only? How?

If you merely want to shift spectral samples up in frequency, the process is fairly simple.  Let's say your DFT comprises the eight complex-valued samples X[0], X[1], X[2], X[3], X[4], X[5], X[6], X[7] and you want to shift them upward in frequency by two samples. The new shifted Y[m] spectral samples would be:

Y[0] = X[6]

Y[1] = X[7]

Y[2] = X[0]

Y[3] = X[1]

Y[4] = X[2]

Y[5] = X[3]

Y[6] = X[4]

Y[7] = X[5]

You'll notice that the shifting could be described as a "circular shift" because of the sample wrap-around behavior. If you want to shift your original X[m] sequence downward in frequency by three samples. The new shifted Y[m] spectral samples would be:

Y[0] = X[3]

Y[1] = X[4]

Y[2] = X[5]

Y[3] = X[6]

Y[4] = X[7]

Y[5] = X[0]

Y[6] = X[1]

Y[7] = X[2]


Z-transform of difference equations and stability of a process According to this paper:

  $y(t)$ is stationary if all of the roots (of characteristic equation) lie outside the unit circle

Here, $y(t)$ is causal.

To me it seems the case is exactly the opposite, the roots should be inside the unit circle for causal process to be stable.

Example, consider the difference equation:

$$y_t-5y_{t-1}+6y_{t-2}=0$$

The characteristic equation is then: $$r^2-5r+6=0$$

Solution for characteristic roots: $$r=2, r=3$$

Solution for the difference equation (assuming 1 as initial conditions): $$y(t)=2^t+3^t$$

Clearly the roots are outside the unit circle (as 1<2,3), and the process unstationary (as it grows with t). To me it seems that the 3rd equation in the paper should have negative exponent for zs (in which case the solution would be the same as mine, even if arrived through the z transform).

Yet, the claim appears in many places so it can't be wrong (I doubt the paper is wrong either). However, for LTI system to be stable the poles must be inside the circle, a other claim that is repeated everywhere. Is the statement just plain wrong or are these things somehow not equivalent?

The paper deals in the lag operator $L$, rather than $z^{-1}$ as DSP people tend to use.

As a result, all the DSP results you know are inverted. So "inside the unit circle" becomes "outside the unit circle", black is white, and ... avoid zebra crossings.


How to design a filter with a certain magnitude response I am trying to design a filter whose magnitude is the same as that of a given signal. The given signal is wind turbine noise, so it has significant low-frequency content. After designing the filter, I want to filter white Gaussian noise so as to create a model of wind turbine noise. The two signals, that is the original and the filtered noise should sound similar.

I am using arbitrary magnitude filter design in Matlab for that (FIR, Order: 900, Single rate, 1-band, response specified by amplitudes, Sample rate 44100 Hz, Design method: firls). The problem is that, although I design the filter using the values from the original signal's magnitude, the filter magnitude fails to follow the magnitude at higher frequencies. Could you please help me with that?

The idea is that I am using a polynomial curve, to fit the shape of the spectrum of the original sound. Then, I filter the magnitude of the generated noise, by multiplying the extracted polynomial curve with the generated noise magnitude. Finally, after calculating the new magnitude and phase, I get back to the time domain.

[x,fs] = audioread('cotton_0115-0145.wav');  % original noise sample
x = x(:,1);                                  % extract one channel
x = x.';
N = length (x);

% fft of original signal
Y = fft(x,N)/N;
fy = (0:N-1)'*fs/N;

% half-bandwidth selection for original signal
mag = abs(Y(1:N/2));
fmag = fy(1:N/2);

% polynomial fitting
degreesOfFreedom=40;
tempMag=20*log10(mag)';  % desired magnitude in dB
tempFmag=fmag;

figure(1)
plot(tempFmag,tempMag,'b');
title('Spectrum of original signal-Polynomial fitting')
xlabel('Frequency (Hz)');
ylabel('20log10(abs(fft(x)))');
axis([tempFmag(1) tempFmag(end) min(tempMag) 0]);

hold on
p = [];

for i=1:4
    p=polyfit(tempFmag,tempMag,degreesOfFreedom);
    pmag=polyval(p,tempFmag);

    plot(tempFmag,pmag,'r');
    pause

    above=pmag<tempMag;
    abovemag=tempMag(above);
    abovefmag=tempFmag(above);

    tempMag=abovemag;
    tempFmag=abovefmag;

end

hold on
legend('signal magnitude','polynomial');

%
loc1 = find(fmag == 0);
loc2 = find(fmag == 22050);
Nmag = length(mag);

M=((Nmag-1)*max(tempFmag))/(tempFmag(end)-tempFmag(1));
freqFinal=tempFmag(1):max(tempFmag)/M:max(tempFmag);
freqFinal=tempFmag(1):max(tempFmag)/length(mag):max(tempFmag);
magnitudesFinal=polyval(p,freqFinal);

figure(2)
plot(fmag,20*log10(mag)');
hold on;
plot(freqFinal,magnitudesFinal,'g');
title('Spectrum of original signal-Choice of polynomial curve')
xlabel('Frequency (Hz)');
ylabel('abs(fft(x))');
axis([freqFinal(1) freqFinal(end) min(magnitudesFinal) max(magnitudesFinal)]);

%%
% noise generation
Nn  = N;
noise =  wgn(1,Nn,0);
noise = noise/(max(abs(noise)));
Ynoise = fft(noise,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;

% polynomial for whole f range
newmagA = 10.^(magnitudesFinal/20);
newmagB = fliplr(newmagA);
newmagB(end+1) = newmagB(end);
newmag = [newmagA newmagB];

%filtering
Ynoisenew = newmag .* Ynoise;
figure(3)
magnoise = 20*log10(abs(Ynoisenew(1:Nn/2)));
fnoise = fn(1:Nn/2);
plot(fnoise, magnoise);

% magnitude and phase of filtered noise
magn = abs(Ynoisenew);
phn = unwrap(angle(Ynoisenew));

% Back to Time domain
sig_new=real(ifft(magn.*exp(1i*phn)));

figure(4)
sig_new = sig_new / max(abs(sig_new));
plot(t, sig_new);

Ysignew = fft(sig_new,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;
figure(5); plot(fn(1:Nn/2),20*log10(abs(Ysignew(1:Nn/2))));

if you are only interested in designing a signal representing the wind turbine noise, you could just generate it from your magnitude. E.g., with your magnitude $A$ over frequency $f$ you can generate a random phase $\phi$ for each frequency and over a time $t$ and then add all together. Here a small code for illustration (I tried to use your varibles. However you need to check orientation etc.):

phi         = rand(size(fmag))*2*pi;
x           = (mag*sin(2*pi*fmag*t+repmat(phi,1,N)));

Hope this helps.


Time duration of a signal to determine the peak frequency I have two questions regarding signal processing:

1) In my analysis I generate a signal in the time domain with a frequency of 0.1 Hz. I apply the FFT and then determine the peak frequency of the signal, which should be 0.1 Hz. But I notice that when I simulate my signal for 200s, the peak frequency is always close to the 0.1Hz, but when I simulate my signal for 20s, then the peak frequency is not always 0.1Hz. Is there a certain rule of thumb for the duration of the signal to have the peak frequency occur in the frequency domain?

2) Currently I am just determining the peak frequency of my signal, as the highest peak in the frequency domain. But sometimes I have a lot of peaks near each other and at times they might overlap. What are some other recommended techniques to determine the peak frequency? I am using Matlab

Thanks!

What window function are you using when you apply your FFT?

If the signal being analyzed is simply truncated on both ends, the window is rectangular.

Lets call the length of time that you generate your signal the analysis window and call the signal you create the "windowed signal".  In this case I'm assuming the window to be rectangular.

One interpretation of the math behind the FFT is that the signal being analyzed is periodic and the complete version of the signal is formed by taking the unaltered input from the analysis window and placing it end to end so it repeats indefinitely as illustrated below

....| Windowed Signal (-1) |  Windowed Signal (0) | Windowed Signal (+1) |...

Depending on where you stop your signal generation, you may introduce discontinuities in the periodic signal that the FFT assumes when you do your analysis.  For example if your signal generation stops at the peak of a sine wave but starts at a 0 crossing, the signal at the point where one windowed signal segment meets the nest will abruptly go from Vpeak to 0.  This discontinuity introduces spurious frequency content in the FFT.  Also if you don't make your analysis window an integer multiple of the period of signal being analyzed, the FFT will assume the signal is periodic in a different period than the signal being analyzed (actually the period assumed by the FFT is always the window length).  this will create frequency lobes that are not harmonically related to the signal you are testing.

As you increase the length of your window, the contribution of the discontinuity and window period diminish, so the longer duration signal will give more accurate results.

The way to improve accuracy of the FFT frequency and/or amplitude resolution is with the use of non-rectangular window functions.  There are a bunch of window functions that are designed for this purpose.  The window function is multiplied by your time domain data (signal).  The effect of the window function is to taper off the amplitude of at the ends of your analysis window (the start and end of the signal you generate) so the edge effects are reduced in the FFT computations.  Different windows have different advantages in frequency or amplitude resolution, so there are trade offs in choosing a window function.

In general, if you want decent resolution in your FFT, you should provide a data length that is least 2x the period of the lowest frequency to be represented and apply a window function to soften the effects of signal truncation.  Note that the longer you make your time window, the greater the frequency resolution of your FFT.  The bins are spaced at 1/T and its multiples where T is the duration of your signal.

There is plenty of information on the internet that discusses the different window functions  used in FFT analysis.  Just Google "FFT window functions".


Real valued FFT in Matlab I'm trying to compute FFT of size (N/2) for sequence of real values of length N. I'm using this recipe.
Here is my code in matlab:

clear;
clc;
x = [1, 2, 3, 4, 5, 6, 7, 8];
n = length(x);
z = [1 + 2j, 3 + 4j, 5 + 6j, 7 + 8j];
Z = fft(z);
Fe = (Z + (conj(Z)))/2;
Fo = (-1j*(Z - (conj(Z))))/2;
for k = 1 : n/2-1
     F(k) = Fe(k) + exp(-1j * 2 * pi * ((k-1)/n)) * Fo(k);
end;
F(n/2) = 0.5*(Z(1) + conj(Z(1)) + 1i*(Z(1) - conj(Z(1))))
fft(x)

Obviously, I get wrong result (at least I think it's wrong). As far as I understand, first four terms of F and fft(x) should be equal.

The part that confuses me is when I compute Fe and Fo I'm using fliplr(conj()) to get Z(N/2) - k)*, I'm not sure, if it's the proper way to do this.

You can do it like this. I think the code speaks for itself:

x = [1, 2, 3, 4, 5, 6, 7, 8];   % some real-valued signal
n = length(x);                  % must be even!
n2 = n/2;                       % assume n is even
z = x(1:2:n)+j*x(2:2:n);        % complex signal of length n/2
Z = fft(z);
Ze = .5*( Z + conj([Z(1),Z(n2:-1:2)]) );        % even part
Zo = -.5*j*( Z - conj([Z(1),Z(n2:-1:2)]) );     % odd part
X = [Ze,Ze(1)] + exp(-j*2*pi/n*(0:n2)).*[Zo,Zo(1)];  % combine
X2 = fft(x);
max(abs(X-X2(1:n2+1)))  % 2.4492e-15


From IQ to pi/4 DQPSK I need to decode $\frac\pi4$-DQPSK packets from an SDR, preferably in C#.

I have the SDR receiver part done. Now I'm trying to learn how to take the IQ stream and get the packet info from it.

I've gone through the math, but don't understand the notation. I've tried GNU Radio, but can't get blocks added on a windows machine. I can't find a library that I can get working. I know this is rather opened ended question, but can anyone steer me in a direction to do this? I've been programming for going on 42 years now, but this is my first foray into DSP, and at my age it doesn't come easy anymore.

If it helps any, I do know a lot about the packets I'm trying to grab. They are transmitted in pi/4 DQPSK at 16kbps or 32kbps (there are two packet types), with a preamble of f337eeb637660672. The packet format is known to me, and is not encrypted.  This is for a project that will be released for free, which combined with my being retired hampers any ability to pay for professional help.

The OP stated he was interested in $\pi/4$-DQPSK (not QPSK), so phase synchronization is presumably not an issue for him.

As far as the actual implementation is concerned, you'll save yourself some time if you become familiar with the bottom of pages 29 (symbol mapping) and 37 (differential detection) in this student paper. Ignore all the old TI DSP-chip specific information and the information on pulse-shaping, which is already built into GNU Radio (if the OP is using that). There is also a reasonable description of a simple timing recovery algorithm. It is serviceable in some scenarios, but not others.