texts:
I am having a problem understanding how the author gets the onset envelope of an audio song in this paper.

The entire section in question has been quoted below:

  First the input sound is resampled to 8 kHz, then we calculate the
  short-term Fourier transform (STFT) magnitude (spectrogram) using 32
  ms windows and 4 ms advance between frames. This is then converted to
  an approximate auditory representation by mapping to 40 Mel bands via
  a weighted summing of the spectrogram values [Ellis, 2005]. We use an
  auditory frequency scale in an effort to balance the perceptual
  importance of each frequency band. The Mel spectrogram is converted to
  dB, and the first-order difference along time is calculated in each
  band. Negative values are set to zero (half-wave rectification),
  then the remaining, positive differences are summed across all
  frequency bands. This signal is passed through a high-pass filter with
  a cutoff around 0.4 Hz to make it locally zero-mean, and smoothed by
  convolving with a Gaussian envelope about 20 ms wide. This gives a
  one- dimensional onset strength envelope as a function of time that
  responds to proportional increase in energy summed across
  approximately auditory frequency bands

My understanding of the procedure of the first part of the procedure is that our 44.1KHz audio signal is resampled to 8KHz and then we perform a overlapping widowed STFT operation on it so we get the Fourier transform coefficients.

But I really have no idea what is happening in the second part of the procedure:

I think we then alter those coefficients in such a way that the new Fourier transform coefficients are now balanced perceptually(ie the coefficients are now what a person's ear would perceive them to be).Then the so-called resultant "Mel spectrogram" is converted to dB, and the first-order difference along time is calculated in each band.Negative values are then set to zero (half-wave rectification),then the remaining, positive differences are summed across all frequency bands.

Why and how would we calculate the first order difference? What use is this? I have no extensive knowledge of the Mel scale and it's so-called "frequency bands" so this part is proving hard to understand.

This signal is passed through a high-pass filter with a cutoff around 0.4 Hz to make it locally zero-mean, and smoothed by convolving with a Gaussian envelope about 20 ms wide. This gives a one- dimensional onset strength envelope as a function of time that responds to proportional increase in energy summed across approximately auditory frequency bands

Why do we need the signal to be locally zero-mean?


texts: Why is it better to sample white noise in frequency domain here?

texts:
We have a channel $c(t)$ that we want to sample. It can be expressed as $ c(t)= g\ast x(t)$ where $x(t)$ is white complex Gaussian noise and $g(t)$ is a function that we know the PSD of, let's call it $G(f)$.

So two options  if we want to generate samples of the channel response:

Let $g(t) = \mathcal{F}^{-1}[G(f)]$, then generate white noise $x(t)$ and let $c(t) = g\ast x(t)$
Since $C(f) = G(f)X(f)$, we instead generate the  white noise $X(f)$, mulitply $G(f)$ and $X(f)$ and let $c(t) = \mathcal{F}^{-1}[C(f)]$

Apart from complexity issues, I wouldn't guess these methods would be too different. But at least for the PSD I tried it for (the Clarke's spectrum bathtub shape), there definitely was a difference at least when the number of samples $N_s$ was small. Sampling in frequency domain was better for small $N_s$. For large $N_s$ it didn't seem to matter.

So my question is, why does it matter? Does it have anything to do with my PSD getting asymptotic around the edges and that this will not be noticed when sampling?


texts: Is it possible to estimate variance of noise for a step answer signal?

texts:
I know there is not possible to find the true noise of a measured signal. The only way to "find" the noise is to estimate the noise. Noise has the mean 0, but the variance varies.

So assume that we have a signal which look like this:

The clean signal is $y(t)$ and the noise is $e(t)$. In this case, $e(t)$ has the variance 1 because it's generated by Octave/MATLAB function 0.05*randn().

So my question is:
How can i estimate variance of this noise from this signal?


texts:
I'm new to signal processing and I'm working with librosa library quite a bit.

I would like to know what is the information that is returned by librosa.load function when I read a audio(.wav) file using it. Is it the instantaneous sound pressure in pa? or the just the instantaneous amplitude of the sound signal with no unit?
How can I use this information to calculate the SPL(no-weigting) of the signal ?


texts: Identifying license plate number of a thief's car from a video

texts:
Motivation:

I recently had my laptop stolen, with lots of pictures of my family, research work and other very important stuff. Of course the theft has been reported, but I think people here are more competent in image processing compared to the local police station :)

Video Recordings

The stealing was captured by a security camera, original .dav files can be downloaded from here:

ch10.dav
ch13.dav

Converted .avi files can be downloaded (or played in browser) from here:

ch10.avi
ch13.avi

Description

In the file named ch10 (full name: HCVR_ch10_main_20160501143700_20160501143947) the thief's car enters the frame from the right at 0:48 . In the file named ch13 (full name: HCVR_ch13_main_20160501143811_20160501144205), in the first 2 minutes, the thief's car enters the frame from the right, makes U-turn, steals my stuff,  makes a U-turn again and go.

What has been tried?

I tried cropping a few frames around the license plate, stretch so that the plate will be on same position and size on every frame, and averaging. The result wasn't better from the original frames :(

How could I get the number?


texts: Impossible to make a scalogram of a signal with 262144 samples

texts:
I have a signal with 262144 samples and "sampled" with 200Msps.

I want to use the cwt (Continuous Wavelet Transform), for making a time-frequency scalogram.

Also, I worked in some scripts in python and Matlab, but it's very hard for the machine for make the process and empty the memory of the computer with 262144 samples.

The codes and examples than I run, only work with 1024 samples... 512 samples, and it's very simple for the computer to work with it. But for 262144 SAMPLES !, a 1D signal, it's impossible to work well and output a great scalogram (It's the target).

The signal has frequencies around 10KHz, 20KHz... 200KHz, 1MHz... 2MHz.

I'm really comfortable with Python language, very powerfull and great for data mining, also Matlab, but not eficient in memory use... I think...

Regards and I waiting any suggestions or code in Python, also toolboxes in Python.


texts: Implementing digital filter by padding with zeros

texts:
I am making a 9th order digital bandpass filter with lower and upper corners of 200 kHz and 40 MHz respectively. I am using this filter to filter a 1D time domain signal which is 64k samples long sampled at a frequency of 150MHz.

I have done some digital filtering before in university, so I know what to expect, but it's been a while.

I have used this site: http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
to generate the filter coefficients and the gain, and I can see how the code works, my question is this:

I start with the value of X[0], then to calculate the value of Y[0] I need values for Y[-1] ... Y[-18] and x[-1] ... X[-18].

I know these do not exist, so I think (from what I remember from university) that I pad with zeros, however, doing a bit or reading I heard it mentioned that padding with zeros changes the sampling frequency.

So how do I go about calculating this new sampling frequency (if indeed required..)?


texts: Auto Correlation vs Cross Correlation vs Convolution and their applications

texts:
I know from wikipedia that auto correlation in done on the same signal while cross correlation is done on different signals.But what does this actually imply in terms of application.I can always apply cross correlation on same signals and get same output.
And in convolution one signal is reversed.Mathematically i understaand the formulas.

But what does these three mean in terms of applications?


texts: Remove background noise from audio file (python or matlab)

texts:
I am implementing a project for infant cry detection and the audio set contains background noises. So for preprocessing i need to remove the background noise from the audio.I am not able to get a proper output for the code in jupyter notebook. For this code the output file does not contain anything even the baby cry is erased.

import numpy as np
import scipy as sp
from scipy.io.wavfile import read
from scipy.io.wavfile import write
from scipy import signal
import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline')

(Frequency, array) = read('cry.wav')

len(array)

plt.plot(array)
plt.title('Original Signal Spectrum')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

FourierTransformation = sp.fft(array)

scale = sp.linspace(0, Frequency, len(array))

plt.stem(scale[0:5000], np.abs(FourierTransformation[0:5000]), 'r')
plt.title('Signal spectrum after FFT')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

GuassianNoise = np.random.rand(len(FourierTransformation))

NewSound = GuassianNoise + array

write("New-Sound-Added-With-Guassian-Noise.wav", Frequency, NewSound)

b,a = signal.butter(5, 1000/(Frequency/2), btype='highpass')

filteredSignal = signal.lfilter(b,a,NewSound)
plt.plot(filteredSignal) # plotting the signal.
plt.title('Highpass Filter')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

c,d = signal.butter(5, 380/(Frequency/2), btype='lowpass') # ButterWorth low-filter
newFilteredSignal = signal.lfilter(c,d,filteredSignal) # Applying the filter to the signal
plt.plot(newFilteredSignal) # plotting the signal.
plt.title('Lowpass Filter')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

write("New-Filtered-Sound.wav", Frequency, newFilteredSignal)

Can i know how to debug this code or get any other code.


texts: help with bandwidth diagrams for digital transmission

texts:
I am taking a undergrad network course that uses Data Communications and Networking 5th edition containing diagrams and legends which I cannot decipher.
Specifically the diagram on the right, has y-axis P for power, then z-axis f/N.  Does this mean the fraction f/N (frequency/bit rate) or does it mean f or N.  Also the start of chapter 4 explains that ratio r is data rate/signal rate or r = N/S.
In the legend it gives r=1 so re-arranging the equation, shouldn't S_ave = N/1?

Also, it explains that Manchester and Differential Manchester take up twice as much bandwidth because of the voltage change during the bit but comparing the right hand side diagram below with the respective right hand side diagram above, how do you determine that it is double?  Is it area under the curve?


texts: Are Total Harmonic Distortion and Noise (THD+N) and THD correlated in practice?

texts:
I've read about THD and THD+N due to interest in audio quality.

Some audio devices (ej. speakers) report %THD or %THD+N. %THD+N is greater than %THD. I've found that sometimes the difference is huge, which makes me think that a significant part of the noise is due to large generation of spurious frequencies (which are not harmonics of the reference input frequency).

I think that in similar circumstances non-harmonic frequencies should generate to the hearer a perceived distortion greater than those due to harmonics frequencies (Am I right?). If this is true, the %THD+N must be much more relevant than %THD. In such case %THD can only be relevant for comparing two audio devices if it is positively correlated with %THD+N.

Q1 Are %THD+N and %THD correlated in practice in the case of common audio devices (speakers, headphones, amplifiers)?

Q2 How important is %THD in the reproduction of the input signal in speakers or headphones?


texts:
I have this lame doubt which has ended up confusing me which is probably trivial. I am modelling a MU-MIMO system with a Rayleigh fading channel. I have created a wanted signal as an OFDM signal, and generated a 3D Rayleigh Channel matrix over all sub carriers. I want to now multiply the transmit vector(OFDM signal) to the 3D channel matrix. Should i convert both the channel matrix and wanted signal in frequency domain and add noise

$y(f) = H(f)x(f)+n(f)$, and then convert it back to time domain
$\textrm{ifft}(y(f))$ and provide an estimate using MRC, ZF or MMSE methods.

$\hat{y}(t) = H(t)’ y(t)$

Or is it just a direct multiplication of

$y(t)=H(t)x(t)+n(t)$

as i understand from most of the literature (y = H x+n)is the usual system model in time domain. If so why isn't it convolved when its combined in time domain, similar to any communication system, where the input vector is convolved with the channel matrix to obtain the output vector in time.


texts: What Is the Intuition of Convolution in The Signal Processing World

texts:
It is known that an LTI system is defined by its Impulse Response and the Convolution Operator.

Could some put some intuition behind?

Thank You.


texts:
Schaum's Outline, Digital Signal Processing, Second edition, 2012, page 101:

Prove that:

$$X(e^{j\omega}) = \frac{1}{T_s}\sum_{k=-\infty}^{\infty}X(j\frac{\omega}{T_s}-j\frac{2\pi k}{T_s})$$


texts: Notch filters: why is $H(e^{j0}) = H(e^{j\pi})$ required?

texts:
I'm studying for the final exams and many example exercises require to design fir notch filters of length 5, given some passband and some stopband frequencies. You've got to find 3 coefficients α0, α1, α2. You need 3 equations.

R(passband ω) = 1
R(stopband ω) = 0

I completely understand why these are here. The final degree of freedom is determined by the equation $H(e^{j0}) = H(e^{j\pi})$, because such thing is true for notch filters. I struggle to understand why this is true.


texts:
I have a device that connect with USB to my computer. This device have digital signal that I can sample it every few seconds. The signal that the device perform, is a digital signal with $X$ seconds of $1$ and $Y$ seconds of $0$ ($X$ and $Y$ are integers and consts).

What I want to do is to estimate if in $Z$ seconds from now, the signal will be $1$. The problem is that I don't know what are $X$ and $Y$ (I don't know how much time $1$ and $0$ are), so I don't know how can I estimate it, and I can sample it every 10 seconds ($X$ and $Y$ can be less then 5).
But if I sample $0$ at some moment, the device tells me in how much seconds $1$ will appear again.

