texts:
Forewords

This question is about methodology references and numerical application. I am posting on Signal Processing because I think this question belong to this place. I am new to the stack, feel free to edit if you feel the Original Post should be improved.

Experiment Setup

Basic experiment setup: inject a dye into a reactor to analyze its dynamic. There are two detector cells (sensitive and selective to dye substance) at in entrance and another at the exit of the reactor. No reaction occurs, the dye is a passive tracer in this experiment.

Two signals are recorded:

$x(t)$, a finite impulse (eg.: we realized a gaussian-like impulse, at the entrance);
and $y(t)$, the response of a system to this impulse (at the exit);

Goal is to compute $g(t)$: the system Retention Time Distribution (RTD).

Signal Processing

For the MCVE sake, I created a synthetic dataset to perform computation on well controlled scenario without noise.

Where $y(t)$ is a convolution such as:

$$y(t) = g(t)*x(t) \xrightarrow{\mathcal{F}} Y(\nu) = G(\nu)\cdot X(\nu)$$

Now, it might resume to get the inverse transform of the transfer function:

$$ G(\nu) = \frac{Y(\nu)}{X(\nu)} $$

But numerically, even for this smooth example, it does not work well. Mainly because $X(\nu)$ has zeros and small numbers that make the operation unstable.

Alternative

There is a processing (adapted from this example) that seems to solve my problem:

$$ G'(\nu) = \frac{X^*(\nu)\cdot Y(\nu)}{X^*(\nu)\cdot X(\nu) + \lambda^2} $$

Where $\lambda$ is related to Signal-to-Noise Ratio:

$$ \lambda^2 = \frac{1}{\mathrm{SNR}} $$

I feel comfortable with this solution, because the operation tends to the former when $\lambda$ goes to zero:

$$\lim\limits_{\lambda \rightarrow 0} G'(\nu) = G(\nu)$$

Then I can estimate RTD by taking the inverse transform of $G'(\nu)$:

$$\hat{g}(t) = \mathcal{F}^{-1}[G'(\nu)]$$

Numerical Example

For the MCVE dataset, it renders as follow:

Derived RTD (after inverse transformation) quite fits the original function:

And error on this method is acceptable for the MCVE:

Even if there are artifacts at the very beginning and end of the recovered signal.

Observations

Setting $\lambda$ to higher values change the error behavior, but the RTD looks fine:

Setting  $\lambda$ to smaller value may affect the shape of the RTD:

Questions

Reading on the internet I had difficulties to find references where the transfer function is the unknown. Most of the articles speaks about deconvolution and estimation of $\hat{x}(t)$ knowing $y(t)$ and $g(t)$.

I also found the Wikipedia page of Wiener Deconvolution which shows a formulae that is very similar to $G'(\nu)$ in my OP. But I could not relate them together. I miss something to get the full understanding on why this methodology is correct, if it is, when performed as I did.

My questions are:

How can I formally justify the expression of $G'(\nu)$? Can it be shown there is a relation with Wiener Filter?
Is there a good reference to apprehend this class of problem?
How is called the oscillating effect at the ends of the signal?

Your help is greatly appreciated.

Code

Below the Python code of the MCVE:

import numpy as np
from scipy import stats, signal, fftpack

# Setup:
fs = 10   # Sampling Frequency [Hz]
dt = 600  # Experiment runtime [s]
l = 1e-4  # 1/SNR in term of amplitude [-]
t = np.linspace(0, dt, dt*fs+1) # Time [s]

# Signals [AU]:
xt = stats.norm(loc=20, scale=5).pdf(t)
gt = stats.lognorm(0.45, loc=-20, scale=100).pdf(t)
yt = signal.convolve(xt, gt, mode='full')
yt = yt[:t.size]/np.sum(yt[:t.size]/fs) # Normalization

# Fourier Transform:
Xs = fftpack.fft(xt)
Ys = fftpack.fft(yt)
fq = fs*fftpack.fftfreq(t.size)
Gs = (np.conj(Xs)*Ys)/(Xs*np.conj(Xs) + l**2) # Transfer function

# Inverse Transform:
gti = fftpack.ifft(Gs)
np.allclose(np.imag(gti), 0) # True
gti = np.real(gti)
gti = gti/(gti.sum()/fs) # Normalization

# Error:
ea = gti-gt
er = ea/gt


texts: On the meaning of s-plane and it's link to a transfer function

texts:
Considering Fourier analysis and let's say I'm walking on the blue frequency axis in the below 3D plot from zero towards infinity:

So each time I encounter a non zero blue bar, I check the frequency at that point and measure the bar's height and say: Oh I found another sinusoidal component of the signal at this particular frequency and at this particular amplitude. This is how I grasp the idea of Fourier analysis.

But when I read about the s-plane I'm extremely confused to make it visual and link it to the actual signal as in the above fashion. I just cannot give a similar meaning what it is that really about.

Now imagine a system's output and input in time domain is given. Vout(t)/Vin(t) is transformed into s-domain and we call this ratio as H(s) which is the transfer function.

Here is a transfer function which is a surface on the s-plane:

Imagine I am walking on the s-plane. And imagine at a point on the s-plane I stop and the point s is a complex number and it is s = ωj + σ. Now if I look upwards towards the surface above me, obviously there will be a point on the surface. Lets call this point x. So what should it mean to me in terms of the signal or the transfer function?

Ok so at this particular point I found "what"? What is "x"? Is it a component of something? And my second question is if H(s) were a signal's transform not a transfer function what would x be in that case?


texts: How to intuitively understand the state space formulation of discrete time system?

texts:
The SS formulation of DT system is given by

$$x[(k+1)T] = Ax(kT) + Bu(kT)$$

$$y(kT) = Cx(kT) + Du(kT)$$

Note: T is the sampling period and often omited

Can someone explain to me why the state equation $x[(k+1)T] = Ax(kT) + Bu(kT)$ is only between the currently sampling period $kT$ and the next sampling period $(k+1)T$?

How can this model be used to describe the state occurring at all times?


texts: MPEG-7 descriptors' implementation in OpenCV

texts:
I'm looking for implementation of MPEG-7 descriptors which are compatible with OpenCV's recent versions and containing "region-based shape descriptor ART(Angular Radial Transformation)".

I encountered an stackoverflow question which had very informative answers:

Link in one answer was containing only the following descriptors:

Color Structure Descriptor (CSD)
Scalable Color Descriptor (SCD)
GoF/GoP Color Descriptor (GoF/GoP)
Color Layout Descriptor (CLD)
Dominant Color Descriptor (DCD)
Homogeneous Texture Descriptor (HTD)
Edge Histogram Descriptor (EHD)
Face Recognition Descriptor (FRD)

Link in another answer was looking too outdated and wasn't containing documentation.

If there is another implementation I would like to go through it.

Thank you.


texts: Image objects dispersion and diversity measures

texts:
I have images like below:

I want to find some measures describing black spots dispersion and diversity.

First measures I use are mass center off all black spots, mean area of spot, stdev of area need more ideas.

I want measures to describe irregularity, differences etc.


texts: Downsampling and Convolution: When Can one Swap Between Them?

texts:
(Code snippet with example below).

In MATLAB, assume I have a long signal vector y (length N) that I wish to convolve with a square wave h (consisting of H 1's). This would give a convolution result of length N+H-1.

Since I have to perform this operation quite often, I thought that I could speed up my code by downsampling the signal y (and also the square wave h) and then perform the convolution. The downsampled signals are referred to as y_ds and h_ds.

However, the result is shifted when compared to the downsampled version of the convolution result (if corrected for scaling differences). So in general I found that:

  conv(y_ds,h_ds) $\neq$ downsample(conv(y,h))

However, the shift is between the two results is not constant, i.e. shifting such that the maxima align, causes the edges to be different.

My questions

Where does this shift originate from?
How can this code be modified such that the downsampled convolution result equals the convolution result of the downsampled signals?

Kind regards

N = 100;
dsfactor = 3;
H = dsfactor*3;
x = linspace(0,10,N);
y = sind(18*x);
h = ones(H,1);

convyh = conv(y,h);

convyh_ds = downsample(convyh,dsfactor)./H;
convyh_ds3 = conv(downsample(y,dsfactor),downsample(h,dsfactor))./(ceil(H/dsfactor));

figure;
plot(convyh_ds,'DisplayName','downsample(conv)','LineWidth',1);
hold on;
plot(convyh_ds3,'DisplayName','conv(downsample)','LineWidth',1);
grid minor
legend('show');
line([0 length(convyh_ds2)],[0 0],'LineStyle','--');


texts: Difference between power spectral density, spectral power and power ratios

texts:
What 'exactly' is power spectral density for discrete signal?

I was always under the assumption that taking the Fourier transform of the signal, and then the ratio of desired frequency range magnitude over the entire frequency range gives the power ratio for that frequency range which is the same as power spectral density. Is that wrong?

Reading a student paper got me confused as it says to compute PSD and then 'absolute and relative spectral powers in desired bands' as well. Are they different? If yes, how does one compute it?


texts: Complex signals is the same as natural signal

texts:
I am working using a sinuasidal signal along with the complex signal derived form it and when I plot the signals spectrum  using FFT I got the same spectrum and it turns out that the signals are the same which I don't get?
import scipy
from scipy import fftpack
from scipy.signal import butter, lfilter,filtfilt
from scipy.signal import freqz
import numpy as np
#from scipy.signal import butter,filtfilt# Filter requirements.
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import pandas as pd

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    print()
    print("low is ",low)
    print("high is ",high)
    print()
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b, a, data)
    return y
freqs = [7,20,30]
sins = ['sin', 'cos', 'sin']
amps = [1,0.5,0.2]
phases = [0, np.pi/4, np.pi/3]
def m_signal(t, freqs, sins, amps, phases):
  ops = np.zeros(t.shape)
  for f,s,a,p in zip(freqs, sins, amps, phases):
    if s == 'sin':
      ops += a*np.sin(2*np.pi*f*t + p)
    if s == 'cos':
      ops += a*np.cos(2*np.pi*f*t + p)
  return ops

fs = 240
T = 2
nsamples = fs*T
time = np.linspace(0, T, num=nsamples)
sig0 = m_signal(time, freqs, sins, amps, phases)

def time_signal(n):
  return np.sin(2*np.pi*n*7) + np.cos(2*np.pi*n*30 - np.pi/4)

# Complex Signal
def signal(n):
  s1 = np.exp(1j*2*np.pi*n*7)/(2*1j) - np.exp(-1j*2*np.pi*n*7)/(2*1j)
  s2 = 0.5*np.exp(1j*(2*np.pi*n*30 - np.pi/4)) + 0.5*np.exp(-1j*(2*np.pi*n*30 - np.pi/4))
  comp = s1 + s2
  return comp

# If fs is 240 Hz and T (time) is 2 sec, then the number of samples is fs*T
fs = 240
T = 2
nsamples = fs*T
time = np.linspace(0, T, num=nsamples)

time.shape

def plot_spectrum(nsamples, x, fs):
  Xf = fftpack.fft(x, nsamples)
  ff = np.linspace(0.0, fs/2, int(nsamples/2))
  # print(ff)
  fig, ax = plt.subplots(figsize=(20,3))
  ax.plot(ff, 2.0/nsamples * np.abs(Xf[:nsamples//2]))
  plt.grid(True)
  plt.show()

sig = time_signal(time)
csig = signal(time)

sig.shape

# Signal on time domain
plt.figure(2, figsize=(20,3))
plt.clf()
plt.plot(time, sig, label='Signal')
plt.grid(True)
plt.show()

sig

csig

# Complexified signal (Magnitude)
plt.figure(2, figsize=(20,3))
plt.clf()
plt.plot(time, np.abs(csig), label='Signal')
plt.grid(True)
plt.show()

# Complexified signal (Phase)
plt.figure(2, figsize=(20,3))
plt.clf()
plt.plot(time, np.angle(csig), label='Signal')
plt.grid(True)
plt.show()

# Spectrum of the signal sig
plot_spectrum(nsamples, sig, fs)

# Spectrum of the complexified signal csig
plot_spectrum(nsamples, csig, fs)

print(csig)


texts: Best (Perceptually / Objectively) Super Resolution Methods Out There?

texts:
I'm curious about the advances in the area of image super-resolution (SR) that have given the best results to date, both perceptually (visually pleasing) and objectively (e.g. PSNR, SSIM criterias). I'm mostly talking about Single Image SR (SISR), but I'm open to input on Multiple Image SR as well.

From what I've researched, Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network looks like the one, at least as far as SISR goes. However, I am unsure as to what alternatives are out there, that come close or maybe even surpass it (both in the field of neural networks and other, more traditional fields).

Remark
If this reaches anyone's interest but my own, perhaps this can be a question that gets new answers as time and advances move forward.


texts: Optimal signal sampling reduction method : Decimation, interpolation or both?

texts:
I have a function with an equation:
\begin{equation}
C(t)=1.6925\left(\exp^{-0.136t}- \exp^{-1.192t} \right) u(t)
\end{equation}

where u(t) is a unit-step function. Then, denoting  the Fourier Transform of C(t) as C(F), we have
\begin{equation}
C(f)=1.6925\left(\frac{1}{0.136+2\pi jf}- \frac{1}{1.192+ 2\pi jf} \right)
\end{equation}

A time and signal representation of the signal is given by the below figures:

We can see that the C(f) is almost zero after 0.8 Hz therefore sampling at 1.6 Hz based on the Nyquist criterion, we can obtain theoretically a good sample of the original signal.

However 1.6 Hz, is about a sampling period of 0.625 seconds, which in my context of biological studies that is way too frequent. I am looking at other avenues to reduce the sampling rate or increase the sampling period and a search in the litterature has yielded decimation as a possible avenue. A quick hack in matlab using the decimate function yields the below figures for a decimated signal after being sampled at 1.6Hz

My end goal is to compute the surface below my original signal within a time interval or its sampled or decimated versions, do you think that decimation can allow me to do that? How can i measure aliasing between the different versions of my signal? Is there a better way to sample that signal and reconstruct it from a limited number of samples?

Here is my current Matlab Code

clear all; clc;close all;

%% Time-domain signal
Ts=0.625; % Sampling time (obtained from the frequency domain plot)
t=0:Ts:60; % in sec
C=1.6925* (exp(-0.136*t)-exp(-1.192*t));
clog= log(C);
figure,plot(t,C);
figure,plot(t,clog);
title('Time domain signal');
xlabel('Time in sec')

%Frequency domain signal
f=0:.1:1.6;
F1=1.6925*((1./(0.136+(j*2*pi*f)))-(1./(1.192+(j*2*pi*f))));
figure,plot(f,abs(F1));
title('Frequency domain signal');
xlabel('Frequency in Hz')
% Via FFT
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz')

%Downsampling a signal and decimation with interpolation
y = decimate(C,10);

subplot 211
stem(0:48,C(1:49),'filled','markersize',3)
grid on
xlabel 'Sample number',ylabel 'Original'
subplot 212
stem(0:4,y(1:5),'filled','markersize',3)
grid on
xlabel 'Sample number',ylabel 'Decimated'


texts: Does a simple photograph contain more information than a complex painting?

texts:
I hope this question is appropriate for this site.

I came across this passage in The Three Body Problem, a novel by Liu Cixin:

  The professor had put up two pictures: One was the famous Song Dynasty painting Along the River During the Qingming Festival, full of fine, rich details; the other was a photograph of the sky on a sunny day, the deep blue expanse broken only by a wisp of a cloud ...
  The photograph's information content - its entropy- exceeded the painting's by one or two orders of magnitude

Representative pictures:

Is this true? How does one explain this counterintuitive phenomenon?


texts: Underdetermined deconvolution of windowed output

texts:
Consider a discrete 'blurred' output $h[t]$ given by the convolution of filter $f[t]$ and signal $g[t]$. This question considers recovering $g[t]$ from a window (subset) of $h[t]$. This causes the problem to be under-determined.

$$h[t] = (f * g)[t]$$

Original signal $g[t]$:

$$g[t] = [1,5,3,5,7,7]$$

Moving average filter $f$ of window size 2:

$$f[t] = [0.5,0.5,0,0,0,0,0]$$

The full 'blurred' output $h$ is given by:

$$h[t] = [0.5,3,4,4,6,7,3.5]$$

Given $h$ the problem of recovering $g$ is well-posed and can be recovered simply using Fourier transforms:

$$g[t] = \mathcal{F}^{-1}\left\{\frac{\mathcal{F}\left\{h[t]\right\}}{\mathcal{F}\left\{f[t]\right\}}\right\}$$

However in practice the recovery needs to be made from a windowed output $h'[t]$:

$$h[t]' = [3,4,4,6,7]$$

The problem of recovering $g[t]$ from $h'[t]$ is ill-conditioned since there are an infinite number of solutions. One approach I have seen suggested is to take $g[1] = h'[1]$. In this case the recovered profile $g'[t]$ is given by:

$$g'[t] = [3,3,5,3,9,5]$$

My understanding is that most approaches are based on regularisation however most sources seem to focus on ill-conditioning resulting from additive noise rather than ill-conditioning resulting from a windowed output. What are the approaches typically used for this problems? Can you apply it to this toy example?


texts: Simplification after sampling an analog signal

texts:
$80\textrm{ Hz}$ Sampling rate produces the $x[n] = 3\cos(5\pi n/4)$ signal, the original analog signal is $x(t) = 3\ cos(100\ π t)$.
This is a solved problem used as an example. In the solution it's also mentioned that the resulting signal can be simplified for $|f| ≤ 1/2$ as, $x[n] = 3 cos(5 π n / 4) = 3 cos[(8 – 3) π n / 4] = 3 cos(2 π n – 3 π n / 4) = 3 cos(3 π n / 4)$.

I do not understand how $3cos(2 π n – 3 π n / 4)$ becomes $cos(3 π n / 4)$.
I know that $cos(2πn-3/4πn) = cos(2πn)cos(3/4 πn)+sin(2πn)sin(3/4πn)$.
And $cos(2π) = 1$ , $sin(2π) = 0$ if it wasn't for the "n" variable then yes i think it would equal $cos(3 π/ 4)$.
Why and how does this simplification happen?


texts: Why are two cameras mounted in paralle in stereo application?

texts: For stereo cameras on the market, two cameras are always mounted side by side and with a displacement that is perpendicular to the cameras’ optical axes. I take this setup for granted. One idea came to my mind whether this is necessary? If two cameras are not parallel and have different focal length, camera calibration can correct the difference. Why are two cameras mounted in parallel? My guess is that two cameras can have a large overlapping region. Am I correct?


texts: Solving linear system and find impulse response

texts:
I have this linear system that is defined by the differential equation
    $$
    y''(t) + 4y'(t) + 5y(t) = 2x(t) + 3x'(t)
    $$
with x(t) the input and y(t) the output. I'm asked to find two things: The output when x(t) = cos(t), and the impulse response

This is what I did: First, I found the frequency response,H(f)
    $$
    H(f) = \frac{2 + 6\pi fj}{8\pi fj -4\pi^2 f^2+5}
    $$

Second, I found the Fourier Transform of cos(t) which equals to:
    $$
    X(f) = \frac{1}{2}\delta (f-\frac{1}{2\pi }) + \frac{1}{2}\delta (f+\frac{1}{2\pi })
    $$

Now, the output basically would be
    $$
    Y(f) = H(f)\ X(f)
    $$
Is it possible to simplify the answer of Y(f)?

For the second part of the question, I found the impulse response using Laplace Transform
    $$
    h(t) = [\frac{-17}{6} e^{5t} - \frac{1}{6} e^{-t}]u(t)
    $$

What do you think about my answer? Is it possible to use Z-transform to find h(t)?


texts:
I have a question that might be stupide !

Supposing that I try to detect object moving in a video or human action.

Many works are based on optical flow computation.

My question is why using OF is the best for motion detection. A tracking of pixels or region of interest will not be as informativa as OF .?


texts: How to calculate $E_b/N_0$ for a 8PSK? My BER is given. But I don't understand how to calculate $E_b/N_0$ as it is inside a Q function.


texts: Relationship between filters, Fourier magnitude, and frequency

texts:
I'm just beginning to study Fourier transforms and I'm having trouble getting this concept down.
I know what a low/high pass filter does to a signal.

What is the relationship between Fourier magnitudes (dB) and frequency (Hz)?
How would Fourier magnitudes get affected when a low/high pass filter is passed on a signal?
How about an all pass filter?


texts: Question on spatial sampling for non-linear and/or non-uniform arrays

texts:
I am a little confused as to how to compute the nyquist criterion for a sensor array that is not linear, and/or non-uniform. For the linear, uniformly spaced case, the spatial nyquist is delimited by finding the half-wavelength that fits between the space of two adjacent sensors. Thus, if the distance between two sensors is given by $d$, then the highest frequency that is allowed is given by:

$$d < \frac{\lambda}{2} = \frac{c}{2f}$$

Where $c$ is the speed of light, and $f$ is the frequency of the signal. And so therefore, the highest frequency allowed such that no spatial aliasing occurs is given by:

$$
f < \frac{c}{2d}
$$

So my question is twofold:

I) What is the spatial nyquist for arrays that are uniformly spaced, but non-linear, (for example, 10 sensors evenly spaced on a circle)?

II) What is the spatial nyquist for cases where the array is non-uniform and non-linear. (for example, some odd shaped number of sensors placed arbitrarily).

I think it has to do with the largest distance between them, but I am not clear on this, because, well, what distances do you care about?

Thanks


texts:
the voss' algorithm generated a 1/f series. What does the series mean?
For example, we use three dice and the result sequence length is 8. Let's set the sequence is 8,15,6,9,14,10,9,7. Then what does this sequence mean? does it mean the sequence of frequency or something else?


texts: FFT of long signal by segments/chunks with discontinuities

texts:
I am processing an EEG brain signal, which has up to 64 data channels sampled at 500 Hz. One of the analyses consists of extracting the ratio of alpha/delta power, where alpha represents the waveforms which frequencies are between 8 and 13 Hz, and delta the waveforms which frequencies are between 1 and 4 Hz.
The signal is first filtered with a 1-15 Hz BP filter (scipy Butterworth 4th order, output 'sos'). Additionally, a common average projector and EOG SSP artifact correction projector are applied. Those are not relevant to this question.
My program compute the alpha/delta ratio at 2 different moments:

A short window (500 ms or 1 second, i.e. 250 or 500 samples)
A long window (4 to 10 seconds)

However, EEG signal is very susceptible to noise, i.e. if the participant moves his eyes, jaws, head, .. the brain signal will be completely masked by artifacts. Thus, a rejection criteria is also applied. It's a simple peak-to-peak rejection, i.e. if the max - min is larger than the criteria on at least one of the channels, this window (epoch) is rejected.
For the short window, I can directly implement this as follows:
if any(np.ptp(data, axis=0) > reject['eeg']):
    continue # skip

With data a channel x samples array (e.g. 64x250) and reject['eeg'] the threshold value.
However, for the large window, I can not disregard the complete window. My idea was to cut down the e.g. 5 seconds signal into 1 second chunks, and apply the rejection criteria on the chunk.
Something like:
data = raw.get_data(picks='eeg') # np.ndarray 64x2500 for 5 seconds
data = data.reshape(64, 5, 500) # I hope this is the correct reshape, not tested.
chunk2remove = np.any(np.ptp(data, axis=2) > reject['eeg'], axis=0)
data = data[:, np.where(chunk2remove == False)[0], :]

However, now I do not know if I can safely apply FFT on the signal. I do not know what the impact of this chunk rejection is, nor what the border effect might be. What is the correct way to apply an FFT to this chunked signal with discontinuities between the chunks?
Currently, the method I use to extract the alpha/delta power ratio is:
fs = 500.
window = 5 # window length in seconds
fft_freq = np.fft.rfftfreq(int(fs * window), 1.0/fs)
alpha_band = np.where(np.logical_and(fft_freq>=8, fft_freq<=13))[0]
delta_band = np.where(np.logical_and(fft_freq>=1, fft_freq<=4))[0]

fftval = np.abs(np.fft.rfft(data, axis=1) / fs)
alpha = np.average(np.multiply(np.abs(fftval[:, alpha_band]).T, weights))
delta = np.average(np.multiply(np.abs(fftval[:, delta_band]).T, weights))

The array weights is simply a [0, 1] weight applied on each channel. As you can see, this is not exactly the power, but as I'm doing a ratio, I did not feel that the squaring was required. Feel free to correct me if I'm wrong.
If you have any comment, additional ideas for improvement, I'm also interested.


texts: Parameters of Butterworth filter that approximates the amplitude response of a Gaussian filter?

texts: For a given narrowband Gaussian filter with a specific central frequency and filter width, I need corners of a bandpass Butterworth filter whose amplitude response is close enough to the Gaussian filter. How should I choose the order and corner frequencies of the Butterworth filter ? The reason is I need to calculate statistics of the initial part of a time-series filtered with a Gaussian filter. However, the Gaussian filter is zero-phase and always shifts energy back and contaminates the initial part of the time-series. With the Butterworth filter, I can have a one-pass filter so that the filtered time-series is causal and there is no contamination of the initial part of the time-series.


texts: How to prove a train of sinc pulses in digital communicaton system are orthogonal to each other?

texts:
Consider a train of sinc pulses:
$$\phi_n(t)= \frac{\sin(\omega_M(t-nT_s))}{\omega_M(t-nT_s)}\quad; n=0,\pm1,\pm2,\dots$$
$\quad$where,$\quad T_s=\frac{\pi}{\omega_M}$
Now ,in order to show sinc pulses are orthogonal we need to prove:
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=T_s \delta_{nk} \quad \dots(1)$$
where, $\delta_{nk}$ is kronecker's delta.
So, i began doing it as follows:
$$T_s=\frac{\pi}{\omega_M}=\frac{\pi}{2\pi f_M}=\frac{1}{2f_M}=\frac{1}{f_N} \quad \dots(2)$$
where, $f_N$ is the Nyquist frequency
$$\phi_0(t)=\frac{\sin(\omega_Mt)}{\omega_Mt}=\frac{\sin(2\pi f_Mt)}{2\pi f_Mt}=sinc(2f_Mt)=sinc(f_Nt) \quad \dots(3)$$
Now, $\mathscr{F}\{ sinc(f_Nt) \}=\frac{1}{f_N} rect(\frac{f}{f_N})$ ,where $rect$ is a rectangular function centred at origin and having width= $f_N$
$$\implies \mathscr{F}\{ sinc(f_N(t-nT_s)) \}=\frac{1}{f_N} \exp(-i2\pi f n T_s) rect(\frac{f}{f_N}) \quad \dots(4)$$
Now we can write:
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=\int_{-\infty}^{\infty} \{ \Phi_n(f) \circledast \Phi_k(f) \} df$$
$$=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\Phi_n(\tau) \Phi_k(f-\tau) d\tau df$$
$$=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s)
\frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(-i2\pi (f-\tau) k T_s) \frac{1}{f_N} rect(\frac{f-\tau}{f_N}) d\tau df$$
$$=\int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(i2\pi \tau k T_s) \{ \int_{-\infty}^{\infty} \exp(-i2\pi f k T_s) \frac{1}{f_N} rect(\frac{f-\tau}{f_N}) df \} d\tau \quad \dots(5)$$
The inner integral of $(5)$ can be simplified as:
$$\int_{\tau -\frac{f_N}{2}}^{\tau +\frac{f_N}{2}} \frac{1}{f_N} \exp(-i2\pi f k T_s) df$$
$$=\frac{\exp(-i2\pi \tau k T_s) \sin(\pi k)}{\pi k} \quad \dots(6)$$
So, $(5)$ can be rewritten as:
$$\int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N}) \exp(i2\pi \tau k T_s) \frac{\exp(-i2\pi \tau k T_s) \sin(\pi k)}{\pi k} d\tau$$
$$=\frac{\sin(\pi k)}{\pi k} \int_{-\infty}^{\infty} \exp(-i2\pi \tau n T_s) \frac{1}{f_N} rect(\frac{\tau}{f_N})d\tau $$
$$=\frac{\sin(\pi k)}{\pi k} \frac{\sin(\pi n)}{\pi n} \quad \dots(7)$$
Now, $(7)$ is even equal to $0$ when $k=2$ and $n=2$
So,where i missed ? any help or suggestions please...