I have basic knowledge in signal analysis and I don't know how to do that.

Is it possible at all? If it is, how can I do it?


texts: Is FFT sub matrix degeneracy a problem in OFDM when there is noise in the frequency domain?

texts:
Say that the discrete Fourier transform (DFT) is used in OFDM. There are a  number of degenerate (singular, non invertible) sub matrices of some DFT matrices. Does this result in any problems?
One example from FFT is where some basis vectors have some elements in common like in this size 8 DFT/FFT. Let us express the DFT as a matrix multiplication X = M*x

Say we have an input vector with a non zero sum of x[0] and x[4] and other elements zero. Output X[0], X[2] and X[4] will have the same value regardless of how the sum is distributed between x[0] and x[4]. This does not look good. It is the same as saying that the submatrices
$$
\begin{bmatrix}M_{00}&M_{04}\\M_{20}&M_{24}\end{bmatrix} =
\begin{bmatrix}M_{20}&M_{24}\\M_{40}&M_{44}\end{bmatrix} =
\begin{bmatrix}M_{00}&M_{04}\\M_{40}&M_{44}\end{bmatrix} =
\begin{bmatrix}1&1\\1&1\end{bmatrix}
$$
are degenerate, having determinant zero. Similar sub degeneracies exist for larger $2^n$ sample sizes. The entire resulting vector X of-course contains all information in the original vector x but the meaning with some applications of OFDM is to transform a signal into a domain with noise.
X is then transmitted, and some noise is added, and an inverse FFT is performed. It has in this example the matrix formulation x = $M^{-1}$ X
$$
M^{-1} = \frac{1}{8}\begin{bmatrix}
1&1&1&1&1&1&1&1\\
1&\frac{1}{\sqrt2}(1+i)&i&\frac{1}{\sqrt2}(-1+i)&-1&\frac{1}{\sqrt2}(-1-i)&-i&\frac{1}{\sqrt2}(1-i)\\
1&i&-1&-i&1&i&-1&-i\\
1&\frac{1}{\sqrt2}(-1+i)&-i&\frac{1}{\sqrt2}(1+i)&-1&\frac{1}{\sqrt2}(1-i)&i&\frac{1}{\sqrt2}(-1-i)\\
1&-1&1&-1&1&-1&1&-1\\
1&\frac{1}{\sqrt2}(-1-i)&i&\frac{1}{\sqrt2}(1-i)&-1&\frac{1}{\sqrt2}(1+i)&-i&\frac{1}{\sqrt2}(-1+i)\\
1&-i&-1&i&1&-i&-1&i\\
1&\frac{1}{\sqrt2}(1-i)&-i&\frac{1}{\sqrt2}(-1-i)&-1&\frac{1}{\sqrt2}(-1+i)&i&\frac{1}{\sqrt2}(1+i)\\
\end{bmatrix}
$$
Some degenerate sub matrices, there are even more, are
$$
\begin{bmatrix}{M^{-1}}_{00}&{M^{-1}}_{04}\\{M^{-1}}_{20}&{M^{-1}}_{24}\end{bmatrix} =
\begin{bmatrix}{M^{-1}}_{20}&{M^{-1}}_{24}\\{M^{-1}}_{40}&{M^{-1}}_{44}\end{bmatrix} =
\begin{bmatrix}{M^{-1}}_{00}&{M^{-1}}_{04}\\{{M^{-1}}^{-1}}_{40}&{M^{-1}}_{44}\end{bmatrix} =
\begin{bmatrix}1/8&1/8\\1/8&1/8\end{bmatrix}
$$
so for noise where X[0]+X[2]+X[4] = constant, x[0] and x[4] stays the same. Doesn't this change the signal to noise ratio?


texts: Haar series approximation of a function?

texts:
I'm trying to code the Haar wavelet approximation of a function.
I'm a complete newbie.

Edit
Below is a pseudo-code of my algorithm for $f: [0,10] \to \mathbb R$, $f(x)=x$.
I picked it here and there on google.

1) Haar wave formula of order $N$  for $x \in [0, \infty)$
$$
\psi_{N,j}(x) = \frac{1}{\sqrt{N}} \cdot\begin{cases}
1 & j =0 \\
+2^{p/2} & a \le x - \left \lfloor \frac{x-a}{L}\right\rfloor L< b \\
-2^{p/2} & b  \le x - \left \lfloor \frac{x-a}{L} \right\rfloor L < c\\
0 & \text{otherwise}
\end{cases}
$$
where
\begin{align}
a &= \frac{q-1}{2^p},
b = \frac{q-0.5}{2^p},
c = \frac{q}{2^p},
L = c-a, \\
p &= \left\lfloor \log_2(j+1) \right\rfloor,
q = j +1 - p.
\end{align}

2) Haar series coefficients
$$
\hat \psi_{N,j} = \int_0^{10} f(t) \psi_{N,j}(t) dt
$$

3) Haar series
$$
f(x) \approx \sum_{j=0}^N \hat \psi_{N,j} \cdot \psi_{N,j}(x)
$$

4) For $N=4$, this gives the graph


texts: Are least square filters, or filters that minimize error energy, the same as least mean square adaptive filters?


texts: How to determine the FFT size of LTE downlink receiver?

texts:
I'm puzzled about the LTE FFT size.

Assuming an LTE base station works at 20M mode, the FFT size is 2048. Then, many users can be served by this base station.

Do all these users work at 20M mode (i.e., the FFT size of all these users are also 2048)?

OR

This 20M base station can serve users working at various modes (e.g., 1.4M, 5M, 10M, 15M, 20M modes)? For example, one user works at 1.4M mode with 128-point FFT and one user works at 10M mode with 512-point FFT can be served simultaneously by this 20M mode base station?


texts: the position of image plane in optical imaging

texts:
In the book Learning OpenCV, there is one figure that show the image plane always is focal image.

The description of the figure from the book is as follow.

We begin by looking at the simplest model of a camera, the pinhole camera model. In this simple model, light is envisioned as entering from the scene or a distant object, but only a single ray enters from any particular point. In a physical pinhole camera, this point is then “projected” onto an imaging surface. As a result, the image on this image plane(also called the projective plane) is always in focus, and the size of the image relative to the distant object is given by a single parameter of the camera: its focal length. For our idealized pinhole camera, the distance from the pinhole aperture to the screen is precisely the focal length. This is shown in Figure 11-1, where f is the focal length of the camera, Z is the distance from the camera to the object, X is the length of the object, and x is the object’s image on the imaging plane. In the figure, we can see by similar triangles that –x/f = X/Z,

From this link https://en.wikipedia.org/wiki/Lens_(optics)

There is thin lens equation. Only object in infinity can be imaged in focal image. Which one is correct?

Edit:
From my previous knowledge, camera calibration is performed for cameras using fixed focal length lenses using methods from the opencv library. In calibration process, the pinhole camera model is assumed and the image plane is f away from optical center. In the book Learning OpenCV, there is one chapter which gives detailed description of calibration process.

Nowadays I learn that the same procedure is used for cameras on mobile phones.
Such cameras have the autofocus functionality. AF is implemented by changing the position of the lens by using VCM(voice coil motor). The link have one figure which shows the process. http://electroiq.com/blog/2011/11/mems-alternatives-for-miniature-auto-focus-cameras/. The focal length of the lens is fixed. If you capture the chess board each time, the VCM will modify the lens to make objects at the center of the image in focus.
Can the calibration procedure apply to such system?


texts: How to Design a Pre-emphasis filter knowing the gain, a zero and a pole?

texts:
I am simulating an FM radio signal in order to use it as a waveform in passive radar; for this purpose I'm following an article:
Analysis and Emulation of FM Radio Signals for Passive Radar

In a certain step, I need to design an emphasis filter for the rdiophonic signal, the filter has these characteristics:

a gain of -5.4 dB ?(is it the gain at half of sampling frequency?)
a zero imposed by international standards (can't find the zero)and a pole - to make it
physically realizable - at 40 kHz. I know the response function of the Pre-emphasis filter

$$ H(j \omega) = 1 + \frac{j \omega}{\omega_0} $$

or in z domain:

$$H(z) = 1-B \cdot z $$

, can't understand where is the pole?
Of course i will use a FIR filter ,but how to find the coefficient $B$.

thanks for all


texts: Is this an error in Oppenheim and Schafer's Discrete-Time Signal Processing (3rd edt)?

texts:
Q: According to Fig 71, y[n]=something+e[n] and hence the transfer function from e[n] to y[n]
denoted as He(z) should be 1. Am I wrong?


texts: How to stitch together MFCCs from multiple frames?

texts:
I am trying to record a sound two times and compare the two recordings. I am playing exactly the same sound both times (for example, play a small .wav file twice and record it each time), so I hoped that my attempt would not rise to the complexity of Speech Recognition...

I have tried computing the DFT (using FFT) and comparing the two power spectrums of the signals (more precisely, comparing their MSE - Mean Squared Error), but the results were discouraging (sometimes for recordings of the same sound the error was greater than for different sounds)...

Then I found out about the MFCCs and computed them for each signal. That leaves me with 12 coefficients for each frame I devided my signal into (details below). How should I combine these groups of 12 coefficients to get some value(s) that describe the whole signal?

(Appending them and extracting the MSE for the two arrays did not return encouraging results either...)

Below are some details about the way I computed the MFCCs and values I chose for computing it. Maybe it's here that I'm doing something wrong:

(sample rate = 22000 Hz)

Pre-emphasis.
Framing with frame size = 512, frame overlap = 200. (because "each frame should be ~20-30ms long" and 512 / 22000 = 0,023 seconds = 23 ms. As "the overlap should be of ~10ms and 200 / 22000 = 0,009s = 9ms)
Apply Hamming window to each 512-size frame.
Apply FFT to each windowed 512-size frame => 512 magnitudes for my signal, from witch I use only the first 256. Domain: Hz from 0 to 22000 (to 11000, respectively)
Compute the Mel Filter Bank: (min frequency = 300 Hz, max frequency = 11000 Hz)

Compute mels from min and max frequences and then compute 26 equally distanced values between these two mel values.
Convert them back to Hz => array of 28 frequency filters.
Compute a filterbank (filter triangle) for each three consecutive values => 28 - 2 = 26 filterbanks.
Pass the whole power spectrum (256 magnitudes) from step 4 through each triangular filter to get a "filterbank energy" for each filter => 26-sized array of energies.

Apply log: ln(each energy) (they are still frequency-domain values)
Apply DCT to the logged energies => 26-sized array of (time-domain?) values.
Take only the first 12 => 12 MFCCs for each 215-sized frame of my signal.


texts: Calculate moving RMS with exponential averaging

texts:
This is moving RMS, now how to add exponential averaging in this
df['signal'].rolling(21*24*60).apply(lambda x: np.sqrt(x.mean()))


texts: Is it possible to reduce the complexity of radix-2 FFT if the input vector contains identical elements?

texts:
my question is about reducing the complexity of radix-2 FFT when the input vector has a specific structure.

For an input vector of x with N elements, the complexity is given by O(N log2 N).

My input vector y with N elements has the structure of:
y = [x1 x2]T

where,

x1 = x2

and x1 and x2 vectors have N/2 elements. Is it possible to reduce the complexity of radix-2 FFT for this y vector? I think it is possible to not use one of the N/2 DFT blocks of the FFT algorithm in this case. Can anyone show me what is the complexity for this vector? Thanks in advance.


texts: Implementation of Fourier Domain Denoising with Hard Threshold

texts:
I just tried the Fourier denoising method with a hard threshold and my code is as follow:

F = imnoise(phantom(128),'gaussian');
B = fft(F);
s = 0.3;
for i = 1:128
for j = 1:128
    if abs(B(i,j))<s
        B(i,j)=0;
    end
end
end

C=ifft(B);
subplot(121)
imshow(F,[])
subplot(122)
imshow(C,[])

Left is the phantom after adding noise and the right is the denoised image from my code.

The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients and the coefficients below the threshold to be 0.

In my code, I just added gaussian noise with default mean 0 and variance 0.01 to my phantom image and then I set the threshold to be 3 times of the standard deviation of the noise, which is 0.3 (this is some good practical threshold I looked at from literature). But why the result is not good? It seems it is still very noisy. Even I change my threshold the result is not good. Why? And I want to my denoised image as close to the original phantom.

Thanks in advance!


texts: Confused about the period of a discretised time-series

texts:
I have managed to confuse myself about the periods of discretised time series. It is all rather embarrassing.

This is one period extracted from my continuous signal (radians on y, time on x):

To use it in my calculations I use a discretised vector of N points. The above signal has a period of $T=0.95s$. The discretised version A has $N=100$ and the angle, over time, is indexed from 1 to 100.

Q1: is the period of my angle vector $[\theta_1,\ldots,\theta_{100}]$; $T=100$?

Now suppose I need to standardise my data as such, incl the index:

% standardise training inputs
m = mean(A);
s = std(A);
x(:,1) = 1/s(2) * (A(:,2) - m(2));  % index
x(:,2) = 1/s(1) * (A(:,1) - m(1));  % velocity

Q2: Is the new period then simply max(x(:,1))?


texts: Lossless algorithm for hyperspectral imaging

texts:
I am looking for Lossless compression algorithm for hyperspectral imaging. The only one that I found was LPVQ.
Do you know other algorithms? I am looking for pseudocode or matlab code.

Thanks


texts:
I had this question in the exam without any further explanation.

Why the time translation invariance is an important feature for time-frequency distributions?

I am writing to ask whether anyone can please explain what is time translation invariance? And why is it important feature?


texts: MATLAB, error while trying to implement a Designed Filter object

texts:
I am using the Digital Filter with properties in MATLAB. I have defined the filter object filt(1) as:

 digitalFilter with properties:
            Coefficients: [20x6 double]
   Specifications:
       FrequencyResponse: 'bandpass'
         ImpulseResponse: 'iir'
              SampleRate: 1000
      PassbandFrequency1: 59.9000
          PassbandRipple: 3
    StopbandAttenuation2: 80
    StopbandAttenuation1: 80
      PassbandFrequency2: 60.0800
             FilterOrder: 40
            DesignMethod: 'ellip'

Similarly, I have other filters under filt object as filt(2), filt(3) up to filt(8), each at different frequencies.

Now, for some reason, I get the error when my program tries to run the line:

Q = filter(filt(1),x);

where x is a 600000 x 1 column vector with only positive and negative real values (resembling a long sine wave).

The error that I get is:

Subscript indices must either be real positive integers or logicals.

Error in digitalFilter/filter (line 870)
    y = obj.FilterFunction(obj.Coefficients,x);

Could anyone please point out the problem?

EDIT:

The designfilt command opens up an interactive interface where you can change and tune your filter requirements. I have optimized the filters, and saved the results as a .mat file. (This whole program used to work, was part of a significant project which we completed succesfully, but after a year doesn't seem to work somehow).

I have attached the filter_8th.mat file. There are various filters here. I have already created these filters before hand because MATLAB takes a lot of time to create these.

In this .mat file, the filter objects are named as D_n1_50, D_n1_60, D_sig50, D_sig60 etc. (easily recognizable from the Workspace.)

The line of code that assigns these filter to the variable filt is, for example: filt = D_sig50. This is being selected in a loop. Which filter will be selected depends on a conditional if-else structure.

The .mat file containing "filters" and the variable x are in this dropbox link:


texts: correlation matrix vs. correlation function?

texts:
Can someone help me understand the difference between a 1-dim autocorrelation function and a 2-dim autocorrelation matrix of a random process aka time series?
My Leon-Garcia textbook defines CX(τ) and RX(τ) as 1-dim functions, yet Haykin's Adaptive Filter Theory defines what looks like the same thing, but as the correlation matrix R instead. I see that the matrix appears to arise when the time series is represented by complex numbers (sinusoids), but is that the only real difference, or is there a fundamental distinction that I'm missing?
And does the correlation matrix R reduce to the 1-dim autocorrelation function when the time series can be represented by only real numbers instead of complex numbers, i.e. in that case the correlation matrix is unnecessary?
thanks,
js.


texts: Why prediction by AR model causes time lag

texts:
Why does prediction by AR model cause a time lag?

Please tell me why theoretically.


texts:
I am working on DSP and finding difficulty to understand the term Normalized frequency often used with DFT & DTFT.

What is normalized frequency in DSP ? and how it is different from analog frequency ?

What is significance to normalize frequency in DSP ?

Why limit of normalized frequency is 2π ?

How FFT deals with normalized frequency ?


texts: Extracting time-tagged pulses from Geiger audio recording

texts:
I am doing some tests with SmartGeiger detector designed for smartphones. The sensor unit sends pulses via phone microphone jack input and the app does the counting and calculates the values.

However, I need to record and process a lot of data and the app is not comfortable for this task (you have to save every value manually). I recorded the data with Audacity - looks this way:

but I need to somehow convert it to the pulses - e.g. pick up the pulse data using some simple threshold filter and count the number of pulses for particular time interval.

One pulse looks this way in the audio data:

Unfortunately the Audacity "Sample Data Export" feature:
https://manual.audacityteam.org/man/sample_data_export.html

seems to only work for short parts of the data. Is there any way to do the conversion?
I only need the temporal location of the pulses.

thanks


texts: Create a color image in Matlab with sinusoidal patterns

texts:
I want to create a color image in Matlab whose components have sinusoidal patterns (okay even if they have different frequencies).

To create a black and white image I did:

fs = 0.08;
W = 256;
H = 128;

img1 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img1(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

colormap(gray(256));
imshow(img1);
title('img1');


texts: How to create/find a dereverberation filter?

texts:
I want to find the coefficients of dereverberation filter in MATLAB. How to do that?
I have found the coefficients of reverberation filter:
numerator  = b = [0.1664 0 0.4084 0 0.3341 0 0.0911]
denominator = a = [1 0 0 0 0 0 0]

Thanks in advance!


texts: A question about the meaning of pole in time domain

texts:
Lets say I have a transfer function $H(s)$ of a system defined in $s$-domain as:
$$H(s) = \frac{1}{s - (-1-j)}$$

So I conclude that the pole on the $s$-plane is where $s = 1+j$. So far so good.

Now does that mean if the Laplace transform of the input to the system is $s = 1+j$ the system goes crazy/infinity/oscillate?
Does that mean if I take the inverse Laplace transform of $1+j$ I can find what input in time domain makes the system unstable?


texts: Find cut-off frequency of a low-pass filter for a given output signal

texts:
Suppose we have an ideal low-pass filter $H(e^{j\theta})$ with a cut-off frequency $\theta_c$ in the range of $0\leq \theta_c \leq \pi$.

I want to know the input signal $x[n]$, as well as the cut-off frequency $\theta_c$, which produces the following output signal $y[n]$.

To solve this I started off with transforming $y[n]$ to the frequency domain by calculating the DTFT:

$$
Y(e^{j\theta})=\sum_{n=-\infty}^{\infty}y[n]e^{-j\theta n}=1+e^{-j\theta}+2e^{-j2\theta}+e^{-3j\theta}+e^{-j4\theta}
$$

Now I can calculate $X(e^{j\theta})=\frac{Y(e^{j\theta})}{H(e^{j\theta})}$. This is only possible for $|\theta| \leq \theta_c$, since outside of this region $H(e^{j\theta}) = 0$ and $X(e^{j\theta})$ would approach infinity. But inside this region $H(e^{j\theta})$ is simply 1 and when transforming $X(e^{j\theta})$ back to $x[n]$ it is just the same as $y[n]$, a behaviour that is expected of a low-pass.

The thing I am not sure about is the cut-off frequency. If I simply plot $Y(e^{j\theta})$, I see that there are frequency components until $\pi$. So I could just say alright, $\theta_c=\pi$. But is it really like that? Because the low-pass is actually just dependent on $\theta$ and not all the multiples of it, which I find in the complex exponentials of $Y(e^{j\theta})$. So how do I calculate $\theta_c$?


texts: Bandwidth of a finite impulse response filter

texts:
In a machine learning application, a model learns to apply a filter h(n) via convolution to a 1-dimensional input signal s(n) (e.g. sampled at Fs = 16 kHz). The filter size (filter length?) is limited to a fixed number, e.g. k = 100, so the model cannot learn a filter that exceeds this duration but it is free to insert zeros and make the filter "shorter".

Does this mean that the filter size limits the duration of the "slowest" signal component that falls into the passing bandwidth and hence Fs/k = 160 Hz is the minimal bandwidth that h(n) can represent (assuming it's a band pass filter)?


texts: How to distinguish between FFT power leak or close frequencies

texts:
This week I was having an interview with a Professor for Grad School. He asked a bunch of questions related to signal processing that lead up to some questions about power leak on the FFT spectrum. So far so good, but then he caught me when he asked something around those lines:
"given just an FFT spectrum, with no prior information about the signal or the sampling frequency, how one can differentiate if two bins close together are a case of power leakage or two distinct frequencies?"
At the time of the interview nothing came to mind so I just answered that it was not possible but since then I was wondering if there is a way of making this distinction. I tried looking online but couldn't find any material that answered this so any help would be greatly appreciated.
Thanks!


texts: Help calculating / understanding the MFCCs: Mel-Frequency Cepstrum Coefficients

texts:
I've been reading bits and pieces online but I just can't piece it all together. I have some background knowledge of signals / DSP stuff which should be enough prerequisites for this. I'm interested in eventually coding this algorithm in Java but I don't understand it fully yet which is why I'm here (it counts as math, right?).

Here's how I think it works along with the gaps in my knowledge.

Start with your audio speech sample, say a .wav file, that you can read into an array. Call this array $x[n]$, where $n$ ranges from $0, 1, \ldots ,N-1$ (so $N$ samples). The values correspond to audio intensity I guess - amplitudes.
Split the audio signal into distinct "frames" of 10ms or so where you assume the speech signal is "stationary". This is a form of quantization. So if your sample rate is 44.1KHz, 10ms is equal to 441 samples, or values of $x[n]$.
Do a Fourier transform (FFT for computation's sake). Now is this done on the entire signal or on each separate frame of $x[n]$? I think there's a difference because in general the Fourier transform looks at all elements of a signal, so $\mathcal F(x[n]) \neq \mathcal F(x_1[n])$ joined with $\mathcal F(x_2[n])$ joined with $\ldots \mathcal F(x_N[n])$ where $x_i[n]$ are the smaller frames. Anyways, say we do some FFT and end up with $X[k]$ for the rest of this.
Mapping to the Mel scale, and logging. I know how to convert regular frequency numbers to the Mel scale. For each $k$ of $X[k]$ (the "x-axis" if you'll allow me), you can do the formula here: http://en.wikipedia.org/wiki/Mel_scale. But how about the "y-values" or the amplitudes of $X[k]$? Do they just remain the same values but shifted to the appropriate spots on the new Mel (x-) axis? I saw in some paper there was something about logging the actual values of $X[k]$ because then if $X[k] = A[k]*B[k]$ where one of those signals is presumed to be noise you don't want, the log operation on this equation turns the multiplicative noise into additive noise, which hopefully can be filtered (?).
Now the final step is to take a DCT of your modified $X[k]$ from above (however it ended up getting modified). Then you take the amplitudes of this final result and those are your MFCCs. I read something about throwing away high frequency values.

So I'm trying to really iron out how to compute these guys step by step, and clearly some things are eluding me from above.

Also, I've heard about using "filter banks" (an array of band pass filters basically) and don't know if this refers to making the frames from the original signal, or maybe you make the frames after the FFT?

Lastly, there's something I saw about MFCCs having 13 coefficients?


texts:
I understand an ideal analog square wave contains sine-wave-components above Nyquist frequency (and towards infinity) and is thus subject to aliasing. So obviously we need an anti-aliasing-filter. Or do we? Why do we need one?

Intuition tells us that a square wave can hardly be aliased. It is so simple and there's no chance for a cycle or half cycle to be "skipped" unless the base frequency was above Nyquist... yet listening to a non low-passed-square-wave after sampling tells us otherwise. Many low pitches are heard, because of aliasing.

How is this possible?

Each successive sine wave seems to counteract the previous one, reducing the ringing at the nulls of the square-wave

Why does superposition hold on such a simple construct?


texts: Filtering random noise from a signal using excel

texts:
I was given an excel file, and I was asked to filter the noise so that signal would be as close as possible to the original one, using also excel.

I kinda understand that I need to do $(x(n)+z)-z$ ($z$ being the noise), but I don't understand how to find the exact formula to find the noise when the noise itself is generated using rand(). We were given option to use DFT FFT FIR IIR, so I guess the solution would be that, but my lecturer didn't really teach us about FIR and IIR (due to corona I guess, not enough time) and the lecture never really cover the part of using DFT FFT to signal with noise. So, how to solve this problem using excel?

Signal given (without noise) is as below:
$$x(n) = \sin\left(2\pi \frac{8}{128}n\right)$$
Information Frequency: 8Hz

Sampling Frequency: 128Hz


texts: Convolution in frequency space with fft and ifft -> num output samples?

texts:
Let's say you had two 1 second sound samples at 44.1khz, and you wanted to convolve them in frequency space.

You FFT each one to get two arrays of complex numbers that each have 44100 items in them.

You then multiply each bin...
for i = 0 to 44099
c[i] = A[i] * B[i]

Then it's time to turn it back into time domain samples.  If you do an ifft you'll get 44100 time samples out, but when you convolve these two audio samples in time space, it'll give you 88199 samples out (44100+44100-1).

Obviously getting 44100 samples out is incorrect.  At what part in the process did I go wrong? (:

Thank you!!


texts: How to convert a spectrogram back to a signal

texts:
I'm converting a signal to a spectrogram, manipulating that (nonlinear stuff), and then want to use the modified audio signal. I can make it so the segments are not overlapping, but am unsure how to deal with the complex part of the signal.

The manipulation I'm doing is increasing the value of the smaller values (but not modifying large ones), effectively lowering SNR. I can't think of how I would do that in the time domain. And I'm unsure whether it would make sense to just do the same thing to the complex values? And then I guess I would have to do ISTFT for each segment, I don't see a function to do that for several segments (a spectrogram) in scipy so I guess I would have to do that manually? Anyone else had to do something like this before?


texts: Temporally locating beats in an audio file

texts:
Imagine you have a long recording of a guitar player playing an arpeggio at a slow speed. To simplify, assume the tempo is 60 bpm or lower. I want to annotate this recording with the exact moments when the guitar player plucks the strings. So:

$$
\mathrm{detect}(recording.wav) = \{t_1, t_2, t_3, \ldots\},
$$
where $t_1, t_2, \ldots$ are the millisecond offsets in the recording of when the strings are plucked.

I have skimmed similar questions related to beat and note detection but they are not really what I'm looking for. Beat detection doesn't work because I'm trying to detect when notes are played not beats. Note detection isn't right either because I'm only looking for when the note is played, not the frequency of it. Probably I'm missing something because this problem seem so trivial.

You can find examples of the recordings I'm talking about here:

https://www.dropbox.com/s/3852tnat73nld2w/beats.mp3?dl=0
https://www.dropbox.com/s/tr986i1uoh11i9g/beats2.mp3?dl=0
https://www.dropbox.com/s/ncosz1d56n8jjst/noisy-beats.mp3?dl=0

The file-format happens to be mp3 but could just aswell be wav.

After some research, I've found that what I'm looking for is called "onset detection." I would have hoped that there would exist robust methods for detecting when the notes are played, since it is so easy for human ears to hear it. But perhaps that is not the case.


texts: Difference between real fft and complex fft with imaginary part of zero in fftw?

texts:
I have a real 2d matrix. I am taking its fft using fftw. But the result of using a real to complex fft is different from a complex ( with imaginary part equal to zero) to complex fft.

real matrix

 0     1     2
 3     4     5
 6     7     8

result of real to complex fft

36 -4.5+2.59808i  -13.5+7.79423i
0  -13.5-7.79423i 0
0  0              0

Code:

int r = 3, c = 3;
int sz = r * c;
double *in = (double*) malloc(sizeof(double) * sz);
fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
fftw_plan p = fftw_plan_dft_r2c_2d(r, c, in, out, FFTW_MEASURE);
for ( int i=0; i<r; ++i ){
    for ( int j=0; j<c; ++j ){
        in[i*c+j] = i*c + j;
    }
}
fftw_execute(p);

using a complex matrix with imaginary part of zero

complex matrix

 0+0i     1+0i     2+0i
 3+0i     4+0i     5+0i
 6+0i     7+0i     8+0i

result of complex to complex fft

36               -4.5 + 2.59808i  -4.5 - 2.59808i
-13.5 + 7.79423i 0               0
-13.5 - 7.79423i 0               0

Code:

int r = 3, c = 3;
int sz = r * c;
fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
fftw_complex *inc = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
p = fftw_plan_dft_2d( r,c, inc, out, FFTW_FORWARD,FFTW_MEASURE);
for ( int i=0; i<r; ++i ){
    for ( int j=0; j<c; ++j ){
        inc[i*c+j][0] = i*c+j;
        inc[i*c+j][1] = 0;
    }
}
fftw_execute(p);

I am after the result of complex to complex fft. But the real to complex fft is much faster and my data is real. Am I making a programming mistake or the result should be different?


texts: Development of Hilbert transform relationship

texts:
In the development of Hilbert transform relationships, Prof. Oppenheim has chosen
\begin{equation}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta
=X_R\left(e^{j\omega}\right)
\end{equation}
in his book "Discrete-Time Signal Processing, 2e" Chapter 11 pp807, Eq. 11-24. I failed to understand how he arrives at it.

Does anybody any idea how has derived, please let me know?

I tried some BUT IT DID NOT GIVE the above relationship:

1:

\begin{align*}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta&=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)2\pi\sum_{k=-\infty}^{\infty}
\delta(\omega-\theta-2\pi{}k)d\theta\\&
=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)e^{j\theta{}n}d\theta\\&
=x_e(n)
\end{align*}

OR

2:

\begin{align*}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta&=\sum_{k=-\infty}^{\infty}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\delta(\omega-\theta-2\pi{}k)d\theta\\&
=\sum_{k=-\infty}^{\infty}X_R\left(e^{j(\omega-2\pi{}k)}\right)
\end{align*}


texts:
I have done fews researches and found out approach such as

multiscale and
region growing,

that could be used of extracting forearm veins.

However, I tried to use frangi filter and Figure 1 is the segmentation result using the best parameter set I could found;

Figure 1

Could anyone suggest what approaches as well the technique to adopt while considering this problem?

I give extra image for further analysis:


texts: How did letter d.) become letter e.)? How do we perform this averaging process? Its clear the averaging the signal causes some noise to be reduced.


texts: [Python ]How can I improve my 1D Kalman Filter estimate?

texts:
I have written the following code to smooth an (almost) linear function:

import numpy as np
import matplotlib.pyplot as plt

# intial parameters
n = 500
sz = (n) # size of array
t = np.linspace(0,10,n)
x = t + 2*t**2 # truth function
z = x + np.random.normal(0,7.5,size=sz) # noise

Q = 1e-4 # process variance

# allocate space for arrays
xhat=np.zeros(sz)      # a posteri estimate of x
P=np.zeros(sz)         # a posteri error estimate
xhatminus=np.zeros(sz) # a priori estimate of x
Pminus=np.zeros(sz)    # a priori error estimate
K=np.zeros(sz)         # Kalman gain

R = 0.1**2 # estimate of measurement variance

xhat[0] = 0.0
P[0] = 1.0

for i in range(1,n):

  xhatminus[i] = xhat[i-1]
  Pminus[i] = P[i-1]+Q

  K[i] = Pminus[i]/( Pminus[i]+R )
  xhat[i] = xhatminus[i]+K[i]*(z[i]-xhatminus[i])
  P[i] = (1-K[i])*Pminus[i]

plt.figure()
plt.plot(t, z,'k+',label='noisy measurements')
plt.plot(t, xhat,'b-',label='a posteri estimate')
plt.plot(t,x,color='g',label='truth value')
plt.legend()
plt.xlabel('t')
plt.ylabel('f(t)')

plt.show()

I have chosen x[0] = 0 because I know that's a good approx and P[0] = 1 because I didn't know better.

Plot of the estimate, noise and the true function:

I know that Kalman filters are not meant for non-linear functions like mine!

I have played around with both Q and R to try to find a better estimation of the true function.My function converges very slowly.  Is there a way to find a "better" Q and R?


texts: Interpolating/Decimating CIC Filter Group Delay

texts:
I have an Interpolating Cascaded Integrator Comb (CIC) filter designed as so:

Rate Change = 100
Number of Stages = 4
Differential Delay = 1

I'd like to calculate the group delay of this filter but I am not sure how.  This question on DSP Related is helpful, but I am not sure my filter will have 0 delay.

Also, Richard Lyons explanation is helpful, but he doesn't directly deal with group delay.

I also have a Decimating CIC filter designed as so:

Rate Change = 25
Number of Stages = 4
Differential Delay = 1

How would I calculate its group delay?


texts: On the transmit weight vector in MIMO-MRC systems

texts:
Recently, I am reading paper [1]. In this paper, the author wrote:

  In MIMO-MRC systems in the absence of interferers, the signal vector at the received $n_R$ antennas is given by
  $$\mathbf{r} = \sqrt{P_0} \mathbf{H}_0 \mathbf{w}_T s_0 + \mathbf{n}, \tag{1}$$
  and the decision variable for detecting the transmitted symbol $s_0$ is obtained by combining with a weight vector $\mathbf{w}_R$. The combined signal at the receiver is given by
  $$ z = \mathbf{w}^H_R \mathbf{r} = \sqrt{P_0} \mathbf{w}^H_R \mathbf{H}_0 \mathbf{w}_T s_0 + \mathbf{w}^H_R \mathbf{n}. \tag{2}$$
  To maximize the signal component, the principle of maximum ratio transmission (MRT) can be used, then the transmit weight vector $\mathbf{w}_T$ can be chosen as $$ \mathbf{w}_T = \frac{\mathbf{H}_0^H \mathbf{w}_R}{\| \mathbf{H}_0^H \mathbf{w}_R \|}. \tag{3} $$

My question is: Why the transmit weight vector can be chosen as
$$ \mathbf{w}_T = \frac{\mathbf{H}_0^H \mathbf{w}_R}{\| \mathbf{H}_0^H \mathbf{w}_R \|}? \tag{4}$$
Thank you very much.

Note

Vectors are expressed in boldface lower case letters, e.g., $\mathbf{a}$, and matrices are boldface capital letters, e.g., $\mathbf{A}$.
$P_0$ is the average received power.
$\mathbf{H}_0$ is the $n_R \times n_T$ channel matrix.
$n_T$ is the number of transmit antennas.
$\mathbf{w}_T$ is the beamforming vector at the transmitter with unit
norm, i.e., $\| \mathbf{w}_T \|^2 = 1$.
The additive thermal noise $\mathbf{n}$ is a complex Gaussian vector with mean vector $\mathbf{0}$ and covariance matrix $\sigma_n^2 \mathbf{I}$.
We assume that the noise power per antenna is identical.
The matrix Hermitian is denoted by $\mathbf{A}^H$.
The Euclidean norm of vector $\mathbf{a}$ is denoted by $ \| \mathbf{a} \| $.

Reference

[1] K. S. Ahn, "Performance analysis of MIMO-MRC system in the presence of multiple interferers and noise over Rayleigh fading channels," in IEEE Transactions on Wireless Communications, vol. 8, no. 7, pp. 3727-3735, July 2009.
doi: 10.1109/TWC.2009.080940


texts:
I am looking for methods/algorithms for signal spectrum shaping. Let's say we have a white Gaussian noise and want to have a desired power-spectral density. The first solution that comes to mind is filtering. I was wondering if there are any other techinques.
I appreciate any hints.


texts: How can I extract a sinusoid from a very noisy signal for computing the phase shift?

texts:
I have done some measurements and I got two very noisy signals. I want to extract a certain sinusoid from the signals and I want to compute the phase shift between them.

I made some simulations and I saw that the noise modifies a lot the value of the phase shift.

How can I get a trustful value for the phase shift?


texts: Create a third octave spectrum from a time signal

texts:
I have to create a third-octave spectrum from a time signal on Octave GNU. I found some code on the net to help me but I don't have all the parts of the algorithm.

I have a .csv file which contains a simple sinus temporal signal. The specifications are:

Duration: 0 to 180s
Frequency: 32768

I calculate the octave bands from 10Hz to 10kHz :

    fMin=10;
    fMax=10000;
    octs=log2(fMax/fMin);
    bmax=ceil(octs/bw);

    %Octave bands.
    fc=fMin*2.^((0:bmax)*bw);    %Center frequencies.
    fl=fc*2^(-bw/2);             %Lower frequencies.
    fu=fc*2^(+bw/2);             %Upper frequencies

After I can display all the third octave filters using Butterworth filters :

    bw = 1/3;
    numBands = length(fc);

    b = cell(numBands,1);
    a = cell(numBands,1);

    figure
    for nn=1:length(fc)
      [b{nn}, a{nn}] = butter(2,[fl(nn) fu(nn)]/(fs/2));
      [h,f] = freqz(b{nn}, a{nn}, 1024, fs);

      hold on;
      plot(f, 20*log10(abs(h)));
    end
    set(gca, 'XScale', 'log')
    ylim([-50 0])

Which gives to me :

The last filters are great but the first ones have missing points.
I don't know how to fix that and, mainly, how can I use this to make my third octave spectrum.
Does someone knows how to do this ?


texts:
In image processing we may use derivatives to help us detect the edges. While at mean time, this operation would also make the data noisier. But why do we have this effect?

My intuition is that if we denote the image data as

$I_{actual}(x, y) = I_{ideal}(x, y) + I_{noise}(x, y)$.

Then apparently when we do $I_{actual}(x_2, y) - I_{actual}(x_1, y)$, the noises will cancel out with each other? But how come they end up adding up to each other?


texts: Comparison between watershed and grabcut

texts:
I know that watershed and grabcut algorithms are both common tools for image segmentation. They both seem to rely on predefined markers to form something like a source-sink system and then uses the gradient of the image to find the boundaries.

The largest difference I observe is that grabcut uses a graph structure where the gradient translates to edge weights, while watershed intuitively works by slowly adding water to the gradient image. But to me, it seems like they would yield similar results?

My question is what are the relative strengths and weaknesses between these two algorithms? When should I use one over another?


texts: ARMA (Auto Regressive Mean Average) Process Representation as AR (Auto Regressive)

texts:
Let's say we have an ARMA (Auto Regressive Moving Average Model) process where the transfer function is a minimum phase system (Namely Invertible).

By Wold's Decomposition it is guaranteed to have MA representation of that system (Infinite order).

My question is could it be also represented by AR (Auto Regressive) system (Finite or infinite order)? If so, how could the coefficients be calculated?


texts: What is the theoretical probability of error for OOK transmission?

texts:
I am trying to simulate an optical wireless communication channel which uses OOK modulation.

Looking at equation 14 in this literature, I found that:

For electrical SNR at the receiver

$$\textrm{BER}=Q\left(\sqrt{\rm SNR}\right)$$

where (according to this page)

$$Q(z)=\frac{1}{2} \text{erfc}\left(\frac{z}{\sqrt{2}}\right)$$

Please note that:

  electrical signal-to-noise ratio is proportional to the square
  of the received optical power Pr.

Next, I implemented $Q(z)$ as follows:

getOokBer (double snr){
    double ber = 0.5 * erfc (std::sqrt (snr/2));
    return ber;
}

At this point I got stuck. I think the next step is to calculate $P_\rm e$ and/or $P_{\rm err}$ which has been discussed here too.

In this paper $P_\rm e$ is given as:

  $P_{\rm e}(k, l)$ -- function that represents an upper bound on the probability that an error is present in the chunk of bits located in interval $l$ for packet $k$.

I found the following implementation for $\rm BER$ calculations which takes $\rm SNR$ in mW. It is not clear to me whether the SNR given in the first equation above is for log scale or linear scale.

All in all, I would like to know how to calculate $P_\rm e$ and $P_{\rm err}$ for OOK?


texts: How can I apply a Gabor filter to a sine waveform?

texts:
I am an ECE student, and I am doing my project on underwater communication.

The main concept of my project is to denoise the signal by using UFB and WPD algorithms, and giving that output to a matched filter. I have already generated the sinewave and added noise and given it to the matched filter.

How do I write the code for UFB (uniform filter bank)?


texts:
I have a bloch simulation that I have implemented in Python. The implementation is very simple as:

import numpy as np

GYRO = 267538030.3797  # in radians / s / T
T1 = 0.3
inv_T1 = 1.0/0.3
T2 = 0.05
inv_T2 = 1.0/0.05
M0 = 1.0

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B) - [M[0] * inv_T2,
                                    M[1] * inv_T2,
                                    inv_T1 * (M[2] - M0)]

Now, during the excitation, the effective field is given by:
So, the pseudocode is:

\begin{align}
B[0] &= \textrm{rf}(t) \cos(2 \pi \omega_{\rm eff}  \ t) \\
B[1] &= \textrm{rf}(t) \sin(2 \pi \omega_{\rm eff}  \ t) \\
B[2] &= 0
\end{align}

Assuming, that the RF excitation is on-resomance i.e. $\omega_{\rm eff} = 0$, I track the evolution of the magnetization of the spin across time and it works and after relaxing for T1 seconds, we have this curve:

So, the transverse magnetization starts to disappear and longitudinal magnetization recovers.

Now, I wanted to simulate a readout gradient which is applied just after the RF pulse. Now, the $X$ and $Y$ values for the effective magnetization is 0 as the RF pulse is turned off. The $Z$ component of the effective field at every time step is given by the product of gradient amplitude and the position. So,

\begin{align}
B[0] &= B[1] = 0 \\
B[2] &= g * p
\end{align}

Now, when I use the same equation to solve this, I get nonsensical results, which is shown by the curve below:

As you can see, the numbers are off the charts and I was expecting an echo at the end. The gradient is a simple trapezoidal gradient and I have checked the units and it all makes sense. I was wondering if my expression for the effective field in the presence of the gradient is correct? I have seen these being represented as rotation matrices as well and I wonder if there is something like that I might have missed.

Update

As per suggestion from M529, I did the following experiments:

# Change the method to ignore relaxation effects
import numpy as np

GYRO = 267538030.3797  # in radians / s / T

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B)

Now I call the method with the different flip angles (from 0-1200) and note the vector length and it does not maintain the length at larger flip angles. So, perhaps there is a numeric issue. I increased the RF bandwidth for computation reasons. Here is the graph of the magnetization at different flips. It is a weird non-linear effect. It first goes down and then climbs back up uncontrollably.


texts: The effects of frequency resolution in FFT on the frequency spectrum

texts:
This is a section for the parameter study of my thesis (basically: why and how I did what I did in the experiments). Currently, I have to choose the fit bandwidth and FFT lines combination for the experiment (for the parameter study, it is just giving an impulse for the generator, or hitting the metal specimen with a hammer). Investigations (actual measurement) say that with frequency resolution under 0.1 Hz, the spectrum obtained would be very squishy. For example, right below is 12800 FFT lines (or 12800 FFT bins) with different bandwidths (0.5, 1, 2, 5, 10 and 20 kHz) and different frequency resolutions (39.1, 78.1, 156.3, 390.6, 781.3 and 1562.5 mHz respectively)

The spectrum with 39.1 and 78.1 mHz frequency resolution (blue and orange) are very "squishy"...
Question 1: Why do we have such squishy spectrum?
For reference: this is non-squishy spectrum

Question 2: At low-frequency range (<50 Hz), smaller FFT bin shows peak better (more data point on the same range of frequency, thus easier to spot the peaks). But why at high-frequency, smaller FFT bins lead to lower peaks (at 300 and 330 Hz)?
Question 3: Is there a "critical point" for frequency resolution? Or the trade-off between frequency resolution vs spectral/time resolution? If yes, is there any source you can point towards for me, or is this more like a rule of thumb?
Thank you for your support
(This question is also posted in Electrical engineering, and per the first comment, it is put here)


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

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


texts: Face temperature detection without thermal camera ( by cellphone for Covid19)

texts:
i like to know what method exist in image processing to have accurate temperature estimation of human body ( especially human face ) without tomography or thermal camera, so i asked about it at Physic SE site, but because of being none clear the range of digital image pressing of temperature estimation which could be used for this propose and the technical details which it could need, i asked the software part of that question here.
SO what  image processing technique do you  know which could be applied for this propose?
Thanks.


texts:
Is there an algorithm which can give the frequency and phase of a three phase signal (three phase voltages) without using zero crossing detection.

I mean an algorithm which is done in open loop.


texts: Wavelet computation with filter bank - differing results

texts:
I'm trying to get a grip on Wavelets. I've read "Wavelets, Their Friends, and What They Can Do for You" which lead me to an implementation of the discrete DWT with filter banks. Basically, I'm using the equations

$$
a_{j,k} = \left[\tilde{H} * a_{j+1}\right](2k)\\
d_{j,k} = \left[\tilde{G} * a_{j+1}\right](2k)
$$

which involve convolution and downsampling to get the coefficients $a$ and $d$. $\tilde{H}$ and $\tilde{G}$ refer to the flipped impulse response of the wavelet, for which I currently use normalized Daubechies 2 (db2), i. e. $H = \{0.483, 0.837, 0.224, -0.129\}$ and $\tilde{H} = \{-0.129, 0.224, 0.837, 0.483\}$. Convolution is defined as

$$
\left[\tilde{H} * a_{j+1}\right](2k) := \sum_{n=0}^{3} \tilde{h}_na_{j+1,2k-n}
$$

and uses a periodically extended signal. For $\tilde{G}$ it is defined likewise.

Now I wanted to transform the very simple signal $x = [1,2,3,4] = a_{1}$ to the coefficients $a$ and $d$. For testing purposes I used PyWavelets, which claims to be tested against Matlab. The result is

$$
a_{0,0}=2.82842712\ a_{0,1}=4.24264069\\ d_{0,0}=-0.51763809\ d_{0,1}=1.93185165
$$

for both implementations, i.e. my own and PyWavelets. However, for a slightly extended signal $x = [1,2,3,4,5]$, the ordering and some values of the coefficients does not match. E.g. for $a$ my implementation yields $[3.725, 5.787, 3.311]$, whereas PyWavelets yields $[3.311, 3.725,  7.106]$.

Now I'm confused and I'm trying to clarify some things, so

Is there a way to test the correctness of the decomposition? E.g. by applying inverse DWT and testing whether the resulting signal matches the original signal?
If 1. succeeds in my case, then is the ordering of the coefficients not relevant if the correct reverse ordering is applied?
Why does the book names Daubechies Wavelet as DB2, although it has 4 taps? PyWavelets does this likewise, but e.g. Wikipedia does this not, anyway, I'm not a fan of wikipedia.

Update

Thanks to datageist for answering question #3.

Regarding the other two, I came up with the idea that my convolution might be wrong, and I tested two other "modes" of PyWavelets, namely zero-padding (zpd) and periodic-padding (ppd). Before that I used the periodic mode (per) which basically also performs ppd but yields less elements, and I had problems with that. So, I changed the signal with the filter and from now on moved the filter through the signal. For ppd I get the following scheme:

$$
\begin{matrix}
x & 1 & 2 & 3 & 4 & 5 & 6 \\ \\
a_3 & 0.224 & -0.129 & 0 & 0 & 0.483 & 0.832 \\
a_2 & 0.483 & 0.832 & 0.224 & -0.129 & 0 & 0 \\
a_1 & 0 & 0 & 0.483 & 0.832 & 0.224 & -0.129 \\
a_0 & 0.224 & -0.129 & 0 & 0 & 0.483 & 0.832
\end{matrix}
$$

Multiplying row elements with $x$ and summing up yields the correct result for each row, and it also matches with PyWavelets (except for rounding issues). The same goes for zpd:

$$
\begin{matrix}
x & 1 & 2 & 3 & 4 & 5 & 6 \\ \\
a_3 & 0.224 & -0.129 & 0 & 0 & 0 & 0 \\
a_2 & 0.483 & 0.832 & 0.224 & -0.129 & 0 & 0 \\
a_1 & 0 & 0 & 0.483 & 0.832 & 0.224 & -0.129 \\
a_0 & 0 & 0 & 0 & 0 & 0.483 & 0.832
\end{matrix}
$$

which also yields the correct results.

However, for per, the results are still wrong, so I'll have to find out what they do there.


texts:
I wish that all those people on that flight were with us, and maybe they still are, as I like to think they are on an island somewhere waiting to us to find them ... however, it does not look good...

Recently, a spy plane searching for possible remnants of Malaysian Flight 370 took this photo:

However they cannot tell what this is an image of, floating in the sea.

I was wondering, what image processing techniques might us in the DSP community utilize here to better ascertain / clean up the image so that we might be able to get a better picture?

The waves of the ocean would seem to contribute to high-frequency components, so perhaps nulling the high frequency fourier co-efficients in the 2D-FFT, and then inverse transforming?

What other techniques might we use? ... perhaps we can segment it while respecting shapes of big features?...


texts:
Suppose I have an instrument which is out of tune. For example lets say, as an example, the notes it plays are as follows: C+20 cents, C#+20 cents, D+20 cents D#+20cents etc. In other words, it plays all its notes 20 cents too high.

I then record this instrument, having played polyphonic music on it.

What algorithm can I use to tune the track so that the notes become: C+0 cents, C#+0 cents, D+0 cents, D#+0 cents. I wish to solve this in the frequency domain and I have implemented an FFT and IFFT in this regard.

UPDATE

Please note that the algorithm should work out the 20 cent offset by itself without me having to supply it. I might later record another instrument with a -35 cent offset. All I know is that the notes are relatively in tune, but the instrument itself relative to other instruments might be out of tune by x cents.

I am interested in human voice, but also other instruments like violin, flute, guitar, piano and analog synth.


texts: Impulse property - needed mathematical proof

texts:
Prove:
  $$ \left. \int_{-\infty}^{\infty}\frac{d\delta(t)}{dt}\phi(t)dt=\frac{-d\phi(t)}{dt}\right|_{t=0} $$

$\delta(t)$ is impulse signal and other one is any general signal.

Can anybody do the proof for me?
Or at-least how to proceed further?


texts: What happens to the audio of a reconstructed signal if we remove too much frequency?

texts: Assuming we have an audio with noise embedded in it, I was wondering what happens theoretically if we remove (filter out) too much frequency when we pass LPF with a cutoff frequency lower than the perfect cutoff frequency needed to reconstruct the signal perfectly? Practically after performing the necessary computations on MATLAB and running the audio, the sound level was somehow lower than that of the perfect reconstructed audio. However, I have no explanation as to why this phenomena occurs.


texts:
I have a question about drawing a magnitude characteristics of H(s).

So my transmittance function is for example:

and in the book the magnitude characteristics is shown as:

And my question is - how to count the value of abs(H(jΏ) ?

I understand that s = 1 is a zeros (local minima), and this function has two poles: -1-j and -1 +1 (local maxima), but I can't count the values highlighted in yellow.


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


texts:
Could someone please explain in what case the phase of a real number is equal to -pi (and not pi)?

I know that for positive numbers, the phase is zero. For zero, we define the phase as zero as well. And for negative numbers, the phase would be pi. But I was reading some script and there it says the phase of a real number is either 0, pi, or -pi.


texts:
In GNU Radio Companion, I have created a very simple flow graph where I generate a 1 kHz, 2.0 volts peak-to-peak sine wave, and feed it into an FFT sink and a time domain sink.

The reference dB in the FFT sink is set to 0, and the reference scale is set to 2.0.

As such, I would expect to see 0 dB on the FFT display, but I instead see −10 dB.

Why is this? Is it a problem in my understanding, or a problem in GNURadio?


texts: Image Processing - find the centre of the estimated circle of eye

texts:
I have been recently working on detecting the direction of eye pointed out and done some basic operations on the image of the eye. the below is the image of boundry of the eye.

Now the problem is how do i detect the centre of the rough circle i have obtained. please help me out by providing me a matlab code.


texts: As mentioned in the title, is it possible to use freqz to find the DTFT of a sine wave? I am confused about what the 'a' and 'b' vectors would look like, since there are only impulses in the numerator.


texts: Show That a 2D Linear Transform $ T \left( \cdot \right) $ Is Homogeneous

texts:
By my understanding, a transform T is homogeneous if T[0] = 0.

Then to prove that a linear transformation is homogeneous we say that:

T[ax(n1, n2) + bx(n1, n2)] = aT[x(n1, n2)] + bT[x(n1, n2)]

What I want to know is the difference between the following two proofs and which one is "more correct" or is there perhaps a even better answer?

let x(n1, n2) = 0

then we know that,

ax(n1,n2)- ax(n1,n2) = 0,
T[ax(n1, n2) + ax(n1, n2)] = aT[x(n1, n2)] + aT[x(n1, n2)] = 0

VS

Simply showing that both sides are equal to 0 if x(n1, n2) = 0 through algebra.

On a side note, can someone also comment on why we would even care if a transformation is homogeneous vs heterogeneous? In reference to a transformation applied onto a 2D signal (such as an image transformation).

Maybe the better question to ask would be to prove if a transformation is heterogeneous, is it linear? (if so can someone help me define what it means for a transformation to be heterogeneous? Like how a homogeneous transform means T[0] = 0.


texts: Determining Power of White Noise Component of an Incomplete Time Signal

texts:
I want to determine the power of the noise component of a time dependend signal x(t) (similar to picture below)

With constant sampling rate I use Parceval's Therorem and calculate the standard deviation with

$$
\sigma^2 = \frac{\text{abs}(\text{FFT}(x))^2}{N}
$$

where N is number of data points and fft matlab dft implementation. For the standard deviation according to white noise I remove the discrete frequencies of the spectrum and use above formula with the mean value of the remaining spectrum.

In some signals there are quite a few NaN data points (~ 10%), so said method can not be applied anymore.

Is there a mathematically correct method to determine the power of the noise component of the signal? Interpolation does not seem to be a good solution, because this would lead to an underestimation of the noise.


texts: What's the Q factor of a digital filter's pole?

texts:
For the analog S-plane, the Q factor depends on the angle of the complex pole $p$ from the horizontal axis, $\theta = \arg p$:

damping factor $\zeta = \frac{1}{2Q}$, and $\zeta = \cos(\theta)$, and $\theta = \arg(p)$ so $Q = \frac{1}{2 \cos(\theta)}= \frac{1}{2 \cos(\arg(p))}$  Right?

But in the Z plane, the curves of constant Q are more complicated:

So how do you get the Q of a complex pole?  To transform an analog pole to the Z plane, you use $z = e^{sT}$?  But when I try to apply this I must be doing something wrong.

Oh, maybe I got it.  Maxima was making assumptions about the pole location being real?

$$Q=-\frac{\sqrt{{\log(|p|)}^2+{\arg(p)}^2}}{2\log(|p|) }$$


texts:
I am understanding that mother wavelets have a center frequency. Wavelets are limited in duration and oscillate then decay, so the center frequency is the mid pulse of the wavelet or so.
But my question is, how does the center frequency play a role in noise reduction for example? what does the center frequency affect? does it have to do anything with the resolution of the denoising? and does the center freqeuncy stay the same in discrete wavelets since we do not scale the actual wavelet?


texts: Designing an IIR to generate a specific data sequence

texts:
I know that it's possible to design an IIR with specific poles and zeroes to create specific frequency responses.

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


texts: How do I recover the signal from an ECG image

texts:
In my project, I have to digitize an ECG image taken with a normal camera (jpeg). For example, I have the following camera captured image:

and I want to get something like this:-

and then the digitized data (x,y points), like in this video on digitization of ECG

I have no idea how to do it, so I searched and consulted several research papers. The general approach of the algorithms is like:-

change to gray level image
delete gridlines
add missing points
convert 2D image to 1D image

I am stuck with the second point, i.e., deleting the gridlines. I looked up some more references to do this and found histogram analysis might be helpful.

Can you please guide me on how to do this (I'm using MATLAB 2010)? Any help would be appreciated.


texts:
How wavelet analysis works as a dimension reduction techniques?
The approximations coefficients at higher level of decomposition are the fewer retained coefficients from the originial ?Is it correct?


texts: About how to increase the FIR filter sampling frequency in FPGA? And what is the tradeoff of increasing the sampling frequency?

texts:
I am trying to implement a FIR filter on FPGA and trying to have a solid understanding of the FIR filter tap delay and sampling frequency.

Does the “one tap” delay equal to “1/Fs (sampling frequency)”? If I have N-tap, the total delay will be N/Fs? If the Fs sampling frequency is increased, the “one tap” delay is decreased?

Is there any trade off/drawback of increasing the sampling frequency? Maybe it will take more processing time for the output y(n) to come out?

Thanks


texts: FFT for distance-limited sum of pairwise products

texts:
I have a sequence of fixed length $N$ input arrays ($x_0$, $x_1$, ..., $x_{n-1}$), and my goal is to compute output arrays as such:

$$d_{00}[l] = \sum_{i = 0}^{N - l} x_0[i]x_0[i + l]$$
$$d_{01}[l] = \sum_{i = 0}^{N - l} x_0[i]x_1[i + l]$$
$$d_{02}[l] = \sum_{i = 0}^{N - l} x_0[i]x_2[i + l]$$
$$...$$
$$d_{10}[l] = \sum_{i = 0}^{N - 1} x_1[i]x_0[i + l]$$
$$...$$
With $l$ having a maximum value that scales linearly in $N$. That is to say, $d_{ij}$ will contain, in it's $k^{th}$ entry, the sum of all pairs in $x_i$ and $x_j$ that are distance $k$ away.

My understanding - Clearly, since each $d_{ij}$ output requires $l \propto N$ entries, each containing $N$ multiplies, the time complexity of producing each $d_{ij}$ output array is $O(N^2)$. These pairwise multiplies hidden behind a summation seem optimal for a Fourier Transform to improve the time complexity to $O(N \log N)$. However, I understand the standard Fourier Transform convolution does not allow for limiting to constant distance, instead allowing easy sum of pairwise indices in the form $\sum_{i=0}^{c} x[i] y[c - i]$.

How can I use this standard convolution, or any variation of the Fourier transform algorithm to make this problem computationally tractable?


texts: Help implementing audio dynamic range compression

texts:
I am trying to implement audio dynamic range compression in JavaScript (not using web audio API).

There is a lot of articles for sound technicians, and some high level documentation, but I couldn't find any helpful reference for actually implementing digital dynamic range compression.

From what I understand, there is at least 3 steps in calculating the rectified signal.

computing the input level
computing the gain to apply to the signal
applying the gain

I process the audio in blocks, so for 1) I was thinking of computing the RMS of one block

Any pointers to a good reference? Or anybody willing to explain me a bit the steps required to implement this?


texts: How to update Kalman Filter with n-th order state translation?

texts:
In Kalman Filter, the hidden state translation is defined by $X_t=F_tX_{t-1}+W_t$, where $X_t$ can be a vector or a single value. This is actually derived from Bayes filter, in which 1-th order markov chain is used with $P(X_t|X_{t-1},...X_1)=P(X_t|X_{t-1})$.

Now, in practice, we actually have that $P(X_t|X_{t-1},...X_1)=P(X_t|X_{t-1},...X_{t-n})$, which is n-th order markov chain, then  how we derive its Kalman Filter with both the prediction and correction step?

I did not find any related documents, and also I did not figure out any tricks to solve my problem based on the standard Kalman Filter (i.e., $X_t=F_tX_{t-1}+W_t$).


texts: Sine wave phase shift from Fourier Transform

texts:
This is probably a really basic question but I'm a little stumped and would appreciate some practical input on how to go about doing this rather than reading dockets of equations semi-related to what I am doing.

So say I have data which is largely reminiscent of a sine wave. I do an forward Fourier transform and low pass filter it to remove some noise. Is it possible to then find the phase or shift of the remaining sine wave with the transformed data in the frequency domain before performing the inverse transform?

I did some reading and it looks like the phase shift information is stored in the imaginary part of the transform?

I've written two algorithms for finding the phase shift in the time domain, but they are clunky and unfortunately unreliable. Any help would be appreciated.

Also, I'm not specifically asking for the answer but keywords to research, or places to find information would be highly appreciated.


texts:
I'm starting to learn about OFDM QAM demapping, and a lot of documentation points to using "soft-output" algorithms that calculate Log-Likelihood Ratios (LLRs).

Why is this method advantageous?

My naive approach would be to just place a grid over the IQ map, and find which grid-square my symbol lands in.


texts:
I have a current signal represented by an array in python, which I would like to add shot noise to. Shot noise is defined as
$$
σ_i = \sqrt{2qI \Delta f}
$$
where $q$ is the elementary charge, $I$ the current and $\Delta f$ is the bandwidth of interest. I guess the simplest way to add the noise would be to to add it in the time domain, using the np.random.normal function in python and simply adding it element-wise to the signal array. The frequency is needed to calculate the shot noise, so I would calculate $\Delta f = \frac{1}{2\Delta t}$ where $\Delta t$ is the time between two consecutive samples in the array (correct me if I'm wrong!).

My question is, how can I add this noise in the frequency domain? I would expect to be able to express the shot noise in $\frac{A}{\sqrt{Hz}}$, and therefore make the magnitude independent of the frequency. This post explains the relationship of the magnitude of white noise in time and frequency domain, but if I use that simple equation ($\frac{\sigma^2 M}{2}$ ) it doesn't become independent on the frequency.

Thank you in advance!


texts: How to obtain mean Power Spectrum from two audio channels?

texts:
I have a source of noise in a duct, on either side of the noise source I have a microphone, which register the pressure, so that I get two pressure signals: $p_1$ and $p_2$. I want to get the acoustic Sound Pressure Level (SPL) spectrum which can be estimated using a periodogram approach, more specifically using the welch procedure. This function is readily available in Matlab as

pwelch

or in python using

scipy.signal.welch

The Question

The question is, since I have two channels of audio, how do I combine them to get a representation of the total SPL spectrum? Or is the mean SPL more appropriate?

Some ideas:

I know that for two sound sources giving rise to $p_1$ and $p_2$ then we have two cases:

Correlated sound: $<p_{tot}^2> = \frac{1}{2}(p_1^2+p_2^2)$
Uncorrelated sound: $<p_{tot}^2> = \frac{1}{2}(p_1+p_2)^2$

In my case I'm guessing I should consider the time series $p_1$ and $p_2$ to be correlated since they come from the same source but are measured on separate sides, and this would imply that the total SPL spectrum would be given by averaging the resulting spectra given by welch on the individual channels. Thoughts?


texts: What is the order of the output of HoughCircles in OpenCV, when having several circles on a frame?

texts:
So I am trying to use stereo vision to calculate the depth of several red balls. With the function HoughCircles, I get an output at one frame from one of the cameras like for example:
[x1, y1, r1
 x2, y2  r2
 x3, y3, r3]

where xi and yi (i = 1,2,3) are the coordinates of the center points for respective ball and ri is the radius.
The problem is that HoughCircles for the other camera at the same capture moment, outputs the coordinates in a different order like:
[x2', y2', r2'
 x3', y3'  r3'
 x1', y1', r1']

where xi', yi' and ri' is NOT equal to xi, yi and ri (where i = 1,2,3) because it is a different camera with a different frame angle than the previous one. However, coordinates with the same indexes, like for example x1' and x1, are coordinates for the same ball, viewed from different cameras (Hope that was understandable).
I want the program to recognize exactly which balls are the same, even though the output-matrices are not in the same order for the respective cameras. And it seems like HoughCircles gives different orders in this matrix for different frames.
How do I control or predict the row-order in the matrices?


texts: Segmentation of small artifacts - Image Processing

texts:
I am working on a task of segmenting small artifacts from a set of images like below:

I am looking for a generic approach, but unsure of what to search?

I tried looking at density flow analysis, but don't know if that fit in.

I am unable to come up with a single pattern in all the images. For example I look out for "loops", but few images have loop plus more lines, like in (d) above.

One property I am certain about, the percentage of in these artifacts will be very less compared to the whole image.

Can you help me with:

Any algorithm you know which will solve the problem?
Any research paper, where a similar problem is solved. It is okay if that is not the topic of research paper, as this kind of work is usually present in per-processing stages.


texts:
BACKGROUND

The bit error probability for BPSK under AWGN is easily derived from tail probabilities of Gaussian distributions and results in

$$P_e = Q\biggr(\sqrt{\frac{2E_b}{N_o}}\biggl)$$

The equivalent bit error probability for DBPSK is given as follows but much more complicated to derive:

$$P_e = \frac{1}{2}e^{-E_b/N_o}$$

A complete derivation for the DBPSK case is here:

http://staff.ustc.edu.cn/~jingxi/Lecture%209_10.pdf

With the same formula and plotted in comparison to BPSK on Wikipedia (https://en.wikipedia.org/wiki/Phase-shift_keying#/media/File:DPSK_BER_curves.svg):

MY QUESTION

I incorrectly thought I could simplify this derivation by extending the simpler BPSK $P_e$ result through understanding what occurs when you multiply two signals with independent noise  (Matt L has provided that here: SNR After Multiplying Two Noisy Signals), since such a product results when performing non-coherent demodulation for DBPSK.

I show this in the block diagram below:

This is the non-coherent structure for DBPSK demodulation. The transmitter is also differentially encoded to minimize error propagation (so that errors always occur in pairs rather than propagate until the next transition).

Here we can see that given an input DBPSK signal with $SNR = SNR_1$, the signal after being delayed one bit period $T$ will also have $SNR = SNR_1$, but the noise component will be independent (assuming AWGN, the noise is one symbol period is independent of the noise in the next symbol period). With reference to Matt L's result linked above, the predicted SNR at the output of the multiplier would be:

$$SNR_2 = \frac{SNR_1 SNR_1}{SNR_1+SNR_1+1}$$

For real signals, the frequency at the output of the multiplier is the sum and the difference of the input frequencies, so in this case the difference is baseband signal of interest while the sum is the double of the carrier that we filter out with the low pass filter (LPF). The power of both the signal and noise components would be effected the same way in this process, so the SNR at the output of the LPF would still be $SNR_2$.

Note for SNR>>1, $SNR_2$ approaches $SNR_1/2$ or 3 dB worst.

Given this, combined with the double error property that a single bit error always results in 2 errors assuming we use differential encoding in the transmitter- I can convince myself that the predicted bit error rate for DBPSK would be as follows (reducing the SNR by 2 when SNR >> 1 and doubling the resulting P_e) but from the detailed derivation this is clearly incorrect. I understand the detailed derivation- my question is not with that but what is the flaw with this alternate approach?

$$P_e = 2Q\biggr(\sqrt{\frac{x}{2x+1}}\biggl)$$

where $x = \frac{2E_b}{N_o}$

It is interesting to note that for higher order M-PSK, this 3 dB result does match (note the difference between QPSK and DQPSK in the plot above). Perhaps this is a clue that real versus complex is a factor?


texts: What causes zero crossings in a QPSK constellation to be a a problem for linearity or PAPR ( Peak to Average Power Ratio)?


texts: Extracting rms delay, delay spread, and CIR from power delay profile

texts:
I have a spread spectrum channel sounder running the following transmitter and receiver:

When running it over the air (~2km), I get the plot of the correlator output:

I can see that there are peaks every pn period. How do I go about processing this data to extract rms delay spread, delay spread, and CIR?


texts: Philosophy of perfect inter-sample interpolation

texts:
How would you define perfect (inter-sample) interpolation, and is it possible?

To quote Armifinn's prior answer:

  "I guess the most important result is that for signals with bandwidth limitation, you can have perfect reconstruction via sinc(⋅)sinc(⋅) convolution; the famous sampling theorem..."

If by perfect interpolation it is meant that an analog band-limited signal is recovered perfectly from digital samples, wouldn't a sinc interpolation lead to some problems? Considering that analog equipment is causal while sinc interpolation is anti-causal, shouldn't a minimum phase low pass filter be used instead?

On the other hand, I understand that the above method would not preserve the phase information of the system. Sampling theorem states that band-limited signals can be captured perfectly. So indeed that would be a contradiction.

Finally, a quick thought experiment: consider an analog dirac delta impulse that is low pass filtered using a perfect causal analog filter so that it satisfies the sampling criterion. The signal is then digitally captured by a DSP engineer. Now, the DSP engineer wants to reconstruct the analog signal in his lab, and after running the perfect sinc interpolation he discovers that the signal in fact started BEFORE the analog dirac delta even occurred, to be more precise an infinite time before the analog signal started. This casts some doubts on how the interpolation should in fact be done, and whether the band-limiting criterion is so easy to satisfy.


texts:
I am a PhD. in pure mathematics.

Could you please illustrate the following statement: the eigenvectors of a graph Laplacian behave similarly to a Fourier basis, motivating the development of graph-based Fourier analysis theory.
I am reading the interesting  paper, but could not get how the  Fourier transform is extended to the graph Fourier transform as illustrated on page 23.  Indeed,  why should the matrix V be considered as the extension of (discrete) Fourier transform?


texts:
This tool (Paul's Extreme Sound Stretch) is able to produce extreme timestrech (50x, or even 100000x). Any song can become an interesting sound texture or ambient music (listen to this example).

The technique used is mentionned as:

  "spectral smoothing" the sounds.

How does an extreme time stretch algorithm work?

(it is opensource, I have studied it a bit, but I haven't been able to see a big picture of the algorithm yet)

I thought about :

sound input ====>  STFT ===> duplicate the frames 50x ===> ISTFT ===> output

but I'm really not sure about this.


texts: Clarification on how bandwidth effects input signal

texts: I know that the system bandwidth is the maximum range of frequencies that a system can handle. Is the reason that the input signal gets distorted as the bandwidth decreases a result of fewer segments of the signal being passed through the system? Does this mean that if I have a bandwidth greater than $B_{h1}$ I would have an output with no distortion?


texts: Column operations in LDPC generator matrix

texts:
I'm trying to produce a generator matrix from a starting low density parity check matrix. There are lots of references on the topic (including this Signal Processing Stack Overflow answer here). They make sense, with the exception of one problem I'm having.

My resource is the book Error Correction Coding: Mathematical Methods and Algorithms by Todd K. Moon. In it, there is a section that describes this that I will replicate below, and then highlight the part I'm stuck on.

  For a code specified by a parity check matrix $A$, it is expedient for encoding purposes
  to determine the corresponding generator matrix $G$. A systematic generator matrix may
  be found as follows. Using Gaussian elimination with column pivoting as necessary (with
  binary arithmetic) determine an $M$ x $M$ matrix $A^{-1}$ so that

  $$H = A_p^{-1} A = [I \quad A_2].$$

  (If such a matrix $A$, does not exist, then $A$ is rank deficient, $r = \mathrm{rank}(A) < M$. In this
  case, form $H$ by truncating the linearly dependent rows from $A_p^{-1} A$. The corresponding
  code has $R = K/N > (N - M)/N$, so it is a higher rate code than the dimensions of $A$
  would suggest.) Having found $H$, form

  $$G=\begin{bmatrix}A_2 \\ I\end{bmatrix}.$$

  Then $H G = 0$, so $A_p H G = A G = 0$, so $G$ is a generator matrix for $A$. While $A$ may be
  sparse (as discussed below), neither the systematic generator $G$ nor $H$ is necessarily sparse.

Here's the problem I'm running into:
I start with $A$ and do Gauss-Jordan elimination in $\mathbb{G_2}$ to get it into reduced row echelon form.

Sometimes the matrix I've been given (I'm not designing it, have to use what I'm given) can't be put into that form with solely row operations.

It requires actually swapping the order of some columns (is that what he means by "column-pivoting" in the reference above?). Everything in me screams that's wrong, but it turns out if I just do the same thing on the decoding side it works out.

But then I'm left with the conundrum that when I go to decode something encoded with this parity check matrix I have to know how the matrix columns were reordered on the transmit side.

Basically I have to re-do the generator matrix, which is by far the slowest part of my process (I'm going to make some improvements there though, using the QC-LDPC structure of the parity matrix).

It seems bizarre to me that a designed system would put in a parity check matrix that requires column reordering to work.

Is there something I'm missing? For what it's worth I also have the example matlab code from the aforementioned book and this is what he does (column reordering and then re-applying that re-orderying on the decode side).


texts: Is a good practice to use one PID in rise time and another PID for steady time in order to eliminate oscillations during steady time?

texts:
I'm trying to control temperature with PID(Kp,Ki, Kd). During the rise time the system use one PID and when target is reached the system switch to another PID with other parameter in order to eliminate the oscillation. I know that I can remove oscillations increasing Kd but in this case it does not work.

Results of real time system

Control temperature with constant PID:

Control temperature using PID1 during rise time (blue color) and when temperature reach the target, the controller is changed to PID2 (Steady state error: black color):

void setTarget(...)
{
    mPid = pid1;
}

void steadyStateError()
{
    float e =  fabs( temperature - target);
    if(e<=0)
       mPID = pid2;
}


texts:
What does 'canonical' mean?
e.g. 'write IIR or FIR in canonical form'

I get the gist of it but maybe someone here has a precise definition. In a non-DSP context it would mean something like 'strictly by the book', and perhaps in a DSP context this would mean 'without considering implementation issues'. Is this true?


texts:
I am studying control system and I have encountered the topic of the inversion of the dynamics. So I have seen that the ideal situation would be

$C(s)=P(s)^{-1}$

but there are some problem with doing that and one of these is the fact that if P(s) is proper, its inversion would be improper. Now, I have also seen that to try to overcome this problem we could place poles at high frequencies in order to make the inverse proper.

One thing I really cannot understand is why do we do this? So why do we need to place the pole at high frequencies and not a low frequencies?

I don't understand why necessarely at high frequencies.

Can somebody help me? Thank's in advance.

[EDIT] My idea is that if we insert a pole we are creating a low pass filter and so it is better to insert it at high frequencies beacuse at low frequency there is our signal of interest, while at high frequrncy there is noise. So if we insert the pole a low frequency we touch also the signal we care about, while if we insert at high frequency we leave it untouched and the filter cuts only noise.

Is it correct?

Thanks in advance


texts: How Does Mean Centering Affect the Result of Using SVD to Compress Images?

texts:
I have been learning about using the Singular Value Decomposition to find low rank approximations to matrices.
I had an image which I converted to a matrix. I regarded each row of the matrix as a 'data point'.

I did two things:

Found the mean of the data points then subtracted the mean from each data point and calculated the SVD of the result
Directly calculated the SVD of the original data matrix

Then I calculated the low rank approximations of the results from 1. and 2. by replacing some of the singular values with zeros.
Finally I converted the results back to images.

Both images looked quite similar. They both looked like how SVD-compressed images should look (according to textbooks).

From a principal component analysis point of view, not mean-centering the data before finding the low-rank approximation shouldn't
work well. It is equivalent to calculating the principal components from the eigenvectors of $\frac{1}{n-1}X^TX$ (where $X$ is the data matrix) instead of the covariance matrix
of $X$. So why did it appear like mean-centering didn't matter? Or does it depend on the type of images?


texts: What are the advantages of barker sequence over marker sequnce?

texts: I read that barker sequence is used for synchronization among other uses.  Barker sequences are rare and they have special properties - this part I understood.  I am missing some thing when I think about how they score over simple marker sequences for synchronization.  You want to know where the information bits start - start with a marker sequence preferably which does not occur in the information.  And you can synchronize.  How do barker sequences score over marker sequences?


texts: Getting the frequency of a 0.33Hz to 3Hz signal

texts:
I read a noisy signal from an accelerometer which is attached to the pedal of a stationary bike and need to find the frequency of the signal/ speed of the person cycling.
I taught of the following solutions and would like to know if I overlooked any or made some incorrect assumptions.

Currently the working solution is a low pass (10-30 point moving average) and then a peak detection with a minimum amplitude requirement and minimum time requirement. However this solution doesn't work very well when nearing 0.33Hz (20 RPM)

I tried using an fft but the result was never accurate enough for the small differences in frequency I am looking for.

I taught about using autocorrelation but how I understand it, I would have to wait atleast a period and then some before I could detect the frequency.

Any other ideas are welcome, or if anybody knows what is used in the commercial version it would be great to know.


texts: Butterworth band-pass filter transfer function

texts:
I have some complex data in the frequency domain and want to multiply by an appropriate function to get complex data in the frequency domain that is appopriately filtered.

I think such a function should be, in terms of frequency:

$$ 1 - \frac{1}{\sqrt{1 + (\frac{fW}{f^2 - f_0 ^2})^{2n}}} $$

where W is the width of the frequency band to pass, and $f_0$ is the center frequency.

Is the above correct, or what would be the correct function to accomplish this?