texts: How do we know that CTFT of the autocorrelation function is the PSD?

texts:
I know that the Fourier Transform of the autocorrelation function is
the Power Spectral Density. But how can we arrive at such a result intuitively?
Is it just a theorem?


texts: Why motion cannot be estimated in DCT domain using block match?

texts:
In all video coding standards, motion estimation is always done in pixel domain using block matching. I wonder why it cannot be estimated in DCT(Discrete cosine transform) domain using block matching directly? If two blocks are similar in value in DCT domain, doesn't it mean that they are similar to each other in pixel domain as well?

All the comments and suggestions are welcome!!!
Thank you!


texts: Higher side lobe levels in linear array beam pattern

texts:
I am working on beam forming of high frequency (>400 KHz) linear array for imaging sonar. I have gathered sensor data of 80 channel linear array sub merged halfway in 8 m deep water tank with a transmitting probe 6 m away from it at same depth. I recorded multiple pulses and perform offline analysis on it. The problem i am facing is very high side lobe levels in a beam pattern developed using recorded data. Can any one kindly guide me what could be possible causes of higher side lobe levels?

Thanks & Regards

Nauman


texts:
I am trying to implement GCC-PHAT in python.

The approach is similar to the following two links:
link1 and link2

It seems the only difference between GCC-PHAT and normal cross-correlation using FFT is the division by the magnitude.

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import rfft, irfft, fftfreq, fft, ifft

def xcorr_freq(s1,s2):
    pad1 = np.zeros(len(s1))
    pad2 = np.zeros(len(s2))
    s1 = np.hstack([s1,pad1])
    s2 = np.hstack([pad2,s2])
    f_s1 = fft(s1)
    f_s2 = fft(s2)
    f_s2c = np.conj(f_s2)
    f_s = f_s1 * f_s2c
    denom = abs(f_s)
    denom[denom < 1e-6] = 1e-6
    f_s = f_s / denom  # This line is the only difference between GCC-PHAT and normal cross correlation
    return np.abs(ifft(f_s))[1:]

I have checked by commenting out fs = fs / denom The function produces the same result as normal cross correlation for a wide band signal.

Here is a sample test code which shows the GCC-PHAT code above performs worse than normal cross-correlation:

LENG = 500
a = np.array(np.random.rand(LENG))
b = np.array(np.random.rand(LENG))
plt.plot(xcorr_freq(a,b))
plt.figure()
plt.plot(np.correlate(s1,s2,'full'))

Here is the result with GCC-PHAT:

Here is the result with normal cross-correlation:

Since GCC-PHAT should give better cross-correlation performance for wide band signal, I know there is something wrong with my code. Any help is very appreciated!


texts: Check vanishing moment property of db4 wavelets

texts:
I am reading the book "a primer on wavelets and their scientific applications" by James S. Walker. In chapter 3, Property I on page 46 says that if a signal $f$ is linear, then the $k$-level fluctuation value $f\cdot W_m^k$ is zero. I am trying to do numerical experiments to verify that. So I take the signal $f=[1,2,3,4,5,6,7,8]$. The wavelet elements of a Daubechies kind 4 are
$$(0.4830, 0.8365, 0.2241, -0.1294) $$. Our scale matrix $A$ will be like
\begin{array}{cccc} 0.483 & 0.8365 &0.2241 &-0.1294 &0 &0 &0 &0 \\
 0 &0 &0.483 &0.8365 &0.2241 &-0.1294 &0 &0 \\
 0 &0 &0 &0 &0.483 &0.8365 &0.2241 &-0.1294 \\
0.2241 &-0.1294 &0 &0 &0 &0 &0.483 &0.8365
.\end{array}
The first 2-level wavelet basis would be
$$W_1^2=\beta_1*A[1,:]+\beta_2*A[2,:]+\beta_3*A[3,:]+\beta_4*A[4,:]$$
where $\beta_1=-0.1294,\beta_2=-0.2241,\beta_3=0.8365,\beta_4=-0.483$ and $A[i,:]$ is i-th row of $A$. What I expect is
$$f\cdot W_1^2==0$$
but the numerical result shows that
$$f\cdot W_1^2=0.3660.$$
Is there anything wrong in this experiment? Thanks in advance!


texts: What is the transfer function $H_{xy}$ of this flow graph?

texts:
I am trying to get transfer function $H_{xy}$ from this graph. I managed to get the other one but I am struggling with this one. Any help or guidance?

Thank you


texts: Is there a way to filter out low volume content in an audio sample?

texts:
Imagine you are listening to a digital file that contains some random speech.

Now, while you are listening to it, you also hear in the background (spatially) some low volume content, either noise, other speakers, maybe even music, etc.

My question is, is there a way to filter this low volume content so I can just
keep the primary dominant component, which in this case would be speech, and sort of attenuate everything else?

The only thing I can think of is some form of VAD or noise gate, or/and some
spectral processing to do a threshold. But these methods seem too harsh on the data and end up usually filtering part of the speech as well.

Is there any other way of doing this or are the methods above pretty much the only way to do it?


texts:
I am getting started with image/video processing. As a start I have a dataset for which the goal is to detect the quality as the video progresses.

Based on my readings I think any video under test would need to be converted into frames and then the quality needs to be determined based on some values calculated out of the frames.

My question is more on the implementation side as I don't understand what exactly I will have to decode (mathematical values) and what values to deduce in order to determine if the quality is good or bad (scale of 1 to 100 may be).

I came across Scikit-Video library which can be useful as I plan to use Python.

Any suggestions on how I should approach this problem?


texts: Genetic or evolutionary algorithms for filter design?

texts:
For what types of filter design problems might genetic or evolutionary algorithms be useful?

What kinds of genetic or evolutionary algorithms are used for DSP problems?

Edit: I expanded the question to include the larger set of evolutionary algorithms, such as differential evolution.


texts: Phasor notation of a complex number with apparently non-constant radius

texts:
I'm asked to find the phasor of this signal:

$$ x_1[n] = 3\sin\left(\frac \pi7 n\right) + 4j\cos\left(\frac \pi7 n\right),\quad (0\le  n \le 20) \,.$$

The phasor is defined as:

$$ G z_{0}^{n}= Ae^{j\varphi}\cdot r^{n}\cdot e^{j\theta*n} = Ar^{n}[\cos(\theta n+\varphi)+ j\sin(\theta n+\varphi)] $$

But I really don't know how to make the conversion! I mean, one of them has a radius of 3, the other is 4. How can I do this?


texts: What is the impulse response of the system with this difference equation:

texts:
The difference equation is written:

$$y(k)-\frac{5}{12}y(k-1)+\frac{1}{24}y(k-2)=x(k)-\frac{1}{2}x(k-1) $$

The system is linear, invariant and causal.
For $x(k)=δ(k)$, the equation will be:
$$g(k)-\frac{5}{12}g(k-1)+\frac{1}{24}g(k-2)=δ(k)-\frac{1}{2}δ(k-1)$$

I calculated $g(0)=1$, $g(1)= -\frac{1}{12}$, $g(2)= -\frac{11}{144}$ and $g(3)= -\frac{49}{1728}$ but i can't find any relation between them.


texts:
I designed a FIR Kaiser windowed low pass digital filter with filter design app, and had MATLAB generate the code for it and return the filter object:

function Hd = filter_lp_100
%FILTER_LP_100 Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 9.2 and the DSP System Toolbox 9.4.
% Generated on: 04-Aug-2017 17:09:08
% FIR Window Lowpass filter designed using the FIR1 function.
% All frequency values are in Hz.

Fs = 500;                % Sampling Frequency
Fpass = 100;             % Passband Frequency
Fstop = 110;             % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0000000001;    % Stopband Attenuation
flag  = 'scale';         % Sampling Flag

% Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fpass Fstop]/(Fs/2), [1 0], [Dstop Dpass]);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);

But this filter does not attenuate the input signal, it even amplifies it!  I expected that PSD for the output becomes smaller, or at most equal to input PSD(for the pass band) and that's not happening in any point. Please see the plots. Blue is input signal power spectrum density, Red is output PSD, and yellow is filter response.

I tried different filters, none of them is performing. My code is very simple, I copied the part that does filtering and plotting filter response:

%% FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%
coef=filter_lp_100;  % MATLAB filter generated function above
out_signal=filter(coef.numerator,1,in_signal);

%%% plotting filter response. f is physical frequencies, fs is sampling rate:
h=freqz(coef,f,fs);
plot(f,10*log10(abs(h)))

Do you see anything wrong with this?


texts:
I am attempting to implement a chirp z algorithm to handle random sized DFTs, and I can not seem to obtain any meaningful results. I have gone over several write ups and "think" I have a handle on what a chirp transform should do, but I wanted to ask here to double check and make sure my understanding of the process is sound. In pseudo code, here is the process I am applying.

*function Chirpz(input x[n], discrete values such that n = 0, 1, ... N - 1)*
*// To keep numbers in perspective I will assume N = 5.*

M = The minimum power of 2 which is greater than 2N - 1 // In our case of N=5 => M=16.

Pad x[n] with 0s so that it has a total length of M.

// Create y[n] by performing the following:

value = 1 - N

for n = 0 to n = 2N - 1

y[n] = e^i(PI / N * value * value)     // W^((k - n)^2 / 2)
value = value + 1

next n

Pad y[n] with 0s so that it has a total length of M.

// Multiply x[n] by W^(-n^2 / 2):

for n = 0 to n = N
x[n] = x[n] * e^-i(PI / N * n * n)

Perform FFT(x[n]) which has M elements using a radix-2 FFT algorithm
Perform FFT(y[n]) which has M elements using a radix-2 FFT algorithm

// With the results of the FFTs now in x[n] and y[n], multiply each value:
for n = 0 to n = M
x[n] = x[n] * y[n]

Perform IFFT(x[n]) which has M elements using radix-2 FFT algorithm

// Multiply x[n] by W^(-k^2 / 2)
for n = 0 to n = M
x[n] = x[n] * e^-i(PI / N * n * n)

// Now the first N elements in x[n] "should" have the result
X[k] = x[n] for n from 0 to N - 1

return X[k]

But my results don't match anything near what is returned by fftw3 given the same input. In fact my results look to be garbage. As I am pretty new to this, I just want to make sure my understanding of the algorithm is sound.

Thanks.


texts: Why MATLAB only allows cpsd and mscohere between vectors with the same length?

texts:
I am studying the applications of finding coherence in the frequency spectrum between two signals for pattern recognition. I have a recording that is way longer than the sampled signal that I want to compare. However, MATLAB only allows cpsd and mscohere between vectors with the same length. I figured that there must be a good reason for that, but I could not find it.

Also, I was digging through the math of the procedure and I figure that I could do a crosscorrelation between the two (although it would be raw, i.e., not unbiased), and then apply the Fourier transform to obtain the correlated frequencies. Do you think it is sound or make sense to do such a thing?


texts: "Perfect waveform" for cross-correlations?

texts:
I am well versed in cross-correlations, as my masters thesis heavily relied on them for music classification and beat detection. So this question focuses more on the signals that I am running an xcorr on as opposed to the actual cross correlation process.

Here is a simplified version of my setup:
Suppose I have the original signal, which I time shift the signal (add a delay), and add noise/distort the signal. Then I run a cross-correlation on the two signals to find this time delay and get some value of R.

Is there some formulation to come up with a signal that will reduce ambiguity errors and optimize my value of R?
I know it will be some non-repeating, pseudo-random signal, but is there a formal formulation for this? Or should I be looking at other techniques?

Perhaps something like this would be perfect: http://tedxtalks.ted.com/video/TEDxMIAMI-Scott-Rickard-The-Wor. You see any issues with this?

I'm not talking about cross correlation processing techniques to improve the performance. I'm strictly talking about the shape of the waveform


texts: About the order of a difference equation

texts: The question is related to a statement in a dsp book.  I am attaching the image of the statement.  In the attached image, is the equation 1.5 a difference equation of order 1 as stated in the book?  I think it is of order 0 since there is no feedback term i.e. no y(n-1).


texts: How can I generate time-variant Rayleigh channel in MATLAB

texts:
One of known models for channel simulation is called Rayleigh Fading channel, it's usually used to simulate wireless channel in time variant/invariant.

As I know, Rayleigh fading channel can be considered as sum of two Gaussian random variables, we can use the below command to generate that channel:

h = 1/sqrt(2)*[randn(1,M+1) + 1i*randn(1,M+1)];

Where M is channel length and h will be the channel.

But what's about if I want to generate Rayleigh channel in time-variant environment?   is that possible?  Please correct me if the above command generates Rayleigh channel or there is another command.

thank you


texts: State of the art in Cancer tumour detection using Digital image processing algorithms

texts:
I am trying to understand how image processing is applied for detection of cancer tumor using image processing. Towards this direction I am trying to get latest information i.e state of the art in this area. I found one book titled

"Image-Processing Techniques for Tumor Detection"

But this is more than 10 years old and I think lots of new developments must have taken place since then.

I will be thankful if somebody can point me to references which are latest on the subject.


texts: Reflecting negative frequency to positive frequency

texts:
I am trying to synthesize sinusoids by using window function in the frequency domain.

It involves:

In frequency domain, shift the window to center around the peak frequency
To generate a DFT frame, sample a few values of the window around the peak as the spectral motif
Inverse-Fourier transform the spectral motif, to generate the sinusoid in the desired frequency

This approach works great in general, except for synthesising low frequencies. Because when shifting the window to the low frequency, the left side of the window will sit in the negative frequency domain.

The figure in the middle demonstrates the issue (T(k) sits in the negative domain)

I found a solution here, it suggested to add the complex conjugate value of the left tail of the window (in the negative domain) to the DFT bin that’s on the right tail of the window (in the positive domain). Which I can’t make sense of, and following this solution creates even more distortion. So I wonder if anyone knows how to do it properly. Any suggestions would be very appreciated!

Some excerpt from the aforementioned solution:

  The reflection about the k=0 axis is due to the specific embodiment
  described herein for synthesizing a sinusoid. For each real sinusoid,
  one peak exists in the positive frequency bins and another peak exists
  in the negative frequency bins. In the embodiment wherein only the
  peak in the positive frequency bins is synthesized, a peak centered
  about a low positive frequency bin spills into the negative
  frequencies (as shown by the plot for Ht(k−bc) in FIG. 3). Similarly,
  a peak centered about a low negative frequency bin spills into the
  positive frequencies. The portion of Ht(k−bc) in the negative
  frequencies that is reflected, or T*(−k), represents the portion of
  the peak centered about the negative frequency bin that spills into
  the positive frequencies.

PS. Some time ago I've raised this question on DSP-related forum. I've got very detailed suggestion from Robert, to ignore the bins lying on negative domain to specify the whole positive frequency, and then complex conjugate to reflect it to the negative frequency, which has improved the problem but it still can’t go down below 80 Hz. So I thought I’ll post again here.


texts:
I am designing some Downsampling filters using the remez exchange algorithm (Matlab FDAtool). I see that some times i get a spike at the first and last coefficients. does anybody knows why is that? my intuition seems to be that this can introduce some artifacts to the signal..
What will be the consequance of multiplying the filter by a window (hamming and so on) to reduce this edge effects?


texts: Understanding $\mathcal Z$-transforms and pole locations

texts:
I am trying to gain a better understanding of pole locations in the $z$-plane of a given discrete transfer function, $H(z)$.  I think I have a pretty good understanding of how to use the $\mathcal Z$-transforms for filter design, but there is something about the core concepts that still confuse me.

I understand that poles are locations in the $z$-plane (complex plane) where a specific value of $z$ results in $\left|H(z)\right|$ approaching infinity.  This happens when the denominator of a $\mathcal Z$-transform expression has a zero at $z=a$.  i.e.

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

would have one pole at $z = 0.5$.  The inverse $\mathcal Z$-transform of this fraction would be

$$
h(n) = \mathcal{Z}^{-1}\left\{\frac{1}{z-0.5}\right\} = 2\cdot0.5^{n-1}
$$

Since $z$ is less than 1, this impulse response will decay as $n$ increases.  So the response is stable (which happens as long as the pole lies within the unit circle).

However, if we look at the definition of the $\mathcal Z$-transform, with $z$ written in complex exponential form, i.e. $z=re^{j\omega}$:

$$
H(z) = H\left(re^{j\omega}\right) = \sum_{n=-\infty}^{n=\infty} h(n)\left(re^{j\omega}\right)^{-n} = \sum_{n=-\infty}^{n=\infty} h(n)r^{-n}e^{-j\omega n},
$$

we can now express the $\mathcal Z$-transform for our pole directly using using $z=re^{j\omega}$ with $r=0.5$ and $\omega$=0:

$$
H(0.5) = \sum_{n=-\infty}^{n=\infty} h(n)0.5^{-n} = \infty
$$

We know this sum is infinity from our original transfer function.

Now here is where my question begins to form. When $h(n)$ is summed against the $0.5^{-n}$ function, we get an unbounded summation, or a pole.  Of course that can only happen if $h(n)$ is infinite (which it is), since the above summation for $H(0.5)$ could not be infinite otherwise.

But if $h(n)$ is infinite such that it sums to infinity multiplied by $0.5^{-n}$, it would clearly also sum to infinity against $0.58^{-n}$ or $0.25^{-n}$ for that matter.

So my question is, if there is a pole at 0.5, why aren't there poles all along the real axis (except for zero)?  Why is there only a pole at 0.5?

To make it explicit, we can use the specific $h(n)$ that was the inverse $\mathcal Z$-transform of $H(z)$, above:

$$
H\left(\frac 14\right) = \sum_{n=-\infty}^{n=\infty} \left[2\cdot\left(\frac{1}{2}\right)^{n-1}\right]\cdot \left(\frac{1}{4}\right)^{-n}
= \sum_{n=-\infty}^{n=\infty} 4\cdot 2^n = \infty
$$

So, this sum makes it look like there is a pole at $z=0.25$ for the $\mathcal Z$-transform of $h(n)$, but the transfer function says there is not!

In other words, the single pole makes sense to me looking at the $H(z)$ transfer function, but it doesn't make sense to me when I take it and use it directly in the $\mathcal Z$-transform expression.

What am I doing wrong?


texts: Generating good synchronization sequences

texts:
As part of my work I am putting together a simple bursty BPSK communication system.  I want to make the demodulation easy, so I am going to put a synchronization sequence at the beginning of each burst to help the receiver lock onto the frequency, phase, and timing.

In order to get good synchronization information and avoid false positives the autocorrelation of the synchronization sequence should be as close to a delta dirac as possible.  There may also be additional constraints, like starting the sequence with all 1's to help the receiver get frequency lock quickly.

How does one go about generating such a sequence?  Is there a systematic way of doing it, or do you use common sense and trial and error?


texts: RNN as posterior probability estimation in speech recognition with HTK

texts:
I'm new to speech recognition and deep learning and in a learning phase.

I'm trying to follow this paper to learn how to use RNN as posterior probability estimation in an HTK environment. The paper proposes RNN-HMM hybrid system, so for the HMM part I need to use HTK platform.

The problem is that I couldn't even start from anywhere. I have a sample code which uses HMM to recognize digits, but I'm unable to solve at which part should I insert RNN to the code.

If there are any ideas, I would be glad.

The system in the paper is shown as follows:

I have codes in python environment and HMM is applied using HTK.
After converting data to MFCC format, I should use RNN, but after using RNN, which steps of HMM should I apply to generate RNN-HMM acoustic system.


texts: Why do mobile phone calls produce a vibration when they are near a speaker?

texts:
I remember that years ago when I had a classic mobile phone (the small one without touchscreen) when I was about to receive a call from someone and my phone was near a speaker, the speaker would start to vibrate quite hard before 2-3 seconds of receiving the call.

I don't have a way to experiment that now since I don't own a small mobile phone anymore and with modern smartphones it doesn't work. Is there any particular reason why?

How is the phenomenon produced after all? I also tried to google about it, but I had no idea how it is called.

Honestly I have no idea where I should've put this question, but searching trough all SE sites I ended up here, but if this isn't the right place please tell me.


texts: Can we decompose a correlated random variable into correlated plus uncorrelated?

texts:
Assume a random sequence
$
Y_k  = X_k  \cdot H_k  + W_k
$
 where $X_k$ is a deterministic quantity, $H_k$ is a correlated random sequence, $W_k$ is Additive White Gaussian random sequence.

Can we decompose/rewrite the above sequence into a form like below?
$
Y_k  = X_k  \cdot \hat H_k  + X_k  \cdot \tilde H_k  + W_k
$

where $\hat H_k $ is deterministic/white and $\tilde H_k$ is correlated random - something like that?


texts: Finding remaining zeros of linear phase FIR filter

texts:
Question given is: if $h(n)$ is a linear phase causal FIR of order $10$ with real coefficients, find the remaining zeros of this filter if the zeros given are

$$q_{1,2} = -2 \pm 2j$$

$$q_{3} = -\frac12 + \frac{j}{2} $$

$$q_{4} = -1$$

I know for sure that $-\frac {1} 2 - \frac j 2 $ has to be another zero, but there are 5 more zeros unaccounted for. How am I supposed to find them?


texts: Online materials on Support Vector Machine

texts:
I am doing my final year project on Image Processing.I want to know about the Support Vector Machine.Please suggest some links for the former.

Thanks and regards in Advance


texts: How do I filter out this horizon from my contour?

texts:
I am trying to identify the T-rex and obstacles (cactus and bird) from Google's no-internet game.

What I have done is:

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

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

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

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


texts: Spread-spectrum systems and multipath fading

texts:
I didn't find a suitable tag for this, but figured it would fit under "Frequency spectrum".

Which of the Spread-spectrum systems listed here are most resistant to multipath fading? Please, briefly tell me why?

The Wikipedia article above lists these (or a combination of at least 2?):
FHSS: Frequency-hopping spread spectrum.
DSSS: Direct-sequence spread spectrum.
THSS: Time-hopping spread spectrum.
CSS:  Chirp spread spectrum.

Are there other distinct techniques the article doesn't list?

--gv


texts:
I am trying to receive and demodulate a 8PSK modulated sinusoidal waveform (real not complex) which is constant envelop and changing only in phase of the continuous sinusoid carrier every Symbol Period.
Here symbol 0(000) to 7(111) are mapped to sinusoidal phase of 0:45:360 degree respectively

For further explanation let me set some parameters

Ts = 1 Hz //8PSK Symbol rate
Fc = 8 Hz //center frequency

Now at the receiver I am multiplying with the same carrier and sampling the resultant at lets say Fs = 64 Hz.

I want to find out the transmitted symbol which is actually phase of the received signal from these discrete samples in two cases:

1st assuming receiver carrier phase is in sync with transmitter carrier phase

2nd assuming receiver carrier phase has some fix delta phase difference.

I tried low pass filtering the samples and plotting which does give different levels as per the actual phase modulated but difference between two consecutive phase (45 degree apart) is not linear hence direct thresholding based on average level in symbol period is not working correctly.

Any pointers are welcome.


texts:
Here is a sample data that I got from the accelerometer sensor:

The first selected area shows the acceleration of the IMU device to right side and the second shows the movement to left side.
I am looking for an algorithm to detect selected areas (as an array of floats) and provide the direction.
I almost don't know anything in this area (signal processing), So I tried to avoid using incorrect keywords, feel free to edit this post.


texts:
Isnt a Step (as a Transfer function -> $ \frac{1}{s} $) and Impulse (trans fcn -> 1) therefore Impulse + $ \frac{1}{s} $ should be equal to a Step in Simulink right?
I made a Simulink model:

But the outputted plots are not the same ?
Why is that? Am I wrong? Should those two blocks be essentially the same?


texts: Pole-Zero Configuration of Glottal Signal

texts:
I am currently studying for an exam in my Speech Processing course this Friday and I stumbled upon this question in one of the sample exams:

The answer to this question was given as:

Can someone explain to me why the pole is drawn in that location and why the reciprocal was taken?


texts: Why is my time domain interpolation via zero-padding in frequency domain wrong?

texts:
Since the process can be applied in either domain to increase the sampling rate in the other domain, I am trying to apply zero-padding in frequency space to recover a 'cleaner' interpolated signal in temporal space.
To do so, I insert zero-valued frequencies in the spectrum at the location of higher frequencies, which is a common practice.
However I don't seem to recover the original signal very well (in black below) after zero-padding (in red).
import numpy as np
import matplotlib.pyplot as plt

# odd dimension for simplicity
n    = 19
npad = 99

x    = np.linspace(0.,4.*np.pi,n)
xpad = np.linspace(0.,4.*np.pi,npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)//2
f_fwd_pad[0:h+1]   = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[h+1:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

fig, ax = plt.subplots(1,2)

ax[0].plot(x,np.real(f),linestyle=None,marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),linestyle=None,marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

Is that result expected? Is there some fundamental understanding that I am missing?


texts: equations to calculate coefficients of recursive IIR Gaussian blur filter

texts:
I am reading the implementation of IIR Gaussian blur filter from this link.

In the source code, the coefficients are generated by the following code

/* Calcualte Gaussian Blur Filter Coefficiens
 *  alpha -> smooting gradient depends on sigma
 *  k = ((1-e^-alpha)^2)/(1+2*alpha*e^-alpha - e^-2alpha)
 *  a0 = k; a1 = k*(alpha-1)*e^-alpha; a2 = k*(alpha+1)*e^-alpha;
 *  a3 = -k*e^(-2*alpha)
 *  b1 = -2*e^-alpha; b2 = e^(-2*alpha)
 */
 void calGaussianCoeff( float sigma,  float *a0, float *a1, float *a2,
   float *a3, float *b1, float *b2, float *cprev, float *cnext)
{
  float alpha, lamma,  k;
  // defensive check
  if (sigma < 0.5f)
     sigma = 0.5f;

  alpha = (float) exp((0.726)*(0.726)) / sigma;
  lamma = (float)exp(-alpha);
  *b2 = (float)exp(-2*alpha);
  k = (1-lamma)*(1-lamma)/(1+2*alpha*lamma- (*b2));
  *a0 = k;
  *a1 = k*(alpha-1)*lamma;
  *a2 = k*(alpha+1)*lamma;
  *a3 = -k* (*b2);
  *b1 = -2*lamma;
  *cprev = (*a0 + *a1)/(1+ *b1 + *b2);
  *cnext = (*a2 + *a3)/(1+ *b1 + *b2);
}

Does anyone know where the equation is from?


texts: How do ADC/DACs deal with out of sync signals near Nyquist frequency?

texts: Let's say a converter is sampling continuously at 48 kHz and the incoming waveform is a sine wave at 24 kHz... How does it faithfully recreate the waveform if the samples are taken at points that are 0?


texts:
I asked this question recently.  From the comments I learned that squelch is not a hard and fast thing, and that many considerations might go into determining a reasonable value, depending on use, the system, etc.

In my system I've decided to average power levels on the frequency of interest.  I use a moving average, so that average might change dramatically over time.  Somewhat arbitrarily, I decided first to use $3\textrm{ dBm}$ up squelch.  So for example, if the power at the frequency is at an average of $-11\textrm{ dBm}$, I let the signal pass at $-8\textrm{ dBm}$.

I was getting way too much coming through, though (which indicates I don't understand much about the signal variance?)  So I changed it to $6$-up $\textrm{dBm}$.  But, still a lot of garbage goes through.

$3\textrm{ dBm}$ means doubling the signal, right?  Is that not sufficient for something like, say, an analogue voice channel?  is $6\textrm{ dBm}$ up still too little, in typical application?  I really have no idea where ballpark / acceptable ranges are even.  I just need a little direction.
I am considering combining my above technique with a static squelch level (ie. $x\textrm{ dBm}$ above average, but must be above $-x.x\textrm{ dBm}$).  Is that also a typical approach?


texts:
It is common knowledge that the duration time samples are collected is inversely proportional to the width of the main lobe in the Fourier domain. The simple example is the rectangular pulse in time whose Fourier transform is the sinc function. The width of the rectangular pulse is inversely proportional to the position of the first zeros of the corresponding sinc function.
As a result, if you desire to decrease the width of the main lobe in the frequency domain, you should take additional samples in the time domain.
I am working on a project wherein I can collect only a certain number of sequential samples, let's say $N$. After I collect $N$ time domain samples, I must wait a certain duration before I can collect another set of $N$ samples. Let's define the number of missing samples in the "gap" as $M > N$. As a result, the time samples might look something like the following:
There are 3 scenarios to consider:

"Sample sets": two sets of samples over the time indices $n \in [0, N-1] \cup [M + N, M + N - 1]$
"Aggregated Samples": single sets of samples over the time indices $n \in [0, 2N - 1]$
"Fully Sampled": single sets of samples over the time indices $n \in [0, 2N + M - 1]$

Considering the scenario 1, shown in the first figure, I have collected $2N$ samples but there is a spacing between them. My goal is to recover the Fourier transform of either scenario 2 or 3. In either case, we should see an reduction in the main lobe width in the frequency domain. However, because of the large gap between sample sets in scenario 1, if I use a nonuniform fast Fourier transform (NUFFT computed in MATLAB) or fill the gap with zeros (which is mathematically identical to the NUFFT in this case), the resulting grating lobes dominate the spectrum. Below is a comparison of scenario 3 and scenario 1. As evident in the figure, the frequencies that are closely spaced are considerably distorted due to the finite size as compared with the ideal case. A similar phenomenon can be observed if we were to compare scenarios 2 and 1.

If I were to recover the frequency domain representation of scenario 3, I essentially would be interpolating between the sample sets.
If I were to recover the frequency domain representation of scenario 2, I essentially would be moving the second sample set next to the first and resolving the phase difference between the set of samples that would exist at $n \in [N, 2N -1]$ and the set of samples I have which were taken at $n \in [M + N, M + N - 1]$.
Some other notes:

The signal I am sampling consists of a sum of sinusoidal signals each of which have different frequencies and are likely to be close to each other in frequency, but the number of different components is typically large $>> N$.
For various reasons, I cannot use an autoregressive (AR) method, MUSIC/ESPIRIT/etc. method, Prony's method, or matrix pencil algorithm
I have attempted various compressed sensing (CS) methods with limited success as my samples are taken with such a specific structure rather than somewhat randomly spanning the entire time duration as common in many CS problems.
I have attempted with some success a fully convolutional neural network approach to learn to fit the NUFFT of scenario 1 to the FFT of scenario 2 or 3.

I would like to get the community's thoughts on the problem and possible solutions. Is there anything I am missing that I could use to solve this problem?
My specific question however is this. I believe the problem appears to be somewhat nonlinear and a data-driven approach may be the best option. In this case, what would be an optimal way to exploit the domain specific knowledge of signal and scenario towards choosing the best deep learning approach? What are some options, preferably in the form of research papers, for taking advantage of the scenario using deep learning? A sparse CNN? A generative model (maybe GAN) for interpolating the gap?
Thanks!


texts: Power and Energy spectral densities of a linear system

texts:
Here is a schematic of my desired circuit.

I want to calculate Power and Energy that the load dissipates during 0< t<8 You can see my attempts:

Therefore

How can find load's energy in Frequency domain for desired time
domain interval (0 < t < 8) ?
Did I obtain power correctly in
    frequency domain?


texts: Recently, I read a paper. In it, it said "a $32^{nd}$ order 100MHz bandpass filter". So, how the filter was designed? Does it mean this filter have 32 taps?


texts:
I am currently working on a product(cant disclose too much details of the actual project) in which some data is being sent via audible(but on the higher freq) sound that is encoded with dqpsk.

However, it gets very inaccurate at times since the receiver and the broadcaster does not have the same sampling windows. I am not trained academically for DSP and therefore having a lot of problem finding out what are the correct terms to search for.

Basically after a lot of debugging and simulation, i found that the most major issue is that the phase information gets more and more inaccurate as the output and the input gets more out of sync... Are there some techniques to find the beginning of the actual signal in the samples? The data isnt contained in just one sample window(they can span through quite a lot of windows)

Thanks! (I would appreciate any sort of help, even if it is just some term i should use to search on google... )

And on a seperate note, when the timing is relatively accurate, I am already able to get the data properly decoded, so I am not really having much problem with the psk side


texts: FIR filters: direct form / transposed FIR

texts:
I am learning about FIR filters and I'm confused.
I am trying to find out different types of FIR filters.

Is direct form and n-tap FIR filter the same?
What does transposed FIR filter do?


texts: How to find derivative of an image along a radial direction with an arbitrary center of rotation?

texts:
I want to find the derivative of an image along a radial direction. For instance, in the image below I want to find the gradient at a point $P$ of the image, along the direction of the circumference of the circle.

Just like we do image intensity derivatives $\frac{dI}{dx}$ or $\frac{dI}{dy}$, all I want to do is to find $\frac{dI}{dr}$ where $r$ is a variable that stores the rotation angle of the image with a particular center.

From what I know we need to follow chain rule of partial differentiation.
$$\frac{dI}{dr} = \frac{dI}{dx} \cdot \frac{dx}{dr}$$

We know $\frac{dI}{dx}$ but how to find $\frac{dx}{dr}$?


texts: Signal Processing in Amateur Radio Astronomy

texts:
I am an amateur astronomer and I took recent interest in Radio astronomy and I am also new to Digital Signal Processing. My team is trying to build a radio telescope. I am interested in the data analysis part. One thing that is predominant in the signals that is received is the noise.

This is the data from the Radio Telescope of The University of Illinois recorded on July 4th 2011 and the telescope was pointing towards the Cygnus constellation:

http://radio.uindy.edu/data/2011/07041059-2.txt

Note: I am interested only in the third column of data

I made a plot of the data using R and that turns out to be something like this:

This is the plot of the Cygnus constellation inferred by the University:

I manually matched the above image with my plot and took out a portion of my plot:

It seemed similar to me. The chances are that my manual matching is wrong. How do I tackle with the situation ?

I have just started learning DSP and please suggest resources and pointers for DSP techniques so that I may apply in radio astronomy.


texts: Radar signal processing - Difference between range gate/range bin

texts:
In learning about radar signal processing, I see two different terms that sometimes seem to refer to the same thing, but sometimes different things:

Range gate
Range bin

Usually, I interpret both to simply mean "time sample", but sometimes I see range gating used to describe just taking a few samples around the estimated target range during a tracking phase.

Is there a more nuanced definition that I should be aware of?


texts: Sampling from a distribution with arithmetic coding

texts:
I'm struggling to implement a sampling scheme where arithmetic decoding is used to generate the samples.

This is based on excercise 6.3 from David MacKay's book, where you compare rejection sampling with arithmetic decoding sampling. For each symbol you want to generate with rejection sampling you have to generate an (uniformly distributed) integer from $1$ to $2^{32}$, rescale it to be in $\left(0,1\right)$ and then if your random number is less than 0.99 emit a 0 (or 1 accordingly).

To generate symbols by arithmetic coding, you feed standard random bits into an arithmetic decoder equipped with the appropriate model.

(Note: I think this works because when we compress data we go from symbols generated under P to symbols under Q - P,Q are probability distributions - and decompression works the opposite way. If we choose Q to be an easy distribution to sample from then by decoding symbols from Q we can generate symbols from P.)

I've implemented the methods in the code below:

clear all;
pdist = [0.99 0.01] ;
counts = [99 1];
N=100;
U = randi(2^32, [1 N]);
U = U/2^32;
rejection_samples(1:N) = 0;

for i=1:N
    if U(i) < pdist(1) ;
       rejection_samples(i) = 0;
    else
       rejection_samples(i) = 1;
    end
end

seq = randsrc(1, N, [0 1;0.5 0.5]);
arith_samples = arithdeco(seq, counts, N);

I'm confused becuase the arith_samples variable sometimes contains a 2 (but I guess this just maps onto (0,1) as the output never seems to contain a 0 - just 1 and 2).

The answers say that the rejection sampling method takes 32 bits per generated symbol, but that arithmetic coding only takes 0.081 = $H_2\left(0.01\right)$. I can't really see that from my code.

Is there some way I could measure 'bits used'?

Also this seems a bit too good to be true - surely there are better ways of generating random numbers than just rejection sampling?

The relevant chapter is found here:

http://www.cs.toronto.edu/~mackay/itprnn/ps/105.124.pdf


texts:
In IEEE 802.11a standards, there are many fields. Of which is the Data field which has 52 available sub-carriers used for transmission of coded bits. If the subcarrier locations are numbered from 32,-31,...0,1,...31.

I have read in the book that

that 48 out of 52 subcarriers are used for transmission of data bits. While 4 pilots are used as pilots for phase and frequency tracking and training. The pilot subcarriers are located at -21,-7,7,21 with values 1 (at -21,-7,7) and -1 at 21.

I am just wondering, usually the preamble (fields before the data field) is used to estimate the channel by sending a pilot symbol (using ALL subcarriers). Then my questions are

1) what is the purpose of sending these 4 pilot tones?

2) I am guessing phase and frequency tracking is different from channel estimations?

3) My last question is do we send the pilot values on the predefined pilot subcarriers ALONG with the data bits? so basically is each OFDM symbol made up of pilot and data tones? This is different than the procedure in preamble where  ALL subcarriers are used to send one pilot symbol, right?

Thanks


texts: How do I convert a two-pole two-zero transfer function from the s-domain to the z-domain?

texts:
I'm trying to convert the following IIR transfer function from the s-domain to the z-domain:
$$
H(s) = \displaystyle\frac{\frac{s^2}{\omega_z^2} + 2\zeta_z\frac s\omega_z + 1}{\frac{s^2}{\omega_p^2} + 2\zeta_p\frac{s}{\omega_p} + 1}
$$

where $\omega_z$, $\zeta_z$, $\omega_p$ and $\zeta_p$ are constants but not known ahead of time (ie. I can't use the actual values when deriving $H(z)$).

I tried using partial fraction decomposition to convert the function to a sum of single-pole filters, with the intention of substituting $s + n$ with $1 - e^{-nt}z^{-1}$ (where $n$ comes from the single pole filter and $t$ is the sample period). This ended up with an unmanageable expression that I couldn't simplify with Maxima.

Is there a more straightforward approach I can take?


texts: Expressing 2D Gabor transform as 1D spectrum

texts:
The Gabor transform (Short-Time Fourier Transform with Gaussian window) is represented as $U(t,\omega)$, where $t$ is the time, and $\omega$ is the frequency.  The author of a monograph (Seismic Inverse Q filtering, pg. 190) is asking for the Gabor spectrum to be transformed from 2D to 1D in the following fashion.

Define $\chi = \omega t$
Take the 2D amplitude spectrum $|U(t,\omega)|$ into a 1D amplitude spectrum as $|U(\chi)|$

The issue that I am struggling with is that $\chi$ is the product of $\omega$ and $t$.

The following Matlab code demonstrates that $\chi$ does not continue to increase monotonically.  Note the "sawtooth" appearance of $\chi$.

So how do I take the 2D signal as a 1D signal?  For each $\chi$, there should be exactly one $|U(\chi)|$.

N0 = 51;
N1 = 1604;
fs = 1.0e3;
f = (fs/2) * linspace(0,1,N0);
omega = 2 * pi .* f;
dt = 1 / fs;
T = (N1-1) * dt;
t = linspace(0,T,N1);

chi = [];
cnt = 1;
for i = 1:length(omega)
   for j = 1:length(t)
        chi(cnt) = omega(i) * t(j);
        cnt = cnt + 1;
   end
end

figure;
plot(chi);
xlabel('\chi')
ylabel('Value')

However, the Gabor spectrum is a 2D matrix, similar to this example, where $\omega$ is on one axis, and $t$ is on the other axis:


texts: How to include phase in a sinusoidal Kalman Filter

texts:
I start with the equation for sinusoidal motion with an offset and differentiate to get the 2nd order ODE describing the motion of the object.

\begin{align}
x &= A\sin(\omega t + \phi) + O\\
\dot{x} &= A\omega \cos(\omega t +\phi)\\
\ddot{x} &= -\omega^2 A \sin(\omega t +\phi)\\
\implies \ddot{x} &= -\omega^2(x-O) = -\omega^2x + \omega^2O
\end{align}

In matrix form this becomes:

$$
\begin{bmatrix}
\dot{x} \\
\ddot{x} \\
\dot{O}
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 & 0\\
-\omega^2 & 0 & \omega^2 \\
0 & 0 & 0 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
x \\
\dot{x} \\
O
\end{bmatrix}
$$

$\implies$ the system dynamics matrix $\mathbf F$ is:

$$
\begin{bmatrix}
0 & 1 & 0\\
-\omega^2 & 0 & \omega^2 \\
0 & 0 & 0 \\
\end{bmatrix}
$$

then by performing $$
A(t) = \mathscr{\mathbf L}^{-1}\left(s\vec{\mathbf I} - \vec{\mathbf F}\right)^{-1}
$$

I get the following discretised $\mathbf A$ matrix:

$$
\begin{bmatrix}
\cos(T_s \omega) & \dfrac{\sin(T_s \omega)}{\omega} & -\cos(T_s \omega) + 1\\
-\omega \sin(T_s \omega) & \cos(T_s \omega) & \omega \sin(T_s \omega) \\
0 & 0 & 1 \\
\end{bmatrix}
$$

This works to filter a sine wave of a know frequency with an offset, it predicts the correct phase, amplitude and even if a change the amplitude of the signal being filtered half way in it corrects for this very quickly and fits to the true signal. However if I shift the phase of the true signal part way through it does not correct for it. How might I include phase $\phi$ in my Kalman filter model such that it can still track the signal with sudden changes in phase?

If I try with my current filter I get the following result:


texts: In edge detection, what value is being compared to the threshold value?

texts:
When we say we set a threshold to the edge finding, what values are we exactly thresholding on? I know a threshold value to limit only certain edges to show up, but what value of the edges is compared with the threshold?

In both the sobel and canny edge detection, they are comparing with a threshold value but I am never sure what values of the edges are being compared.


texts: I was thinking, if a realtime system is one in which the processing time is shorter than the sampling period, how does a feedback filter such as a biquad perform the computation without waiting for, say, in the case of the biquad, two additional samples?


texts: How to modify an LTI differential equation to be acausal or anti-causal?

texts:
I'm trying to wrap my head around causality in LTI-systems. Considering continuous time only, I'm happy with the fact that the system is causal iff the impulse response function $h(t)=0$ for $t<0$.

What I don't understand is how this is embodied in an LTI differential equation. Specifically, what makes the positive and negative time directions fundamentally different? (a point which is especially clear if considering e.g. length as the independent variable, rather than time)

I'm hoping to shed some light on this by considering the LTI equation

$ \frac{dy}{dt} = u(t)$

For this, or a more complex LTI system, I could set the RHS to $\delta(t)$ and then use Laplace transforms to solve for $y(t)=h(t)$. In this specific case $h(t)$ is the Heavyside step-function, i.e. $h(t)=0$ for $t<0$ and $h(t)=1$ for $t>0$. Thus this is a causal system.

If I wanted to modify the above differential equation to make the system either acausal, or anti-causal, what would I have to do?


texts: How we can represent a bandpass signal equation in matlab

texts:
The general form of bandpass signal is given by
$$s(t) = \Re\{S_b(t)\exp(j2 \pi f_ct)\}$$
How this can be represented in matlab.


texts: Why there is static noise and voice distortion after convolving it with impulse response?

texts:
I have been working on the decorrelation of audio signals while following this paper :
The Decorrelation of Audio Signals and Its Impact on Spatial Imagery
So far I have generated impulse responses in such a way that magnitude is 1 in all frequency spectrum, while random phase generation between -π to π. I correlated these signals and I am getting the required correlation.
Now moving on, when I convolve this impulse response with my audio signal, results in the case of decorrelation are good but the audio becomes noisy and distorted. I am attaching the reference wav file and decorrelated wav file.
Some points regarding computation:

All audio data computations are in float range [-1.0:1.0].
I am using Python to perform all the calculations
The input reference file is in Short (PCM 16 bit) format which is converted to float range [-1.0:1.0].

I have saved convolved_decorrelated_output.wav in both formats PCM 16 bit short and Float as mentioned above. Still getting the same results. So there are no issues regarding data type conversion as per my knowledge.
Any help will be appreciated.
Regards,
Khubaib Ahmad
P.S: Link to download wav files:
https://www.researchgate.net/post/Why_there_is_static_noise_and_voice_distortion_after_convolving_it_with_impulse_response
Impulse Response:

Convolution_code:
resultant_chunk = list()
for chunk in data_chunks:
    resultant_chunk.append(np.convolve(chunk,impulse1,mode="same"))

resultant_chunk = np.hstack(resultant_chunk)
# diff = len(resultant_chunk) - len(data)|
diff = len(resultant_chunk) - len(recorded_data)
resultant_chunk = resultant_chunk[:-diff]
len(resultant_chunk)


texts: How to estimate the mean of a clipped signal?

texts:
An accelerometer is used on a high vibration environment (a bicycle) to measure  gravity. The signal clips, so data higher than a clipping value is clamped to that value.
So, simply filtering the signal gives the wrong median mean value when the vibration is high.

I want the mean value over a time much longer than the vibration frequency, so a few second filter is fine. I do not care about the vibration waveform itself. Just want the average value. How to remove the error caused by the clipping?

Bad idea one:
I had suggested that the signal's noise distribution be modeled, the amount of clipping predicted based on measured vibration amplitude, and that used to correct the data. Maybe the distribution is Rayliegh, maybe Gaussian, maybe... So maybe that first idea is not a good one.

Possibly brilliant idea two:
But how about just measuring the distribution on the spot? Make a histogram with several bins, accumulating a count in a bin each time the sampled data is in a bin's range.

Then, fit a curve over the top 3 or 5 highest bins, and the position of the peak is the signal median the modal mean, where the signal mean value would be, if it were not clipped.

Is this crazy? Good? Done before and has a name?


texts: What is the principle behind my voice echoing perpetually when 2 separate devices are used for playback and recording on video calls?

texts:
I had to use 2 separate devices for a video call, desktop for video/audio output and smartphone for recording, since the desktop microphone was not working properly. Unfortunately for me, when I started speaking, my first few words got echoed perpetually, and my following words were ignored by the microphone. What could be the reason behind this phenomenon?
More importantly,
Why does this phenomenon not occur in video calls if the same device is used for recording and playback? How is this phenomenon avoided?


texts: Convolve a minimum phase wavelet with delta peak

texts:
let's assume I have two time series: one with a delta peak at position $x$ and another time series with a minimum phase wavelet that starts at time $y$ with a maximum peak at position $z$.

Now, when I convolve both time series, is the beginning of my minimum-phase wavelet at position $x$ or is the maximum amplitude of the minimum-phase wavelet at position $x$ ?

I know, that when I convolve a zero-phase wavelet, the maximum amplitude is at the position of the delta function.

I am puzzled..........


texts: Using the flip-and-drag method to draw a convolution

texts:
So, we are asked to graphically convolve (using the flip-and-drag method) the following pair of signals: $$x_1(t) = \delta(t+1)-2\delta(t-1) \ \ \ \ \ \ \ \ \ \ \ \ \ x_2(t) =  \begin{cases} -1-2t, & -\frac{1}{2}\leq t\leq 0 \\ -1+2t, & 0\leq t\leq \frac{1}{2} \\ 0 & |t|>\frac{1}{2}\end{cases}$$

where $\delta(t)$ is the dirac delta function.

I don't know how to draw using Latex, so I am going to try and describe you what I've drawn so far using words (as best as I can). So, I've successfully been able to draw $x_2(t)$ (which gives me an inverted triangle whose summits coordinates are : $\left(-\frac{1}{2}, 0\right), \left(\frac{1}{2}, 0 \right)$ and $\left(0, -1 \right)$).

My problem starts when I have to draw $x_1(t)$. After doing some research, I have found that we represent the dirac delta function using "arrows", but I still don't know how to draw $x_1(t)$ .

Also, if we assume I was able to draw it (i.e I would have a graph with arrows), how am I supposed to draw the convolution (using the flip and drag method)?


texts: What is the betwwen [i,j] and [u,v] in explinations of correlation and convolution (picture attached)?

texts:
I don't know how to write 'summation' symbol here hence posting the picture.

Can someone explain to me the difference between i,j and u,v in this explanation of correlation and convolution?

I know this discussion must have taken place somewhere but it is difficult to look for such discussion by keyword 'correlation' and convolution'.

I have already referred to the following posts:
1. Difference between correlation and convolution on an image?
2. What is the difference between convolution and cross-correlation?

Adding the source for reference: https://www.youtube.com/watch?v=C3EEy8adxvc
Thank you in advance.


texts: Amplitude decrease in FFT at higher frequency for multisine

texts:
I am adding sine signals with different frequencies but same amplitude=100.
Sampling frequency is  50MHz.
I am not getting 100 in frequency domain for frequency 787KHz.
At 1K, 3K I am getting 100 amplitude,
but at 787KHz I am getting 75.9 .
All the signals have 100 V amplitude.
Any reason for this?


texts: Getting a single feedback peak from a basic 4 stages phaser

texts:
I try to set a very basic 4 stages phaser (two 2-pole all-pass filters), with a particular feedback path where only one peak is raised at center frequency when the feedback is above 0.

The following amplitude responses describe what I try to get: in the first measure (light green) we get two notches (no feedback). In the second measure (purple) we get two notches and a single peak at center frequency (feedback).

My chain is very simple: the input is sent into two 2-pole all-pass filters having the same center frequency, the output is then mixed 50/50 with the input. I've also added a feedback path, where the output of the last all-pass filter is feed back to the input of the first one (input + state * feedbackAmount).

Instead of getting a single peak I get two (as described in the following picture). I've also tested negative feedback it does not help. Any idea?

Edit1:
My signal flow leading to two peaks when the feedback is close to 1:

Edit2:
My implementation may be wrong indeed, to be sure I've also run a test in Ableton Live 8. Noise straight to the phaser.

First: no feedback:

Second: full feedback, we get two peaks too. (Most phaser I've tested provide the same result):


texts:
I'm down sampling voice audio by first performing an FFT, then only taking the parts of the result that I need, and then performing an inverse FFT. However, it's only working properly when I'm using frequencies that are both power of two, say down-sampling from 32768 to 8192. I perform an FFT on the 32k data, discard the top 3/4 of the data and then perform an inverse FFT on the remaining 1/4.

However, whenever I try to do this with data that doesn't line up properly one of two things happen: The math library I'm (Aforge.Math) using throws a fit, because my samples are not a power of two. If I try to zero-pad the samples so they become power of twos, it get gibberish out on the other end. I also tried to use a DFT instead, but it ends up being insanely slow (this needs to be done in real time).

How would I go about to zero pad the FFT data properly, both on the initial FFT and the inverse FFT at the end? Assuming I have a sample at 44.1khz that needs to get to 16khz, I currently try something like this, the sample being 1000 in size.

Pad input data to 1024 at the end
Perform FFT
Read the first 512 items into an array (I only need the first 362, but need ^2)
Perform inverse FFT
Read the first 362 items into the audio play buffer

From this, i get garbage out at the end. Doing the same thing but without having to pad at step 1 and 3 due to the samples already being ^2, gives a correct result.


texts: I'm having a difficult time understanding why the matrix for LTI systems is a Toeplitz matrix. I can see why $h_{n,m} = h_{n' + q,m' + q}$ given that $n' = n - q$ and $m' = m - q$, and $$\sum_{m'= -\infty}^\infty h_{n' + q,m' + q} x[m'] = y[n']$$ and $$\sum_{m= -\infty}^\infty h_{n,m} x[m] = y[n],$$ but the leap to the equivalence $h_{n,m} = h_{n + q,m + q}$ escapes me.


texts: Is there a way to further optimise my FFT?

texts:
I've written my own FFT but wanted to know if there was perhaps a way to combine the Cooley and Turkey methods. So far I've got N/2 log N real multiplies for transforming a real signal.

My current implementation has N log N multiplies (well a little more because I've not optimised out the case of the complex sinusoid at 0 radians), but I have an optimisation that reduces the number of multiplies by 50%.

Note: this is absolutely not homework, I've written this algorithm because I found the explanations of the Cooley/Turkey method far too convoluted, and at the same time I wanted the ability to select which bins I processed in both the FFT and IFFT as well as being able to output the IFFT to multiple channels. So I just decided to figure out my own radix-2 while at the same time understanding how and why the transform works (which I have).

ORIGINAL QUESTION (edited for clarity) ---

Example code: kfft.c

I wanted to know whether there is a similar FFT algorithm to mine.

I basically perform FFT by factorising calculations by recursively subtracting the second half of a signal from the first, meaning x[n] only needs to be multiplied by the complex sinusoid for only half of the duration, resulting in N/2 log N real multiplies (shown below). The summed signal is for multiplying with complex sinusoids that are multiples of the sampling period. For example groups {2, 6} are multiples of 2, and therefore have a periodicity of 4 samples.

// For bins { 1, 3, 5 and 7 }

x1[0] = x0[0] - x0[4]
x1[1] = x0[1] - x0[5]
x1[2] = x0[2] - x0[6]
x1[3] = x0[3] - x0[7]
// Calculate sum (required for the next level of bins)
x0[0] = x0[0] + x0[4]
x0[1] = x0[1] + x0[5]
x0[2] = x0[2] + x0[6]
x0[3] = x0[3] + x0[7]

// For bins { 2 and 6 }

x2[0] = x0[0] - x0[2]
x2[1] = x0[1] - x0[3]
// Calculate sum (required for the next level of bins)
x0[2] = x0[0] + x0[2]
x0[3] = x0[1] + x0[3]

// For bin { 4 }
x3[0] = x0[0] - x0[1]
x3[1] = x0[0] + x0[1] // <-- This is the 0hz bin (sum)


texts: How to narrow the frequency band of a wavelet

texts:
I have a ricker wavelet with a dominant frequency of 15 Hz. The fourier transform shows its frequency band is almost to 50 Hz. How can I narrow the frequency band of this 15 Hz ricker wavelet?

I have added a picture of the wavelet and its fftshifted frequency band (The x axis is not the frequency samples).


texts: Basics of waveform, how to define the zero point

texts:
Audio is a sequence of waveforms which have peaks, troughs, 0, just like below.

How to define the x axis (y = 0) in the picture, how to choose the 0 point in the waveform since there are so many points on the waveform?

Why not move the waveform above or below the x axis(y=0)?

This question seems a little stupid. But I still hope someone can help with it.


texts: Given the 4 corners of two QR codes on a photographed sheet of paper, can we get a deskewed rectangular image?

texts:
I made a photo of 2 sheets of paper with a phone camera, each of them cropped at the:

(x1, y1) top left corner of top QR code
(x2, y2) bottom right corner of bottom QR code

I hoped that the x, y coordinate of checkboxes would be nearly the same. In fact it's not, because of the angle of the camera, skewness, etc.

Question: How to have a better deskewing by using 8 corners:

the 4 corners of top QR code  (detected automatically by this method): (x1, y1), ..., (x4, y4)
the 4 corners of bottom QR code: (x5, y5), ..., (x8, y8)

?

Note: links of the 2 images: https://i.stack.imgur.com/SR8hz.jpg and https://i.stack.imgur.com/XsxrT.jpg


texts: Help in understanding the formula of Signal-to-Noise-Ratio (SNR) - Part 1

texts:
Question 1:

Consider an Autoregressive model :
\begin{align}
y[n] &= y[n-1] + x[n]\\
z[n] &= y[n] + w[n]
\end{align}
where $y$ is the output observation, $x$ is a random input and $w$ is Additive White Gaussian Noise (AWGN). In general, I have seen in many research articles and tutorials that the SNR = $\frac{E[y^2]}{\sigma^2_w}$

How to arrive at this formula and what is the meaning of Expectation? Is it calculating the mean?


texts: Does it really take 6 minutes and 20 seconds to properly fade-out a single precision audio signal?

texts:
Update

This question is purely theoretical and should not account for any psychoacoustic effects, or practical implementations.
Perhaps the question should be rephrased to "How long will it take to apply the smoothest possible fade, in a single precision system"?

A sine wave is generated within a 32-bit float system, with a sample rate of 44.1 kHz. The task in hand is to fade out this sine with as little digital distortion as possible (this is the meaning of 'properly'), without using dither. The maths below resolves to approximately 6 minutes and 20 seconds to do so. Is this correct?
Axioms

Digital audio within a 32-bit float system is represented with the values -1.0 to 1.0
All the values within this range are normalized
The sine wave peaks at 1.0 and -1.0
The audio resolution is $2^{25}$, made of 23 bit mantissa + 1 implied (normalized) bit + 1 sign bit.
The possible gain coefficients are $2^{24}$, with values ranging between 0.0 to 1.0

Fading out
Any fade will create distortion due to quantisation of the gain values.
For instance, a fade from unity gain (1.0) to silence (0.0), if happens between two consecutive samples, is similar to 1-bit reduction (this will generate high and low frequency distortion correlated to the phase of the sine). If the fade is performed over 4 samples, there will be 4 gain steps; 8 samples mean 8 gain steps; and so forth.
It follows that the smoothest fade possible is that where between each sample the gain changes by the minimum possible amount, which is $1/2^{24}$. This gives $2^{24}$ gain steps.
$2^{24} = 16,777,216$
Meaning 16,777,216 samples are needed for the fade. If the system sample rate is 44,100 Hz:
$16,777,216 / 44,100 = ~380.4$ seconds.
Which equals approximately 6 minutes and 20 seconds.
This seems unreasonable. So where is the calculation wrong?

Just to further demonstrate what distortion I am talking about.
The original sine and its frequency analysis:

With a long fade:

With a medium fade (50 samples):

And with a 1-sample fade:


texts:
In general any form of enropy is defined as uncertainty or randomness. In a noisy environment, with increase in noise, I believe that entropy increases since we are more uncertain about the information content of the desired signal. What is the relationship between entropy and SNR? With increase in signal to noise ration, the noise power decreases but this does not imply that the information content of the signal increases!!
The information content may remain same, so does that mean entropy is unaffected?


texts: Discrete Fourier transform in a multidimensional space

texts:
I want to measure the frequencies at which a point oscillates in a multidimensional space, let's take the example of a point on a 2d-surface.
For now, I naïvely split the signal in two, along the cartesian coordinates, and compute the two FFTs separately. This approach could work but seems to be arbitrary and may loose relevant information. Is there a canonical way to perform the frequency analysis in multidimensional space ?

It is unrelated to 2D FFT. I take a vector of size 2*T and 2D FFT is concerned by plain T^2 data, like images.

Here is an example where halving the signal along the polar coordinates seems better than along the cartesian coordinates.

The data used, a point evolving in 2D,

And here is the FFTs for

The cartesian coordinates x and y (a bit noisy)
The polar coordinates theta and rho (very clear)

edit after Mohammad's comment about velocity

And here is the FFT for the velocity vector (very noisy, and looses information in high frequencies)


texts: Why is circular convolution used in DSP? Why not linear convolution?

texts:
Why are we using circular convolution in DSP?
What's the main solid reason for the use of it in digital
processing?
Why does the concept of circular convolution come more often than
linear convolution?


texts:
I am doing offline signature verification for my project. While researching online I found that tha 1D log Gabor method can be used to find Hamming distance of an image. I got stuck in understanding 1D log Gabor. Wikipedia have some information. My question is how can we find the value of $f$, $f_0$ and $\sigma$ of an image given in Wikipedia 1D formula ?

Here is the link where I got idea about log Gabor can be one method for verification.


texts: Psudocode implementation of true envelope

texts:
I am trying to implement the algorithm described in this paper, I'll quote the relevant portion: http://recherche.ircam.fr/equipes/analyse-synthese/roebel/paper/trueenv_dafx2005.pdf

  Let $V_i(k)$ be the cepstral representation of the spectral envelope at iteration $i$, that is the Fourier transform of the filtered cepstrum, and further initialize the iteration using $A_0(k) = \log(|X(k)|)$ and $V_0(k) = -\infty$. The algorithm then iteratively replaces the current target amplitude spectrum according to
  $$A_i(k) = \max(A_{i-1}(k), V_{i-1}(k))$$
  and iteratively applies the cepstral filter to the updated target spectrum $A_i$.

Could someone write a simple pseudocode version of this algorithm?

I think I'm getting confused by exactly what things are named. If I understood correctly, the cepstrum is the inverse fourier transform of the logarithm of the magnitude of the fourier transform of a signal, and, again, if I understood correctly, I was able to implement a smoothed spectrum by zeroing out the samples to the right of the cepstrum, and then taking the fft and exponentiating. However, when I tried to implement this "true envelope" algorithm I was completely lost.

EDIT: as suggested, this is what I think I should do.

v = array of -9999999
a = log(mag(fft(x))
repeat a couple of times:
    for k from 0 to a.size():
        a[k] = max(a[k], v[k])

    //Maybe? Not really sure.
    v = ifft(log(mag(fft(a)))
    v[from 50 to v.size()] = 0
    v = fft(v)
    //Do I also have to do something to a?


texts:
I have a question concerning methods to find different aspects of "similarity" of signals. So, first, I'd like to characterize the signals to be dealt with:

We have a number of close to rectangular signals between values 0 and 1, starting and ending with some distinct flanks; furthermore, the signals are subject to some noise. Some of the signals are more like peaks which would not reach 1 (2 flanks meeting in the middle).

We have two sets X and Y of such signals containing n vectors each. The vectors in Y are ideally identical to the vectors in X, but time shifted. This is for the ideal case though, in the real application, the signals are not likely to be identical, but we suppose to find some similarities between them.

My question is now: how could i find the time shift between the two sets? First, I need a method to somehow estimate the similarity of the signals in sets X and Y; afterwards, one can find the time-shift.

I thought about cross-correlation first. The issue is, that those peak signals deliver bigger cross correlation coefficients when correlated with long rectangular signals. So, could you think of solution by normalizing the cross correlation coefficients with some other characteristic value of the signal (I didn't find a meaningful one yet).

After that, we thought that probably dynamic time warping techniques could help - I'm still familiarizing with that one. Do you think this could be a good way to go?

Generally: do you have another technique that could be used and helpful for this problem? I'm not an expert in signal processing (I'm just a mechanical engineer), so don't be too cruel to me.


texts: Feature maps for a Convolutional Neural Network

texts:
I hope this is the right place to ask this, so here goes:

I am currently trying to implement a convolutional neural network in C++, but since I have no formal education in signal processing, image classification, or anything of the sort, I am a bit confused as to the filter part of it. So I understand the general structure, but I don't understand where you get the convolution filters from. Do you just use a general set of filters that are supposed to identify things like edges. Or do you train the filters?

So say I have 64x64 images of either a face or not a face. If I understand  it correctly, you pass this image into a set of convolutional filters, and pass that output into something that downsizes the new images. Then you keep doing this until the images are suitable for direct classification. Am I getting something wrong here? And how do I know which convolutional filters to use? Is there some sort of standard set of filters that detect features?

Any help would be appreciated!!!


texts: Shazam / SoundHound Music Analysis - Recommendation for Reading

texts:
I need to explore some area around music frequency analysis.
A good example is the way Shazam / SoundHound app works (mapping frequency domain, comparing, creating some vector map, handling noises and errors).

I am looking for a while for a practical book or articles about this. By practical i mean with a simple language/code examples, etc.

No Sigmas, and integrals for infinity that takes you back to university.

It's really hard to find any music/audio analysis practical material.