texts: Does the integrator in a PID controller removes the disturbance?

texts:
I am reading "Small Unmanned Aircraft: Theory and Practice" (page 7). We have the following Transfer Function (TF):

$\phi(s)=\frac{a_{\phi 2}}{s(s+a_{\phi 1})}(\delta_a(s)+\frac{1}{a_{\phi 2}}d_{\phi 2}) $

where $\delta_a$ is the output signal to the actuator and $d_{\phi 2}(s)$ is an unwanted disturbance.

Initially, the disturbance is neglected, and only a PD controller is employed, and the resulting 2nd order TF (with $k_p$ and $k_d$)is compared with the general second-order TF:

$H(s)=\frac{\omega_n^2}{s^2+2\zeta \omega_ns+\omega_n^2}$

From which the value of $k_d$ is obtained (and $k_p=\frac{\delta_{a,max}}{e_{max}}$, where $e_{max}$ is the maximum error).

Then, the disturbance is considered, and an integrator is added to get a full PID. The resulting TF (with $k_p$,$k_d$ and $k_i$) will be different, actually it will become 3rd order. However, the book does not update the values of the gains after the integrator has been added, and somehow claims that the integrator will basically nullify the disturbance, so no change is needed for the gains.

Can anyone explain this to me?


texts: Causal systems - Paradox? Confusion of definition?

texts:
For a system such that $y[n]=\cos(x[n])$, I'd like to determine whether the system is causal or not.

Definition 1: A system is causal if and only if the output $y[n]$ is a function of a linear combination of inputs $x[n-k]$ such that $k \ge0 $.

We can see that the system is causal by this definition because $y[n]$ only ever depends on $x[n]$.
Definition 2: A system is causal if and only if the impulse response $h[n]=0$ for all $n<0$.

The response of the system to $\delta[n]$ is $h[n]=\cos(\delta[n])$.

$$h[-1]=\cos(0)=1 \ne 0 \quad \therefore\text{The system is not causal.}$$

Have I hit a paradox? Am I missing something here?


texts: Expression to represent alternating pulse train?

texts:
I have a question regarding an alternating pulse train signal, $p(t)$, and a system which looks like this:

Immediately, I'm supposed to be able to extrapolate that
$$
P_t = P_1(t)-P_1(t-\Delta)
$$

How do you obtain this?
How would you get the Fourier transform of this pulse train from there?


texts: How to find the phase lag between a sinusoidal input force and the response?

texts: I have a system that is forced by a sinusoidal function with specific amplitude, phase and frequency. I get a deformation of sinusoidal form as well, and I would like to find the phase lag. I use FFT and isolate the phase lag numerically, it says it is pi/2 however when I look at the plots, it is not pi/2 at all. What can I do? Thank you. PS:Can provide FFT code if needed.


texts: Unable to stabilize 8-PSK signal with self designed costas loop

texts:
I am trying to make a signal processing block which implements the costas loop to stabilize an 8-PSK signal. Here is the flowgraph with "costas_loop" being my block.

My workflow is the following -

in0 = input_items[0]
out = output_items[0]
# <+signal processing here+>
feedback = np.ones(in0.shape, dtype=np.complex64)
on_first_mul = np.ones(in0.shape, dtype=np.complex64)
op_thresh_imag = np.ones(in0.shape, dtype=np.float64)
op_thresh_real = np.ones(in0.shape, dtype=np.float64)
a = np.array([op_thresh_imag])
b = np.array([op_thresh_real])
in_iir = np.zeros(in0.shape, dtype=np.float64)
out_iir = np.ones(in0.shape, dtype=np.float64)
out_vco = np.ones(in0.shape, dtype=np.complex64)
in0 = in0/math.sqrt(2)

for i in xrange(0,self.iter):

  # Multiply input signal with feedback of that iteration
  on_first_mul = in0*feedback
  # Phase error detector
  real = on_first_mul.real
  imag = on_first_mul.imag
  a = self.threshold(imag)
  b = self.threshold(real)
  in_iir = np.arcsin(imag*b - real*a)

  # IIR filter implementing (1.0001 - z)/(1 - z) (old style)
  # self.prev_input[i] contains last input value from previous in0 chunk
  in_iir_delay = np.concatenate([[self.prev_input[i]],in_iir[0:-1]])
  # Update this for next in0 chunk
  self.prev_input[i] = in_iir[-1]
  out_temp = in_iir*1.0001 - in_iir_delay
  # Workaround for adding y[n-1]
  out_temp[0] += self.prev_output[i]
  out_iir = np.cumsum(out_temp)
  self.prev_output[i] = out_iir[-1]

  # VCO implementation 1/(1 - z) (old style)
  out_iir[0] += self.prev_phase[i]
  in_vco = np.cumsum(out_iir)
  self.prev_phase[i] = in_vco[-1]
  real_part = np.cos(self.k_factor*in_vco)
  imag_part = np.sin(self.k_factor*in_vco)
  out_vco = real_part + 1j*imag_part
  feedback = out_vco

Also, here is the rest of the code -

def __init__(self, samp_rate, iter):
  gr.sync_block.__init__(self,
    name="costas_loop",
    in_sig=[np.complex64],
    out_sig=[np.complex64])
  self.samp_rate = samp_rate;
  self.iter = iter;
  self.call = 0

  self.prev_input = np.zeros(self.iter, dtype=np.float64)
  self.prev_output = np.zeros(self.iter, dtype=np.float64)
  self.prev_phase = np.zeros(self.iter, dtype=np.float64)

  self.k_factor = -5/samp_rate

def threshold(self, in0):
    output = 0
    if in0 <= math.cos(7*math.pi/8):
      output = -1
    elif in0 > math.cos(7*math.pi/8) and in0 <= math.cos(5*math.pi/8):
      output = -1/(2**0.5)
    elif in0 > math.cos(5*math.pi/8) and in0 <= math.cos(3*math.pi/8):
      output = 0
    elif in0 > math.cos(3*math.pi/8) and in0 <= math.cos(math.pi/8):
      output = 1/(2**0.5)
    elif in0 > math.cos(math.pi/8):
      output = 1
    return output

This code unfortunately "almost" works. I get a stabilized 8-PSK output in the beginning,

Unfortunately, after about 160 calls to the costas_loop's work(), the 8-PSK output starts to move and destabilizes. What could be going wrong?


texts:
I've been completely stuck on a portion of my assignment for a few days now. After plenty of searching around, I have been unsuccessful in discovering information that leads me to the correct solution. That said, this is for OpenCV in Python, using Numpy for matrix calculations.

The problem statement:

  Construct the derivative of Gaussian kernels,  and  by convolving the above two kernels: =∗; =∗

Additionally, here's a screenshot of the problem in full form:

What I Have Attempted:

import cv2
import numpy as np

def my_Normalize(img):
if len(img.shape) != 2:
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) #convert img to grayscale
img_float = cv2.normalize(img.astype('float'), None, 0.0, 1.0, cv2.NORM_MINMAX)
return img_float

def my_DerivativesOfGaussian(img, sigma):
#the two 3x3 Sobel kernels
#Sx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], np.float)
#Sy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], np.float)

sobelx = cv2.Sobel(img,cv2.CV_32F,1,0,ksize=3)
sobely = cv2.Sobel(img,cv2.CV_32F,0,1,ksize=3)

#the Gaussian kernel, taken from example program 9
halfSize = 3 * sigma
maskSize = 2 * halfSize + 1
gKern = np.ones((maskSize,maskSize)) / (float)( 2 * np.pi * (sigma**2))
xyRange = np.arange(-halfSize, halfSize+1)
xx, yy = np.meshgrid(xyRange, xyRange)
x2y2 = (xx**2 + yy**2)
exp_part = np.exp(-(x2y2/(2.0*(sigma**2))))
gKern = gKern * exp_part

#convolve the two kernels
Gx = cv2.filter2D(sobelx, -1, gKern)

IxN = my_Normalize(Gx)
cv2.imshow('Gx Normalized', IxN)
return Gx, Ix

img1 = cv2.imread('TestImg1.jpg')
Ix, Iy = my_DerivativesOfGaussian(img1, 1)

cv2.waitKey(0)
cv2.destroyAllWindows()

Here is a side-by-side comparison of the result I was given for a test image (left image), compared with the result I obtained from the above execution (right image):

Where I am confused:

The resulting image that was obtained is close to the target image, except blurred slightly. The blurring leads me to believe there's an error with the Gaussian kernel. I also obtain a drastically different image when I switch the order of the Gaussian and Sobel kernels in the cv2.filter2D() function. However, this is counter-intuitive for me, as the association property states that convolution should be equivalent, regardless of the order in which the convolution occurs.

Specific Questions:

Is the cv2.filter2D() function able to convolve two kernels?
Assuming the answer to ^ is false, do you know of a function that is able to perform such an operation?
Where does my logic/understanding begin to fail?  (i.e. what topics/concepts am I completely not understanding?)

Thanks so much for your time, any help is appreciated.  This is the first time I have been absolutely stumped on an assignment, and can not think of a way to solve this problem without brute-forcing the convolution formula via nested loops.

Here are the images supplied for verifying correctness of the convolution, in their native formats:

Base input test image

Normalized Ix

Normalized Iy


texts: Finding dependence between two sparse signals

texts:
I am new to signal processing and would really appreciate your help. happy to provide more context/details.

In my research, I have a system that has many sensors $X_1, \ldots, X_N$ and using this information, the system makes a prediction on a value, $A$. For simplicity, let's assume that $A$ is binary, and please keep in mind that this is a real system I am trying to improve so that means that I can tell whether $A$ was right or wrong (the system can be checked against "ground truth").

Now each sensor makes measurement sporadically (not necessarily randomly, but definitely not some obvious pattern) for reasons that are too complicated to go into right now. This means that $X_1$ for example can resemble this: $[\textrm{2, 5, 1, na, na, 1, na, na, 89, $\dots$}]$ where each index of this vector is a time step, and $\textrm{'na'}$ represents times when the sensor did not take a measurement. This means that the vector actually contains two sets of information: the value of each individual measurement, and when the measurements were done. There is a third value (but maybe that's not important for now) which is about how much effor each sensor put into collecting this information.

Now my problem is that I want to measure the dependence between $X$'s. i.e. given a specific vector of measurement from, say, $X_4$, what is the probability that $X_{10}$ was influenced by $X_4$. Specifically, I am looking for a number (could be a probability, or a correlation coefficient) that would say that $X_{10}$ was influenced by $X_4$, and I want this number for all pairs of $X$'s.

What I am trying to do is to factor in these influences into the system that is trying to make a decision. Since I know the ground truth, I can tell how good my "dependence measure" is.

I have zero background in signal analysis and I thought that you all might know something here. I think this is an interesting general problem and I would really appreciate some help.


texts: State Space conversion of Sinusoidal model

texts:
i would like to learn how to convert  sinusoidal model into state space form which has  following equation

our model consist of  sum of periodic components with additive of white noise, given by following form

$y(t) = A_1\sin(\omega_1 t) + A_2\sin(\omega_2 t+) + ... + A_p\sin(\omega_p t) + z(t)$ where  $z(t)$ is additive white noise and

$\omega_i=2*\pi*F_i $

i was reading special  document  Tracking a Sine Wave

i somehow understand step when we have  single sinusoidal components, but what about  when we have several ones? please help me


texts:
We know that for two signals $x[n]$ and $h[n]$ such that :
$$
y_{1}[n]=(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]
$$
We can deduce that :$$y_{2}[n]=x[n+2]*h[n]=y_{1}[n+2]$$ and $$y_{3}[n]=x[n]*h[n+2]=y_{1}[n+2]$$ but can we know what :
$$
y_{4}[n]=x[n+2]*h[n+2]
$$
gives in terms of $y_{1}[n]$? I am unaware if I might be asking a really stupid question or am I not seeing it? Is it $y_{1}[n+4]$?


texts: What is the significance of the linear phase in the phase response of an M-point moving average filter?

texts:
I have plotted the magnitude and the phase responses of an M-point moving average filter, the following are the plots when M = 10:

We can observe that corresponding to every lobe in the magnitude spectrum we have a line in the phase spectrum - a linear phase response is observed corresponding to every lobe. (See the below picture)

I tried to observe the magnitude and phase responses for different values of M and observed a similar pattern.
Is there any reason behind this pattern? Is there any significance for this?


texts: Why is the value of $p_i\log(p_i)$ in entropy formula considered $0$ when $p_i =0$? I get that it is the limiting value, but does that mean $p$ is never equals to zero, but only tends to?


texts:
$T(x[n]) = ax[n] + b$

$ T(\alpha_{1}x_{1}[n] + \alpha_{2}x_{2}[n]) = \alpha_{1}ax_{1}[n]+ \alpha_{2}ax_{2} [n] +b $

$ \alpha_{1}T(x_{1}[n] ) +  \alpha_{2}T(x_{2}[n]) = \alpha_{1}(ax[n]+b) + \alpha_{2}(ax_{2} [n] +b) $

therefore they are not equal so the system is non-linear?

or would it be something like this:

$ T(\alpha_{1}x_{1}[n] + \alpha_{2}x_{2}[n]) = \alpha_{1}ax_{1}[n]+\alpha_{1}b + \alpha_{2}ax_{2} [n] +  \alpha_{2}b $

$ \alpha_{1}T(x_{1}[n] ) +  \alpha_{2}T(x_{2}[n]) = \alpha_{1}(ax[n]+b) + \alpha_{2}(ax_{2} [n] +b) $

and therefore they are equal so the system is linear?


texts: How to find f given that f(some_image) = processed_image

texts:
I have a before/after image pair with pixel transformations only (for example, grayscale, or a hue shift. No blur, border, sharpen, etc.). I'm trying to find a function such that for every pixel (r, g, b) in the original image, f(r, g, b) = (r', g', b'), (r', g', b') being the pixel in the "after" image.

I was thinking about a polyfit function, but I couldn't find one that works in 3D.

I'd appreciate any help with this, or another method of achieving what I need.


texts:
Having only dealt with digital filter scarcely, the question dawned to me when I used the firls function in matlab to design an equalizer with a certain gain response.

In general, can we prescribe an arbitrary shape as the equalizer filter response and hope to obtain a FIR/IIR filter that matches this shape within a set accuracy?

Filter Order Rule of Thumb explains the relationship between filter order and the required ripple, rejection and transition bandwidth. But my question is slightly different. The equalizer I'm designing doesn't have a stringent requirement on rejection or transition bandwidth, what matters is the in-band amplitude/phase matching the prescription. For example, the prescribed shape could be predominantly a simple straight line with a slope (to equalize linear gain slope), or more practically on top of that it could have slight perturbations here and there of higher orders. I'm thinking if a polynomial of a certain order can model the prescribed in-band shape, then the FIR would have require the same order. Question is, has this been studied before?


texts: Thesis Topics for DSP using Matlab/Simulink/LabVIEW

texts:
I am searching for a good thesis topic on DSP with the use Matlab/Simulink/LabVIEW for my master, however, when I contact my department (we only have 1 DSP professor), he asks me to try to propose one but I cannot think of any.

I wonder if someone can help me with some ideas


texts: What are the possible ways to find spots with most concentrated contrast on an image?

texts:
I'm thinking of pixelize + threshold, but I need a more formalized way of finding the darkest spots on the photo, with tweakable parameters.
Something I can think of is the parts with lowest spatial frequency and least luminous pixels together, i.e. the spots with darkest pixels spread over a largest area (on average) should win but small black dots should be considered too, with some sort of adaptive weighing.
desired result; left image transformed to right image:

UPD: The required result is, in fact, a metric. Which should be a ratio of locally cumulated contrast over a known area.
Local points of negligibly small size should not count, no matter how contrast they are; bigger areas with mediocre contrast features should win [because their cumulative luminosity power is higher]. The closest thing I can think of is a 2D probability distribution function; kernelized, smoothed out, spatially-oriented.


texts: What is the time it takes to measure the SINR in LTE?

texts:
I would like to know how much time it takes the receiver to evaluate the SINR?

Scenario:
Assume that a base station (eNodeB) would like to transmit to one user. Let us say that it sends a pilot signal of L bits to that user. Then the user calculates the SINR (assuming that there are a set of N other interfering base stations) and compares it to a certain threshold and then feeds back 1 bit to the base station to tell it if the SINR is above or below the threshold. How much time it takes for this process?

In general, how the standard LTE works to evaluate the SINR at the receiver?


texts: Can deep learning be applied to video compression?

texts:
H.265 is based on incremental gains over time. The paradigms have largely remained the same. What are the odds of deep learning providing huge performance gains? Here is the main metric:

Perceptual quality remains the same at 2x improvement in bit-rate.
Assume unlimited compute power at encoder.
No need for speed.

Of course, I understand this is a really ill-posed problem because there are so many factors involved. But if it can be done, what would give the performance gain? Is it the following?

Larger CTUs instead of the 64x64 one in H.265. Perhaps even go towards object-based coding.
Better block matching algorithms that map longer dependencies.
Perceptual loss for block matching instead of MAD or MSE.
Etc...


texts: How to find PDF value from histogram plot in MATLAB?

texts:
I have the following histogram plot in MATLAB:-

I have calculated the occurrence for each container (approximate,) probability for each container (approximate,) and the width of each container from the histogram plot.
I am not sure how can I use these parameters to find the estimate value of  probability density function (PDF.)
Could someone please tell me which formula should I use to calculate the PDF using these parameters?


texts: Image processing and the Fourier Transform

texts:
I am trying to understand how and why the fourier transform is used in image processing / computer vision. Below is what I have gathered so far. Would my understanding of it be correct? If not, could somebody explain it to me in simple, plain english? Or, does anybody have anything to add to it? Last but not least, could somebody explain the "discrete fourier transform"?

The fourier transform decomposes an image into its sine and cosine components. Put simply, sine and cosine are waves starting at a minimum and maximum respectively. In the real world, we can't tell whether a wave that we observe started at a maximum or minimum point, and therefore we can't really distinguish between the two. Therefore, sine and cosine are simply referred to as sinusoids.

When applying the FT to an image, we transform it from its spatial domain into a "frequency domain", which in essence is the image represented in terms of its variation in colour and brightness over time (well, not time, but space. That is, over a number of pixels).

EDIT: Why would I use the Fourier Transform? And what are its benefits over other methods? For example, one application in literature is in shape recognition or noise elimination. In basic terms, how could one go about shape recognition using the FT?


texts: Concepts Surrounding First Two Samples in Difference Equation

texts:
I'm attempting to construct a difference equation in matlab for the function u(n) and find the auto-correlation matrix R.

I'm following the technique as is defined by my text below:

The system I'm trying to model is also defined by the text below. My a1 and a2 are different.

I have it defined in Matlab as is shown below and included results for my auto-correlation matrix. Let me clarify that I know that technically R is wrong because I'm not computing the expected value. I hadn't figured out how to do that yet and am attempting to tackle what I perceive to be the larger problem first.:

a1 = 0.1;
a2 = -0.72;

m_max = 101; % Sample Length

v = randn(1,m_max); % Create m samples of Gaussian white noise
                % with 0 mean and unit variance

u = zeros(1,m_max); % Temp matrix for input vector

for n=3:m_max
    u(n) = v(n) - a1*u(n-1) - a2*u(n-2);
end

n = [1:1:m_max]; % Reconstruct n matrix to plot

u = transpose(u);
u_H = u';

R = u*u_H; % Corrleation Matrix

figure(1);
plot(v);
axis([1 m_max -5 5]);

This runs without any problem, but the issue is that my first two values are always zero. This makes sense in context of the code because I made the zero, but I know conceptually my correlation matrix R, as is shown by image 1 above,should not be 0 for r(0). Because my for loop thought must start from 3 (because u(n-2) is not defined for n<3) I'll always have 0's in my R matrix. I suspect that I would need to know initial values but those aren't provided to me by the problem. Am I thinking about this incorrectly? Any help would be much appreciated. Thanks!


texts:
I've received my transmitted signal that I sampled it on 20480000 frequency , transmitted frequency 868MHZ , bitrate 100KHZ , I see the "bits" that Im transmitting in my plot in matlab it looks like:

so what Im now going to do is doing LPF in the cutoff frequency = bit rate and that's in order to be able to do zero crossings, so my output signal of LPF must be like "sinusoidal", I mean by that "smearing the ups and downs that we see in the photo above to look like continues sinusoidal.
Any help how can I do that in matlab? thanks
attaching down a photo that I showed what I mean by "smearing" by dark marker..

eyediagram:


texts: Relationship of input-outputs in Hartley transformation based on few elements in a vector

texts:
This question is related also to that question HERE, but with further details and another issue I need to solve.

Assuming we have a vector $X = [x_0,x_1,\ldots,x_N]$, its Hartley transformation is equal to $H(X) = \Re(FFT(X)) - \Im(FFT(X))$. where $H$ denotes the Hartley transformation, and $FFT$ is the discrete Fourier transformation. Which means that $y = H(X)$. My question : Is it possible to fix some values in the vector $X$, for example,  $X_{1:4:N}$ to be any values in order to have their correspondent values $y_{1:4:N}$ equal to  $[1,1,\ldots,1]?$. What supposed to be the valued of  $X_{1:4:N}$ to have the results of $y_{1:4:N}=[1,1,\ldots,1]?$

$NP:$ The values in both vectors $x$ and $y$ in the other index of $1:4:N$ are random data.


texts: minimum requirement of SNR for fft algorithm

texts:
i will like to ask that does fft algorithm requires any minimum snr for frequency estimation

please also refer some technical documentation in this regard


texts: Channel impulse response in freqency domain

texts:
I have a doubt about DFT matrix and channel impulse response.

Let $h$ is a channel impulse response , $F$ is a DFT matrix. IF I take a product. I got:
$$ F* h=\tilde{h}$$

$\tilde{h}$ is a frequency response of the channel, isn't? or channel impulse response in frequnce domain


texts: Frequency Modulation to arbitrary Audio signals

texts:
I understand how FM synthesis works, but I was wondering how to FM arbitrary samples or part of samples. I.E any frame of audio samples.

Further more the modulator can itself be an arbitrary frame of samples itself.

How is this done in an algorithmic manner?


texts: Difference between frequency components and harmonic components - Fourier

texts:
What is the difference between frequency components and harmonic components? The first concern the continuous domain of frequency, while the second concern the discrete domain of frequency ($f_{k}=kf_{0}$ with $f_{0}$ foundamental frequency)?

Thank you for your time.


texts: Detecting the rotation of digits in an image

texts:
I am trying to find the rotation of a given set of digits from the image. For example:

There is no additional background information, only the digits are given in the input.

I have tried using methods as described: 1 , which uses the Tesseract library to solve for the rotation, but it does not work well in this case as there are no multiple lines (of text) involved in this scenario.

For this purpose, it can be assumed that the decimal point and the "mm" are always present in the inputs. Are there any other methods or any additional features that I can exploit to solve this problem?


texts: Understanding the concept of semi-blind receivers

texts:
In the book "Signals and Images: Advances and Results in Speech, Estimation, Compression... ", an overview of using tensor analysis in signal processing is given. On page 243 it is explained that an algorithm for semi-blind receivers can be obtained by combining two tensors.

Can someone provide me with an explanation for semi-blind receivers? What is it?


texts: What are the steps for applying a 50% Overlap to an FFT?

texts: I'm currently working with audio where I apply an Hann window before running an FFT on a sample window of 4096 samples. I've been told that since I'm applying a Hann window I also need to overlap each FFT by 50%, however I can't seem to find a clear explanation of the steps involved in applying an overlap. Could anyone give me an idea of how to do so?


texts: Recovering eigenvalues from singular values

texts:
I am dealing with a problem similar to principal component analysis. Aka, I have a matrix and i want to recover the 'most efficient basis' to exaplin the matrix variability. With a square matrix these are the eigenvectors, weighted by the eigenvalues.

Originally, I was dealing with square matrices, and I used eigendecomposition to recover the eigenvectors, as explained above. Now however, I am dealing with rectangular matrices and using SVD to recover the efficient basis, ie A=USV' where the vectors of U are the recovered basis weighted by S, the singular values.

In my particular application, the sign of the eigenvalues/singular values makes a difference.

Here is my question: with eigendecomposition and square matrices, the eigenvalues will be positive/negative. With SVD, the singular values re contrained to be the absolute value of the eigenvalues, ie s_i = |lamba_i|.

Is there anyway to recover the 'sign' of the eigenvalues through SVD?

Thanks!


texts: Running the BM3D Denoising Algorithm in Python

texts:
I am looking for a strong baseline in image denoising and therefore wanted to have the BM3D algorithm in my benchmark.

These 2 python implementations:

https://github.com/ericmjonas/pybm3d
https://github.com/liuhuang31/BM3D-Denoise

have been unmaintained for a long time (and the second one is documented in what I think is chinese). I have experimented a bit with the first one without satisfactory results (see https://github.com/ericmjonas/pybm3d/issues/11 with non satisfactory fixes like clipping).

This package : https://pypi.org/project/bm3d/#description , doesn't have a documentation or source code easily findable.

Then I found this: https://docs.opencv.org/master/de/daa/group__xphoto.html#ga2fc5a9661c1338a823fb3290673a880d, in the docs of opencv but there is no indication as to how to use it and no examples in python. I saw this SO question but it's for a C++ implementation and I don't know how it would translate in Python because I am not familiar at all with opencv.


texts: segmenting region of interest (ROI) from background

texts:
I am dealing with images from biological slides and want to segment out the ROI (cells) from the background. I have tried k-means clustering, otsu thresholding, adaptive thresholding, edge detection etc. Running out of ideas now.
Which other segmenting techniques can be tried for such images? Keeping in mind that the color of cells and background can change
click here for image


texts:
I have continuous data with very high amplitude periodic noise. Data is sampled at 500 Hz. The period of the noise does change with time but within a short interval (~ a minute) it is constant. Top picture shows an interval where there is very little noise. Second picture is the noisy data. In this example the noise period is 112 samples, that is 0.22 second. Third picture shows the frequency spectra of the noisy data up to 25 Hz (sorry, freq. axis slightly shifted) and the last one is a section of auto-correlation function. Sometimes the noise is more complex (multiple periodic), but the periodicity is still there. I have tried to remove the noise by making a reference signal and some adaptive filtering in a windowing scheme. Though it more or less works, it takes too much time. I have many channels of streaming data. Not all channels are affected and between those that are, the noise doesn't correlate (directional sensors). The filter does not need to do an excellent job. Just so to allow a real signal to be detected.

Appreciate your help


texts: Fourier transform of even/odd parts of a complex signal

texts:
Why does Oppenheim state the following properties:

\begin{align}
\mathcal F\big\{x_e (t) \big\} &= \Re\big\{ X(j\omega) \big\}\\
\mathcal F\big\{x_o (t) \big\} &= j \Im\big\{ X(j\omega) \big\}
\end{align}

where $x_e (t)$ and $x_o (t)$ denote the even and odd parts of a signal $x(t)$, respectively, only for a real signal $x(t)$?
I don't understand how those two properties could break down if $x(t)$ were complex in general.


texts: When sampling theory condition fails, what should I do?

texts:
I met with a problem: if I sample audio data at 11kHz (yes, that is slow!) for 10 seconds and apply a Fourier Transform to the data, at what frequency will a 15 kHz tone appear?

So the original signal is 15kHz? which is greater than half of sampling rate 5.5kHz, right? and thus alias appears. And then, in this case, how should I find the appears one?

Thank you!


texts:
I have the following question:

a) For an N-tap FIR filter, sketch the direct form realization. How many multipliers per sample? You may choose 8.
b) For an FFT calculation of an N-tap filter sketch the realization. (N=8, FFT size is 256)

The answer to a) is given here, and I think 8 multipliers are needed per sample, assuming an 8-tap filter.
With respect to b) I am confused. I don't know how to form a realization of the FFT of a filter. Could somebody help me out with b). Thanks!


texts:
Is there a specific algorithm/formula for computing the bandwidths of filter banks across a frequency range where each bank has equal power/volume/intensity?  (The bandwidths will be different for each band, but should still all sum up to the total bandwidth).

(For example, a section split into 5 bands across 2kHz to 5kHz. At a given instant, the algorithm determines that each band needs to be 2dB, with the first band having a bandwidth of 200Hz (ranging from 2kHz to 2.2kHz), the second band having a bandwidth of 400Hz, the third band having a bandwidth of 700kHz, the fourth band having a bandwidth of 700Hz, and the fifth band having a bandwidth of 1000Hz)

My current plan of attack on this problem is something like number partitioning over a DFT vector, but I could be dead wrong.

I'm aware that the bandwidth will shift as volume/intensity changes, and I'm also aware intensity is loosely defined in this scenario, so I'm guessing some kind of envelope follower will be needed.


texts:
I intend to whiten the spectrum of a time series using a point-wise normalization in the frequency domain.

In my first test, it looks pretty good.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

The spectrum is flat after whitening.

But when I windowed the spectrum, sth strange happened.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);

k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m:-1:m-k+1) = y1(m:-1:m-k+1) .* w;
y1(n-m+3:n) = y1(m-1:-1:2);

x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

As shown in subplot 3, x1 and x2 are the same (the difference is zero; see red line), but their spectrum in subplot 4 (black and red) is different.

After IFFT y1 to x1 and FFT x1 to y2, y2 differs from y1.


texts:
I am currently a bit confused.
I am currently trying to understand the MFCC extraction process, and don't understand how the Mel filter banks only provide a scalar output from each filter bank. What does this scalar value mean? I know the filters aren't infinitely narrow, thereby only allowing one frequency so what about the frequencies around the center frequency?


texts: Do you really have to always satisfy Nyquist with all signals if your purpose is to reconstruct a signal?

texts: If not, when can you not comply with it?


texts: How the noise in time series can be reduced?

texts:
I have asked the question in the following link: Are time series data always contain noise?
I wonder if we could use variational methods to denoise the time series data just as in computer vision?
Thanks in advance...


texts: How to Improve Least Means Squares (LMS) / NLMS Filter Performance?

texts: Are there ways to increase computational performance of a normalized least squares (NLMS) filter? Multidelay block frequency-domain (MDF) filters have been proposed to do this, but they also take away from convergence speed and accuracy, because they only update the estimated impulse response once every block, not once every sample. Are there any other methods?


texts: I am a beginner to Digital Signal Processing. I took this course of DSP from Coursera in which I am stuck with this question. Since, I am new, please be descriptive. Any help will be appreciated. Thanks!


texts: Fourier Transform of a Pulse Train (Frequency Comb Simulation)

texts:
I am trying to take a Fourier transform of a femtosecond pulse train and plot the frequency spectrum.
I followed this question and am able to create a pulse train and plot both the time and frequency domain.
My question is, why do the pulses overlap in the time domain graph, since their width is on the order of femtoseconds, and the distance between them (T) is much larger?
Is it an issue with my time vector or am I interpreting the graph incorrectly?

%Using a Gaussian carrier envelope for the e-field, I show a comb of laser
%frequencies defining frep, and its representation in freq space. The
%sampling frequency is 10 times the freq (f=10*fs, fs=fr so fr=f/10).

%Define variables
fr=120e6;             %rep rate, units: Hz
fs=10*fr;             %sampling frequency, units: Hz
T=1/fs;               %period, units: s (~fs) (sampling time)
tau=100e-15;           %proportional to pulse duration, units: s
%tau=(pulse_duration)/sqrt(2*log(2));
wl=2*pi*fr;           %carrier frequency, units: Hz
ceo=0;                %offset freq., can change later

%Define electric field
L=1500;               %sets axis length
t=(0:L-1)*T;          %time, units: s

Es=zeros(length(t));
for n=0:100           %n is the number of pulses (depends on L)
    for i=1:length(t)
        Es(n+1,i)=exp(-(t(i)-n*T)^2/tau^2)*exp(-1i*wl*t(i));
    end
end
E=sum(Es,1);

%Fourier Transform to frequency
F=abs(fft(E));        %Fourier transform the vector Es
FT=F(L/2:L-1);           %Take out half due to symmetry
FTn=FT/max(FT);        %Normalized y-axis
f=(1/L)*(L/2:L-1)*fs;    %Frequency for plotting

%Plots versus time and frequency
figure(1)
plot(t,real(E))
xlim([0 max(t)/2])
xlabel('Time (s)')
ylim([-1.2 1.2])
ylabel('Amplitude')
figure(2)
plot(f,FTn,'k')
xlabel('Frequency (Hz)')
xlim([min(f) max(f)])
ylabel('Normalized Amplitude')


texts: Meaning of a formula about transient frequency bins

texts:
I'm trying to understand the following formula in this paper:

Sadly, the paper does not define most of the symbols it is using, so my guesses so far are:

$T_{tr}$ refers to the sampling interval

i.e. $\frac{1}{44100} s$

$\phi$ refers to phase

specifically: $\phi(k_{tr}, (n-1)h)$ refers to the angle of the complex vector representing frequency $k_{tr}$ for the FFT from time point $(n-1)h$ to $nh$

so basically this is saying: out of all FFT frequency bins, only include those frequencies $k_{tr}$ where above equation holds true

I tried googling around to see if there's a standard definition of "transient frequency bins" (like there is for FFT frequency bins) this is referencing, but no success. I also don't understand

Why would one compare radians $\phi$ with a unit of time $T_{tr}$?

Is this maybe a stand-in for 0 to account for sampling errors?

Why isn't $T_{tr}$ just called $T$ as it seems to be a property of the signal and unrelated to the transient bins?

(If this is not the sampling interval, then how would one obtain $T_{tr}$?)

How does this formula define $K_{tr}$?


texts: Bandpass image filtering for psychological experiment. How to choose parameters?

texts:
For a psychological experiment I would like to reproduce this image filtering in python. The paper describes it as

"we decomposed the original faces into six bands of spatial frequencies of one octave each — at 2.81, 5.62, 11.25, 22.5, 45 and 90
cycles per face"

I was thinking of using a Difference of Gaussians, but my question is, how to choose the appropriate value for both sigmas? How to relate the sigma values to the pixel width of the face?
Would I be better using a different approach?
Thanks


texts: Carrier recovery with the Costas loop : PLL pull-in range

texts:
I have created a channel model in GNU Radio (a modified version of the default GNU Radio Channel Block). Its performance under AWGN (without timing, phase or frequency offsets) aligns properly with theory (for BPSK). I'm currently evaluating the effect of carrier frequency offset on BER by means of the flowgraph below. The block "debug_bpsk_phase_recovery" is a modified version of the default GNU Radio Costas loop (provides a way of configuring the loop damping factor). The curves below were obtained by running the simulation with loop bandwidths $B_{n} = 0.1$ and 0.05 and two different data rates (250k and 500k). The frequency offset is $0.01R_{s}$.

Out of curiosity, I'm trying to find out how much frequency offset the Costas loop can handle. To my understanding, the maximum frequency offset is related to the PLL pull-in range. From theory, the PLL pull-in range is $\left( 2\pi\sqrt(2)\times \zeta\right)\times B_{n}$ (where $\zeta$ is damping factor (which I set to 1.0 for critical damping (some texts use $\sqrt(2)/2$). This evaluates to $0.88\times R_{s}$ and $0.44\times R_{s}$ for $B_{n} = 0.1$ and 0.05 respectively ($R_{s}$ is symbol rate).

However, as seen in the table below (SNR = 6dB), the loop fails to handle offsets above $0.04R_{s}$ (Bn = 0.1) and $0.03R_{s}$ (Bn = 0.05), which is much much less than theory. Any ideas of why is this happening?

As I understand (from the Mengali book), course carrier frequency (like FLL block in GNU Radio) can reduce offsets to within 10% of the symbol rate. I therefore tried to include the FLL block. This didnt help at all. As a matter of fact, the FLL introduced BER degradation even for small frequency offsets (~0.01) that the Costas loop was able handle on its own.
Why does the FLL show this behavior?

The book also suggested the use of fine carrier recovery (tracking?) (clock-aided, decision-directed) before phase recovery. Any suggestion on algorithms that I could use for this are welcome.

Thanks in Advance,
M.


texts: Calculating the phase shift between two signals based on samples

texts:
I am comparing two signals in MATLAB Simulink for finding the phase between them. To do this I am inspired by using the code found here.

I have two vectors of the same size which are a collection of samples of the two signals (sampling is more than fast enough). They are sine-signals with mostly the same frequency. To find the phase between them the code below is used

signal_one; %A vector with values of signal 1
signal_two; %A vector with values of signal 2

dot_product = dot(signal_one,signal_two);
norm_product = (norm(signal_one)*norm(signal_two);
phase_shift_in_radians = acos(dot_product/norm_product)

This code works satisfactory, and I am able to get the phaseshift from it, I just don't understand why. Can anyone show me why this gives the phase shift? A mathematical demonstration or a reference to an equation would be great.


texts: How to compare the syllable extracted from the audio

texts:
I need to recognize just one bird. I am trying to compare the syllables extracted from the sound samples. As a result of the spectrogram function, the power spectral density is taken from the sound. They are 2 dimensional values. I like to compare this test signal with a known signal which is another 2 dimensional value. Could you please suggest what would be the best function to use for comparison. Please keep in mind that the syllable of the test bird could be time shifted or they closely match the original bird. I am attaching the sample spectrogram and its syllable extracted here after filtering..

Regards


texts: A [Very] Rough Approximation to "White Noise" with Sine Waves

texts:
This might be a silly question, but I'm pretty rusty on my signals. The problem is the following:
I'm using a program which does not have a random signal generator, but I can sum up sines; and I need to have a "white noise" signal.

What would be a good way to "approximate" a white noise signal with sines?

I know that I'll need an infinite amount of sines... but understand that I can't do that, I'm limited to the program limitations. So my specific questions are the following:

1) Should the frequencies of the sine waves be below or(and?) above the Nyquist frequency?

2) Say I've chosen I want "white noise" in frequencies 10Hz-100Hz, and I can only use 10 sines. Should I equally space the frequencies inside this range? Or choose them at random?

3) Now that I have my frequencies, what magnitudes should I assign these sines if I want to have a specific "Gaussian" distribution of my "white noise" signal values?

Hope this question doesn't cringe upon anybody that likes "pure" white noise!


texts: Perform 2D convolution on a 2D unit step function

texts:
I am stuck with the following 2D convolution:
Consider the sequence $f$ defined by the following:
$f[m,n] =
  \begin{cases}
       1 & \text{0≤ m ≤ n} \\\
       0 & \text{otherwise}\\
  \end{cases}$
Determine the convolution of $f$ with itself.
Basically, I have to find the 2D convolution of $f[m,n]$ with itself. I am unsure how many cases I’ll need because this convolution is infinitely large. Also, I’ve unsure of what my equations should be. I was told that my answer should be analytic and not graphical so I need to find a formula.Any help?


texts: I'm looking to do my senior capstone project in the field of medical tomography, what is a good book to get acquainted with the field?

texts: I'm a junior currently pursuing a mathematics degree with aspirations toward medical school. I have a senior project I have my entire last year to devote time to, and I've decided that I want to tackle the challenge of medical signal processing. I'd like to write a program that reconstructs a 3d-model of the object from the cross-sectional images provided to it, very much like a CT scan. However, I have no idea where to start. You can assume that I have a good level of mathematical knowledge or the drive to acquire what is needed, so don't be afraid to throw some academic papers my way! I'd really like to dive in to my research as soon as possible, but I can't quite figure out where to start. What materials would you recommend to someone starting out, such as myself?


texts: Values of angular frequency of discrete-time sinusoidal signal

texts: Why are the values of angular frequency of discrete-time sinusoidal signals in between -pi and pi?


texts: Two different answers while doing Inverse Z-transform

texts:
Given, $$X(z) = \frac{z}{3z^2 - 4z + 1}$$
Question 1

I need to calculate inverse z-transform for ROC $|z|>1$

When I try to calculate inverse $z$-transform using partial fraction of $X(z)$ and using long division method, I seem to get two different answers.

Let us use partial fraction first. Since $X(z)$ is already in its proper form, the partial fraction evaluates to:
$$ X(z) = \frac{\frac 3 2}{z-1} - \frac{\frac 1 2 }{z-\frac 1 3}$$
On rearranging,
$$X(z) = \frac 3 2 \frac{z^{-1}}{1-z^{-1}}-\frac 1 2 \frac{z^{-1}}{1-\frac 1 3 z^{-1}}$$
Taking inverse z-transform,
$$x(n) = \frac 3 2 u(n-1) - \frac 1 2 (\frac 1 3)^{n-1}u(n-1) \tag{1}$$

Now on using long division method, the first few terms of the division are:
$$\frac 1 3 z^{-1} + \frac 4 9 z^{-2} + ...$$
This suggests that the signal in time domain is:
$$x(n) = \{0, \frac 1 3, \frac 4 9 , ...\}$$

But when I substitute $ n = 0,1,2,...$ in equation (1), I get the sequence:
$$x(n) = \{0,1,\frac 4 3,...\}$$

Why am I getting two different sequences for the same function $X(z)$?

Question 2

The notes my teacher provided me had this exercise done by finding partial fraction for $\frac{X(z)}{z}$ rather than $X(z)$ itself. As the given function is already in the proper form (order of denominator is greater than the order of numerator), is it incorrect to calculate partial fraction of $X(z)$ itself?

This has eaten almost three hours of my time. Any help would be appreciated. Thanks.

EDIT

For when the ROC of $X(z)$ is $|z| < \frac 1 3$, $x(n)$ should be anti-causal. While shifting $x(n)$ due to presence of $z^{-1}$ in the numerator, which one among the following two is correct?

$$x(n) = - \frac 3 2 u(-n) + \frac 1 2 (\frac 1 3)^{n-1} u(-n) \tag{2}$$
or
$$ x(n) = - \frac 3 2 u(-n-2) + \frac 1 2 (\frac 1 3)^{n-1} u(-n-2) \tag{3}$$


texts: Better way to turn time series spikes/high frequencies into discrete signals?

texts:
I am just wondering if there is a better way to transform time series data spikes/high frequencies into discrete signals?
Currently, I have time series data that I am applying a median filter to loop through and filter large spikes.  I am just using the standard deviation to find the lower and upper regions of my stationary signal.  Once everything is removed but the spikes, I take the sign of the value to get a discrete signal.

Currently it works but I’m wondering if there is a better method that can be utilized...


texts: time frequency localization using wavelet transform

texts:
I am currently doing analysis on Photoplethysmograph (PPG) data and I want to know the frequency (heart rate) at every time point so a windowed FFT might not be the best option. I am looking at wavelet to generate frequency and time information. I have been working with matlab example code however I have trouble determining what is the best wavelet to use for this application. I do now know much about wavelet. Does the wavelet used being dictated by the shape of my signal. In my case PPG signal which looks like this:

Or some other consideration is needed when choosing wavelet? How critical is the choice of wavelet?

Thanks,
Kelvin


texts:
I am currently using overlap-and-add method for a noise cancellation project. Everything works well so far except that when it comes to the real time, I need to overlap the input signals for FFT input. For example, I am using 2048 points FFT with 1024 overlapping samples from previous frame(50% overlapping rate). I am illustrating here with numbers. Each number represents half frame(1024 samples). So the real time input signal goes like 1-2-3-4-5-6-7-8-9-10, but I want to buffer the signal into 1-2-2-3-3-4-4-5-5-6-6. Is is possible to do that?

I have tried to delay the input data audio stream by half data frame samples(1024). So at the same time, I have to data streams like below:

 1-|2|-|3|-|4|-|5|-|6|-|7|-|8|-|9|-|10|           --input data stream
   |1|-|2|-|3|-|4|-|5|-|6|-|7|-|8|-|9|-|10|       --delayed data stream
    ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑

at any given clock cycle, I would want to read like the arrow indicated to feed the input for FFT. But since FFT and the input are operating at the same clock frequency, I can't really get the data for FFT input.

thanks
yong


texts: real time signal analysis of output from an audio device

texts:
I'm writing a test program that can detect pops/cracks/noise (unknown frequencies and duration) in the audio output of an audio device. The purpose of this program is to reduce human error and make the test more systematic.

The test consists of playing a pure sine tone (say, 1kHz) by the test device, and recording the analog output with a recording device. The minimum test duration is 8 hours and any artifacts present in the output will be logged and saved for future access.

Currently, my approach is as follows:

Sample Rate(playback device): 48000Hz

Sample Rate(recording device): 48000Hz

Windowing function : Hanning window

FFT size = total samples

Record 10 seconds of audio samples.
Apply window, run FFT on the recorded samples.
Record the next 10 seconds of samples.
Apply window, run FFT on new samples.
Compare both FFTs. If the frequency bins have similar power (to a set episilon), then the audio signal is clean. Otherwise, log the timestamp of this occurence.
Repeat steps 1-5 until 8 hours have surpassed.

A 1kHz sine tone, along with a 0.01s 3.5kHz sine tone that repeats randomly, is used as a test signal for the program. Unfortunately, my results have a high frequency resolution and it seems difficult to detect what is considered noise.

I am not well-versed in DSP and was wondering how I can improve the accuracy of my results?


texts: Contourlet transform for image compression

texts:
Before today I had no idea there is any transform an such; I was only aware with DCT and wavelet transforms but today one of my teacher who is doing research on image processing project had told me that contourlet transform is way better than wavelet transform.

It is better because wavelet transform extract frequencies in 3 directions (diagonal, horizontal, vertical) while contourlet transform is capable of extracting information in any direction.

There are various books on wavelet transform but I couldn't be able to find a book that mainly focus on contourlet transform, is there any book you would recommend for contourlet transform?

P.S: I'm not talking about research papers.


texts: Can ICA be applied, when the number of mixture signal is less than number of source signal?

texts:
I am referring to the the following paper : Non-contact, automated cardiac pulse measurements using video imaging and blind source separation

In the above article, the authors are able to extract cardiac pulse signal out from RGB components. I try to visualize the process as follow.

R' = R + cardiac pulse
G' = G + cardiac pulse
B' = B + cardiac pulse

R', G' and B' are the colour components observed by camera. R, G, B are the colour components for a person, by assuming that he doesn't have any cardiac pulse.

It seems that we will be having 4 sources (R, G, B, Cardiac pulse). We are now trying to obtain 1 of the 4 sources (Cardiac pulse) from 3 mixture signals (R', G', B'), by using ICA.

Does it make sense? Am I missing some techniques? Or, am I making a wrong assumption on the process?


texts: How to Find the Frequency Response of a Communication Channel from Input and Output Symbols in MATLAB?

texts:
I want to find the channel frequency response of a digital communications system. I have the functions of the input symbol (a triangle) and the output symbol - a distorted triangle: $$\dfrac{1}{1 + \left( \frac{2t}{T} \right)^2 }$$
I want to plot the frequency response of the channel as if it was a digital filter.

How can I do this in MATLAB?

I tried using symbolic variables, computing Fourier transforms of the input and the output signal using fourier command, and then finding the numerical values using subs and double commands, with no success.


texts: Are Chebyshev Type 1 filters very unstable at very high orders?

texts:
I am currently experimenting with a digital implementation of a Chebyshev Type 1 filter of mine, for audio purposes. I have written the code myself so I was wondering if the behavior I've noticed is due to bad code or something else. My implementation works perfectly fine up to an order of 30. I know that this is ridiculous but that's not the point, it's just for educational purposes.
Anyway, when using a filter of order 30 or higher part of the frequency response that should be flat has resonant peaks. For a lowpass filter, it meant boosted low frequencies for a highpass filter it resulted in boosted low-mids. The highpass filter did not behave this way when the cutoff was above the resonant peaks, by the way, and changing the ripple didn't change the trend.
Is this normal and a result of the filter getting more and more unstable? Or does this indicate that there is some part of the code not working properly?
I am using double precision where I can, but the framework I am working in wants the sample values with 32-bit precision so I can only do the filter coefficients and filter memory with 64-bit precision and have to reduce precision when writing to the audio buffer. Maybe that's an issue?

Edit: To add to my OP, I am using cascaded biquads to realize the filter.


texts:
Can someone explain waveshaping to me?

I only know it shapes waveshapes by passing through some functions, but I don't yet understand e.g. what the plots (such as the following in Melda Production MWaveShaper VST-plugin) mean.

Both axes have dB values so I read the axes to be input and output. However, how is the distortion "visualized" in the graph? How does one know how the graph alters audio?


texts:
It seems like (controlled) noise can sometimes be good: dithering in an image, or a softer sound in audio, or white noise alone. Here's a paper on noise improving detectability. My intuition also says that noise in a financial trading algorithm could act as a sort of hedge.

I have few examples — only a vague philosophical instinct. Is there any general-purpose literature or theory on why and how noise can be good? Or other/better examples?


texts: how can i perform a registration to 2D images using ''rigid " transformation in matlab?

texts:
i need to perform registration to 2D image using"rigid transformation" in  matlab .i need this for measuring the similarity between my images .any one know how this can be done using matlab?

an example on my images :

regards


texts: What is the definition of half duplex?Is this system,$A$ give $x$ to $B$,then $B$ give $x$ to $C$,called half duplex?

texts:
I saw the definition of half duplex in wiki,and it said

half-duplex (HDX) system provides communication in both directions, but only one direction at a time (not simultaneously). Typically, once a party begins receiving a signal, it must wait for the transmitter to stop transmitting, before replying.

However,i saw a paper: http://jultika.oulu.fi/files/nbnfi-fe2019040110636.pdf ,its system model is a transmitter transmit the $x_1$ and $x_2$ information to the user1(U1),then U1 transmit $x_2$ to the user2(U2),and it said this is half-duplex.

So in fact,the definition of half duplex is that device can transmit and receive,but no need to transmit and receive at the same time?

I mean,$A$ give $x$ to $B$,then $B$ give $x$ to $C$,this is also called half-duplex?I thought this is called simplex duplex before,because i thought the half-duplex is $A$ give $x$ to $B$,then $B$ give $x$ back to $A$.

So i want to ask that $A$ give $x$ to $B$,then $B$ give $x$ to $C$,this is also half-duplex?


texts:
I have a source (that I can't share) that states that bilateral filtering isn't shift-invariant (that is, it depends on the input signal), hence we cannot use convolution and the Fourier convolution theorem (that is, convolution is the spatial domain corresponds to multiplication in the spectral one).

Why isn't bilateral filtering shift-invariant? Why exactly can't we apply convolution and the convolution theorem to the bilateral filtering of an image? Why wouldn't a Gaussian filter also be signal dependent?

I understood that bilateral filtering combines domain and range filtering, that is, when updating a pixel of an image, it weights the contribution of other pixels in the image using both their Euclidean distance (to the pixel being updated) and their intensity (in case of grayscale) value. However, I don't get why we couldn't also use the convolution theorem to filter an image using BF.


texts: How to identify the beginning and the end of a signal using MATLAB

texts:
I have multiple signals as shown in this image. I'd like automatically detect time for the earliest detection of the signal (approx 0.7e4 in S1-B3-B (orange)) and last detected occurrence of the signal (approx 4.2e4 in S2-B3-B (purple)).

Is there a smart way to that? I'm quite new to signal processing.


texts: The effect of the mixed signals probability distribution on the ICA performance

texts:
Does the mixed signals distribution affect the ICA performance? I mean the ability of ICA to get the sources.

Assume that the mixed signals follows a Gaussian distribution and of course the sources do not follow Gaussian distributions. Then, this implies that we have large independent components or sources, so will the ICA be able to get all of that components?
On the other hand, if the mixed signals follows a Poisson distribution, should that mean something regarding the ICA method? I mean should we expect that the ICA will not be able to extract the sources of that signals in an efficient way?


texts: How to plot even and odd parts of log-Gabor filter?

texts:
I am reading a paper inwhich there is a plot of log-Gabor quadrature pair which are even and odd.

I'm wondering how these two are plotted in time-domain, since log-gabor is real valued in frequency-domain and I cannot use inverse Fourier of real/imaginary part to result in even/odd parts in time-domain!

Does anyone have an idea?


texts: Is a radix-4 implementation faster than a equivalently well coded radix-2 FFT?  And if so, why would it be faster?


texts: Half of fundamental frequency in Fourier transform

texts:
I understand why there are peaks in integer multiples of the fundamental in the Fourier transform. But why is there often a peak in half the fundamental?

I am testing for the guitar specifically. I am doing an STFT on an guitar audio monophonic signal and then get the 4 strongest frequencies in each note. I always get the fundamental in those 4 and of course some multiples of it but also half of it.


texts: Justification for Squared $ {L}_{2} $ Data and Smoothness Term as an Error Bound

texts:
Often in variational methods (and not only) we have an energy that is of the form:
$$E(u) = \frac{1}{2}\|f-u\|^2_2 + \frac{\alpha}{2}\|\psi(u)\|^2_2,$$

where the first term is referred to as the data term, and the second as the smoothness term. I understand that the squared $L_2$ norm is especially appealing since it results in a simple expression when we are looking for the minimum. I am looking for the formal motivation for designing such an energy. More specifically, I would like to see it as resulting from some bound upon the actual error. Let $f$ be a degraded/noisy image of the original image $g$, and $u$ be the image that we are looking for. Is it possible for the above to be interpreted as a bound of $\|g-f\|^2_2$? Let us consider an even simpler problem (where $h$ is a sufficiently smoothed version of $f$):
$$E(u) = \frac{1}{2}\|f-u\|^2_2 + \frac{\alpha}{2}\|u-h\|^2_2$$
If I start from $\|g-f\|_2$ I can get the following bound:
$$\|u-g\| = \|u-f+f-g\| \leq \|u-f\| + \|f-g\|$$
$$\|u-g\| = \|u-h+h-g\| \leq \|u-h\| + \|h-g\|$$
Combining the two inequalities with some weight $\lambda \in [0,1]$ I get:
$$\|u-g\|\leq (1-\lambda)\|u-f\| + \lambda\|u-h\| + (1-\lambda)\|f-g\| + \lambda\|h-g\|$$
Since $f,g,h$ are fixed, the minimization is only over the first two terms:
$$\min_u (1-\lambda)\|u-f\| + \lambda\|u-h\|$$
And this can be rewritten through $\alpha$ by setting $\alpha = \frac{\lambda}{1-\lambda}$. My issue is that these are still not squared norms. If I try to square both sides, this results in a different minimization energy. The other way to get the squares involves writing out something of the form:
$$\|u-g\|^2 = \|(u-f) + (f-g)\|^2 = \|u-f\|^2 + \|f-g\|^2 - 2(u-f)\cdot(g-f) \\ \leq \|u-f\|^2 + \|f-g\|^2$$
The last inequality holds only when $(u-f)\cdot(g-f) \geq 0$ though. Said otherwise, the angle between $(u-f)$ and $(g-f)$ needs to be smaller than 90 degrees which is not something that can be guaranteed I believe.
A looser bound is:
$$\|u-g\|^2 \leq \|u-f\|^2 + \|f-g\|^2 + 2\|u-f\|\|f-g\|,$$
but it also involves non-squared terms.

Here's something else I noticed about the non-squared energy:

$$\alpha\|u-f\| + \beta\|u-h\|, \alpha, \beta > 0$$

will always have as solutions $u=f$ or $u=h$ for $\alpha<\beta$ respectively $\alpha>\beta$. When $\alpha=\beta$ any solution of the form $(1-\lambda)f + \lambda h$ is admissible. It is thus clear that the above is not very interesting as a minimization energy.

Cross-posted to: https://math.stackexchange.com/q/3692290/463794
I expect to get a different perspective there.


texts: Convolution with an impulse response does not give the same result as the frequency response

texts:
Context

I am trying to design a linear-phase FIR filter with the following frequency and phase responses:

Designing the filter

Given its characteristics (peak filter at 1kHz with 9dB gain and Q=7), I got the following IIR biquad coefficients:

$$b_0 = 1.0889938082113808$$
$$b_1 = -1.8806408829603927$$
$$b_2 = 0.80787505079350063$$
$$a_0 = 1$$
$$a_1 = -1.8806408829603927$$
$$a_2 = 0.89686885900488144]$$

I then computed the impulse response of that biquad filter and truncated it to 1024 taps. I got the following filter:

To make the phase of my FIR filter linear, I applied an FFT on its impulse response, I put the imaginary part to 0 and came back to time domain by doing an inverse FFT.

The final taps I got can be found there.

The problem

Now, when I plot the frequency/phase responses of this filter using Numpy/Scipy or Matlab, I get exactly the same plot as the one of my expected filter (see figure 1). But, when I apply this filter on a pink noise signal by convolution and measure the transfer function, I get this instead:

I am pretty lost here, I would assume that convolving my noise signal with my filter and measuring the transfer function would give me the same curve as the one I computed theoritically but it seems not to be the case (that is, in my mind, the definition of what convolution does).

Some additional information:

I have verified that the problem does not come from my implementation of the convolution operation, by testing with many different implementations; if you still believe the problem come from there, you can try yourself with the taps I provided.
I trust the software I used to measure the transfer function depicted in the last plot, and anyway I could hear that the transfer function is not the one I expect.
To plot the transfer functions from Python, I used the following code:

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

def plot_transfer_function(b, a=1):
    """Plot the frequency and phase response of a digital filter."""
    w, h = signal.freqz(b, a)

    fig = plt.figure()
    plt.title('Digital filter frequency response')

    # Frequency response
    ax1 = fig.add_subplot(211)
    plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), 20 * np.log10(abs(h)), 'b')
    # Make the plot pretty
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Frequency [Hz]')
    plt.axis([20, 20000, -12, 12])
    plt.grid()  # show the grid

    # Phase response
    ax2 = fig.add_subplot(212)

    angles = np.unwrap(np.angle(h))
    plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), np.degrees(angles), '#aa0000')
    # Make the plot pretty
    plt.ylabel('Angle (degrees)', color='#aa0000')
    plt.axis([20, 20000, -180, 180])
    plt.grid()

    fig.axes[0].set_axis_off()  # disable main figure axis

    plt.show()

plot_transfer_function(TAPS)  # "TAPS" is an array containing the 1024 taps

So, what is wrong with my reasoning here? Why do I not measure the same frequency/phase response when calculating it and when measuring it after a convolution?

Additional question: for my FIR filter to have linear phase, I thought it should be symmetric (or anti-symmetric). Well, when I plot my 1024 taps, they appear to be pretty symmetric, excepted the first which does not have its equivalent at the end of the impulse response. Why is it so?


texts: Can any one tell me how to design wavelets from splines using matlab. whether we can make wavelets from higher order splines or only with B splines


texts: a neural network approach for FIR filter

texts: I am trying to write a code for a neural network to do the digital filtering on some signals. Is there any neural network model for digital filtering?


texts: How to understand Kalman gain intuitively?

texts:
The Kalman filter algorithm works as follows

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

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

  Predict

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

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

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

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

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

Thanks and regards!


texts: Developing a classifier for a periodic music

texts:
I am trying to develop a classifier for Indian Taals, for which I would be focusing only on the Tabla.
Now once a extract a period of the Tala(a repetitive pattern), I want to develop some feature vectors for a Tala, which can be used to differentiate two Talas.

I am not sure what they should be. As of now, what I am planning to do is to find the beats of the Tala by beat tracking, or the onsets, and then take the median of the spectrogram values at the beat times to be my feature vector for a particular tala.  Also, there is the concept of Bayan/Dayan in Tabla which have different pitches. I am not very sure how to use this information into my feature vector.

I would also want to include the time(temporal characteristics) into my feature.

Any help would be appreciated.


texts: What could be the discrepancies between theory and practice when implementing a proportional controller?

texts:
I implemented the following proportional controller circuit in practice.

The transfer function of the plant I uses is
$\frac{20}{(s+2)(s+3)}$.
Before doing the experiment, I plotted the root locus plot of the system.
Then after doing the experiment from the plots, the damping ratio was found using the 2nd order equations. During the experiments the gain of the circuit was changed such that different responses where obtained such as an under damped response.  These were found by the finding the value of the peak overshoot from the plots.  After finding the damping ratio, the new poles of the system were found and they were compared to the root locus obtained before.
They did not match exactly.  I know that there could be practical errors such as mismatches in the resistor values in the ciruit, error when finding the overshoot amplitude from the graph etc,   but could there also be some theoretical errors when using this system.

Since I am using the general equation of the second order system, could this have an effect, I know that the plant is second order but when considering the whole closed loop system, does this still apply?


texts: Estimating the input to a system from a system state using EKF

texts:
[ Cross-posted from: https://math.stackexchange.com/questions/164169/estimating-the-input-to-a-system-from-a-system-state ]

I have a system for which I have obtained a non-linear time-varying state-space representation. For this system I am able to measure one of the states. I would like to estimate the input from this.

In my parameter estimation notes I have found a scribbled aside "[the extended Kalman Filter] can be used to estimate the system input if the output is known". So far I've not found any references or sections in my notes that describe the use of the EKF in this way. Can anyone please suggest a reference (or provide) an explanation of how to perform this task?

The system is:
$$x(t)=
\begin{bmatrix}
x_1(t)\\
x_2(t)
\end{bmatrix}; \quad \dot{x}(t)=
\begin{bmatrix}
x_2(t) \\
\alpha(t) x^2_2(t) + \beta(t) + \gamma(t) u_1(t)
\end{bmatrix}; \quad y=x_2(t)
$$
Where: $\alpha(t)$ and $\gamma(t)$ are constant parameters; $\beta(t)$ is a time-varying parameter; $y=x_2(t)$ is a known (measured) time-varying state; and $u(t)$ is the unknown system input. I have made a simplifying assumption that $\beta(t)$ is constant initially. After proving against this simplified system I want to move to a system where:
$$x(t)=
\begin{bmatrix}
x_1(t)\\
x_2(t)
\end{bmatrix}; \quad u(t)=
\begin{bmatrix}
u_1(t)\\
u_2(t)
\end{bmatrix};
\quad y=x_2(t)\\
\quad \dot{x}(t)=
\begin{bmatrix}
x_2(t) \\
\alpha(t) x^2_2(t) + \beta(t) + \gamma(t) u_1(t)
\end{bmatrix}; \quad \beta(t) = f(u_2(t))
$$
Where: $u_1(t)$ is the unknown system input ($u(t)$ in part 1); $u_2(t)$ and $f(u_2(t))$ are a known input and a known function.  If someone can point me at a reference for the first part (constant $\beta(t)$) or give a brief explanation of how to approach this then I hope to be able to work to a solution for the second part.


texts: Fourier series expansion of $\exp[-j2kA\sin(\omega t)]$

texts:
I've been going through this paper where we exploit the well known Fourier expansion of the model signal as shown in the image below.

I've never come across this well known fourier expansion before. Can someone point me towards right direction to fill in the missing steps in above image and also the use of Bessel functions got me confused a bit. Any references would be of great help. Thanks.


texts:
So I was thinking about the problem of shadows in images. From an image processing point of view, would you classify cast shadows in a scene as an additive phenomenon or as a multiplicative phenomenon on top of the shading?

From what I can understand, for the case of non-lambertian image formation equation, the specularity is formulated as additive because it is dependent on the view direction. But for the case of shadows, is there any similarity to that? Or should shadows be formulated as just as a scaling term over the shading?

Any views or pointing towards proper papers or materials would be highly appreciated. Thanks!


texts: How to apply interpolation algorithm for FFT generation?

texts:
How to decide interpolation  algorithm before applying FFT?
I need to find FFT for multi-sine signal .
Before finding the FFT, need to apply interpolation algorithm.
How to decide whether I need to apply interpolation algorithm ? Which algorithm do I apply ?
Any reference for my queries?
Samples are 50000.frequency resolution is 1Khz. If I pad zeros samples in time domain, number of samples  will be 65536.


texts: How to plot the poles and zeros of this Bandpass Filter in MATLAB?

texts:
So I'm trying to design a band pass filter in MATLAB (with a center frequency of 10kHz and a sampling frequency of 44kHz).

I have calculated the transfer function but I'm not sure how to enter this transfer function into MATLAB to generate a plot of the frequency and phase responses of the filter, and the poles and zeros.

It could also be that the transfer function is calculated wrong? I haven't seen the denominator have so many terms before so I'm not sure..

This is my first time posting on SE, so any help is appreciated.

PS: I kinda always thought that the poles of a band pass filter would be pi/2 always, so I guess I was wrong?


texts: How to prevent (or remove) aliasing when a cross-correlation between spectra is performed

texts:
This is my problem: I'm implementing a tool which cross-correlates a spectrum (fluxes as function of wavelength, obtained directly from an instrument) with a template (fluxes as function of wavelength, obtained by stacking observations). In order to perform such operation I have to dilate the template by multiplying the wavelength by a factor $(1+z)$ where $z$ is a floating value whose range is manually set (typically could be something like $[2.0, 2.5]$, always less than $7$). After the dilation, I have to resample this modified template so that it shares with the spectrum the same sampling (an upsampling).
This tool is actually working with $z<4$, for values of $z$ larger than this threshold I found there is a sort of aliasing effect.
Be honest, I'm not sure that it is aliasing, but looking the cross-correlation function, it shows the presence of something that looks like an aliasing effect: in the figure there are some results, in each plot, upper panel is the spectrum, bottom panel the cross-correlation function (the blue line is centered on the peak of this function, while the red line is the "true" value; but these lines are not important)

Now, these results are confusing me since I don't understand how to manage these oscillations in a spectral domain. How can I prevent this mechanism? Or how can I remove these "frequencies"?


texts: How can we combine 2D image with 3D image for the same sense?

texts:
We have 3D point cloud data captured by a laser scanner and need to register it with 2D images to get further information. In our case, the 2D image and 3D points are not overlapping and therefore we need to preprocess them so as to extrinsically calibrate the 3D scanner with the 2D image sensor (camera).

For the same modality, such as 2D-2D, this is done with good results using standard techniques, but for different data space, we still face challenges in this task.

These kind of images can be used in auto-driving cars and other applications.


texts: Envelop detection using matlab hilbert function

texts:
I am trying to detect the envelope of a high-frequency signal (similar to below)  using Matlab Hilbert function

This is my output

It generally follows the envelope, but there is a high-frequency ripple in the envelope, see zoomed version below

Here is my code:
t =out.tout';
x =out.Itxmod';
y = hilbert(x);
env = abs(y);
figure;
plot_param = {'Color', [0.6 0.1 0.2],'Linewidth',2};
plot(t,x)
hold on
plot(t,[-1;1]*env,plot_param{:})
hold off
xlim([0.001 0.002])
title('Hilbert Envelope')

Question: How can I get a smooth envelope the signal without this high-frequency ripple?


texts: Is there a FFT algorithm which is both in-place and constant geometry?

texts:
I work in crypto where the NTT (Number Theoretical Transform) which is just a FFT in finite fields is used. I wanted to know if there exists an implementation of the FFT which is both in-place and constant geometry? I will explain the terms below.

An in-place algorithm will not need more space than the space used to store the array. In FFT algorithms, this means the two indices accessed in the innermost loop are also the two indices where the output is written.

A constant geometry algorithm will always use the same pair of indices, at all different stages of the algorithm. For an intuitive explanation, take a look at this.

As much as I have seen, there seems to not be an algorithm which is both. I would like to know if one exists. Or if there is a reason if it doesn't exist.


texts:
Edit: Please treat code here as python-esque pseudo-code; it would syntactically fall closest to MATLAB or Python+numpy+scipy users.
If I have a signal in frames, and I want to put it back together - I just add each frame back into into a growing vector at a distance of hop from the end. More specifically, I use hop * frame_index.
Now I understand that this is naive; I believe ideally, I would satisfy a cola condition based on the type(window), len(window), and len(hop).
If I have frames, and just want to put them back together (and in particular I don't need the spectral benefits of other windows - I don't plan to go spectral) - Am I right to understand that a rectangular window always satisfies cola, so, that's all I should need here.
BUT, for a rectangular window, I would expect some gain in power for the n_overlap samples in each window - because the redundant samples get doubled. I thought this is why you needed a hann or hamm etc.
So my big question here is:
Is this understanding flawed? and if not, then how would I go about performing 2 or 3 time varying operations on a signal that each work better for different window and hop sizes.
Side-Quest: When I look into reconstructing from frames, I always come back to putting an STFT back together, or OLA convolution. Is there some convolution approach that just puts the signal back together perfectly from its frames - It makes sense that something like this could exist, but I don't fully understand how.


texts: To recover the original signal when compressive sensing is used, we get an error, to avoid the error, we need to modify the l1eq_pd.m file. How to go about it? Where can I get this file in Matlab?


texts: what is probability density function (PDF) of periodic signal

texts:
Anybody knows what is the best pdf for periodic signal ?

I have signal with $480$ samples $\left(16000*0.03=480\right)$ but I could not find the exact pdf for that. sometimes it appear like 2 mixed Gaussian.


texts: Shift property of the Short-term Fourier transform (STFT)

texts: Properties of STFT (Short-term Fourier transform) say that it preserves time shift up to modulation. Does it mean it is sensitive to time-shift in signal?


texts: Frequency estimation on chaotic waveform

texts:
The target signal is generated by the following analytically solvable hybrid chaotic oscillator

\begin{equation}\label{eq:1}
\ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2)(u-s)=0
\end{equation}

where

$\omega$ determines the symbol period as $T = \frac{2\pi}{\omega}$.
$\beta \in (0,T^{-1}\ln 2)$ affects the fluctuating amplitude.
$s \in{\pm 1}$ is the random bipolar symbol sequence driving the system.

Its typical continuous waveform $u(t)$ for $f=100 \, \mathrm{Hz}$ is as follows:

I am interested in estimating the basis frequency $f$ of the signal. However, the amplitude spectrum of $u(t)$ has no obvious feature that can be used to estimate the frequency. Hence, I calculate

$$z(t) = |u(t)| - 1$$

The amplitude spectrum of $z(t)$ is as follows

Obviously, we can see a spike at $f = 100 \, \mathrm{Hz}$. But it is hard to calculate the accurate expression of $z(f)$ since $u(t)$ is a kind of non-stationary signal.

The STFT time-frequency spectrum is shown below

As shown in the STFT results, in some time intervals, no obvious frequency component can be found. However, in the part of approximately periodic waveform, the frequency component at $f = 100 \, \mathrm{Hz}$ is strong.

Is there any theory that can support my method of frequency estimation from the FFT of the whole $z(t)$?

Here is the MATLAB code for generating the chaotic waveform from an analytical expression (not provided above) and my current experiments on frequency estimation.

clear;close all;

finit = 100;%set frequency
% finit =2000;
% finit =2;
om=2*pi*finit;
be=om/(2*pi)*log(2)-0.01*rand(1,1);
T = 2*pi/om; %Period
N=100;% number of symbol
tall = N*T;%total time
h=T*0.01;%time step
t=0:h:tall-h;
ut=zeros(1,length(t));
u_in = 2*rand(1)-1;
% u_in =0.99;
for i=1:N
    if u_in>=0
        s_in = 1;
    else
        s_in = -1;
    end
    t_f = (i-1)*T:h:i*T-h;  %Time of each frame
    n_f = (i-1)*length(t_f)+1:i*length(t_f);
    ut(n_f)=s_in+(u_in-s_in)*exp(be*(t_f+(1-i)*T)).*...
           (cos(om*t_f)-be/om*sin(om*t_f));
    u_in=s_in+(u_in-s_in)*exp(be*T)*cos(om*i*T);
end

figure
plot(t, ut)
xlabel('Time t/s');
title('u(t)')

y = abs(ut) - 1; % Observed signal

%% FFT spectrum
L = length(ut);
Ts_v = t./(0:L-1);
Fs = ceil(1/Ts_v(2));
deltaf= 0.1;
NFFT = floor(Fs/deltaf);
fs = Fs;
f = Fs/2*linspace(0,1,floor(NFFT/2+1));
Y = fft(y,NFFT)/L;
Yabs = 20*log10(abs(Y(1:floor(NFFT/2+1))));
figure;
plot(f,Yabs);
xlim([0,max(f)/4]);
ylim([-100,0]);
xlabel('Frequency f/Hz');
title('FFT Spectrum');
%%  Wavelet Transfrom Time-Frequency spectrum
wavename='cmor3-3';
totalscal=256*2;
Fc=centfrq(wavename);
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs);
coefs=cwt(y ,scals,wavename);

figure;
subplot(211)
plot(t,y,'linewidth',1)
xlim([0 1])
ylim([-1.5 1.5])
subplot(212)
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
ylim([0 200]);
% colorbar;
xlabel('Time t/s');
ylabel('Frequency f/Hz');
title('Wavelet Time-Frequency Spectrum');
%% STFT Time-Frequency spectrum
window_len = 10*2*pi/om/h;
noverlap = 50;
nfft = NFFT;
[~,F,T,P]=spectrogram(y,hamming(window_len),noverlap,nfft,fs);

figure;
subplot(211)
plot(t,y,'linewidth',1)
xlabel('Time t/s');
title('|u(t)|-1')
xlim([0 1])
ylim([-1.5 1.5])
subplot(212)
surf(T,F,abs(P),'edgecolor','none')
axis tight;
view(0,90)
ylim([0 200]);
xlabel('Time t/s');
ylabel('Frequency f/Hz');
title('STFT Time-Frequency Spectrum');
```


texts: The higher the sidelobes of a signal frequency spectrum, the greater the "cross-talk" between transmission channels

texts:
Why is that 'The higher the sidelobes, the greater the "cross-talk" between channels' when transmitting signals via hardware electronics?

where,

Cross-talk:  a signal transmitted on one circuit or channel
of a transmission system creates an undesired effect and/or interference in another
circuit or channel.
Sidelobes: The sidelobes (spectral leakage) introduced to signal frequency spectrum after
application of window functions to the transmission signal.


texts: Subtracting one FFT from another for noise reduction

texts:
I am using an FFT/iFFT library in c#, and working with arrays of audio & complex data.

What I am hoping to do is take a small amount of 'background' audio signal, and use it as a fingerprint to subtract that background signal content from another 'real' signal.

The fingerprint (F) piece of signal is an FFT transform of audio 1 window size in samples. The 'real' (R 1..n) signal is a series of FFT transforms of many windows.

this is essentially how I understand things after trying to figure it out with a text editor, some basic maths and google ...  If anyone can point me at a step by step that a software dev with some dsp interest and basic math (intermediate maybe) can understand I would appreciate it :)

I initially, expecting problems, just subtracted the real part of each bin in F from the real part of each correpsonding bin of each window of R. The iFFT of the results produced a noisy signal around the window boundaries.

Then I (think) understood I need to be subtracting the magnitudes represented, not just the real parts, - and that this was in a polar coordinate system when doing the math. So I thought, just subtract the magnitudes - the phases shouldn't change?  Then came to realise (maybe?) the phases were not aligned between the fingerprint F and each real R signal window. So now I wonder can you align the phases of the bins in the fingerprint with those of each of the real signal windows, or vice versa, before subtracting the magnitudes? and then iFFT wont be noisy?

This is a wordy question, I understand. I think what I am really asking is,  can someone help me straighten out my understanding of what's happening/required maths wise? I think I have a good mental model of whats going on, just struggling to understand how to manipulate it?

I understand each bin represents an average of a band of frequencies - hence also phases.. I think.. though this confuses me then as to how iFFT can reproduce the signal, though I think I have that pinned down to sample rates.

I think I understand the concept of 'windowing/convolving' the signal, but that seems like a technique to hide the problem (unless you could convolve in a way that essentially also aligns phase?).

After here I get a bit lost....


texts: DFT without knowing the exact value of sampling rate?

texts:
I have a sequence of (real) numbers that represent the magnitude of a certain natural event. I know that the samples are equispaced in time, but not the exact value of the spacings. So does that mean I do not know the sampling rate? And, how should I proceed on to doing a DFT on this dataset?

P.s. Let's say the dataset represents a historical record of the heights of high-tides on an unknown planet with one or more natural satellites. I know that the data points are equispaced, but I don't know the exact value of the spacing.


texts: Software PLL tracking of carrier frequency in bandlimited transmission

texts:
My understanding is that a PLL(either in hardware or software) can recover the carrier frequency(or be used to correct the frequency offset) of a system. My question is how does this work when the transmitter is sending a pulse-shaped signal?

For a an approximate sinusoid the operation of the PLL makes sense but how does tracking work for carriers transmitting data?

As an example, suppose you're receiving a raised-cosine pulse, what do you do?


texts: Explanation for the quantization matlab function

texts:
I found this code on quantization on page 314 in the book Modern Digital and Analog Communication Systems by B.P.Lathi and Zhi Ding. I can understand only half of the code. Please, someone, help in explaining the code.

Code:-

function [q_out, Delta] = uniquan(sig_in, L)
% L = number of uniform quantization levels
% sig_in - input signal vector
% Function Output
%     q_out - quantized output
%     Delta - quantization interval

 sig_pmax = max(sig_in);              % finding the positive peak
 sig_nmax = min(sig_in);              % finding the negative peak
 Delta = (sig_pmax - sig_nmax)/L;     % quantization interval
 q_level = sig_nmax + Delta/2 : Delta : sig_pmax - Delta/2;  %defining q-levels

% Till here I understand the code. I am not able to understand from here.
% Mainly not able to understand what is sigp doing

 sigp = (sig_in - sig_nmax)/Delta + 1/2;   % convert into 1/2 to L+1/2 range
 qindex = round(sigp);                     % round to 1,2,3..... L levels
 qindex = min(qindex,L);                   % Eliminate L+1 as rare possibility
 q_out = q_level(qindex);                  % use index vector to generate output

end

Ebook link - https://www.academia.edu/36135973/Book_Modern_Digital_And_Analog_Communication_Systems_4th_edition_by_Lathi


texts:
I have a signal source that I can digitise and read.  It produces mainly a 50Hz sine wave, but there is distortion as well.  I'd like to estimate the total harmonic distortion for this signal.  I've looked at some on line resources that seem algebraically complicated.  Is there some simpler estimate I can make?  I have the following spectrum of the signal:

and a waveform graph (black line):

I don't have access to any specialised hardware, but I have Scope (sound card based oscilloscope) and Audacity like software...


texts:
Let STFT[frame, freqbin] be the short-time Fourier transform of some audio data.
(here frame is the index in time, and  freqbin is the index in frequency).

Very often, STFT[k, f0] will be close to STFT[k+1, f0]...

... then it could be interesting to store, for each frequency bin f0, only the differences :

STFT[0, f0], STFT[1, f0] - STFT[0, f0], STFT[2, f0] - STFT[1, f0], STFT[3, f0] - STFT[2, f0], etc.

(such processing is known as delta encoding) ; these differences will be small on average, and this could lead to a reduction of the average bit usage for storage of the audio data.

Can this idea be turned into an efficient audio data compression ?

(I don't speak here about psychoacoustic effects that could lead to use different levels of quantizations for each frequency / or frequency masking ; I know that such ideas are used in MP3 compression).

Edit: as mentionned by @pichenettes, the graph above showing continuity between adjacent frames is the "magnitude plot". The "phase / angle plot" is a bit less continuous :


texts: Avoiding clicks with changing biquad coefficients

texts:
I have a DSP application that filters a signal with five cascaded biquads whose coefficients are determined by the Audio EQ Cookbook formulae (Thanks for this work, Mr. Bristow-Johnson btw). Each filter's coefficients can be changed individually via RS232 from a matlab app, however doing so produces audible click noises.

To prevent clicking, my first idea was to lowpass filter the individual coefficients to smoothen the transition. A friend of mine told me that this could cause instability. He is way more experienced in the field than I am, also I still lack the proficiency to prove/-disprove his claim or finding out which circumstances can cause instability in this case, so I take his word for granted.

However, his solution (Smoothing parameters instead and then calculating the correct coefficients) doesn't work for me. Calculating intermediate values in Matlab and transmitting them is not an option because of bandwidth constraints, calculating them on the DSP is not an option because of limited CPU time.

One solution would be to employ two parallel filters and crossfade between them. But I don't particularly like it, because it doubles my computational load, also two changes in rapid succession would lead to audible clicks anyways.

The next idea I had was to employ Gold-Rader structure (see below) for the biquad denominator. Since this structure uses the pole radius, which is always kept below 1 (as long as I filter it with a filter that doesn't overshoot...), it can't become unstable... except for the fact that the radius is multiplied with $sin(\phi)$ and $cos(\phi)$ before introduced in the signal flow graph. So if I'm transmitting the calculated products $r \cdot sin(\phi)$ and $r \cdot cos(\phi)$, I'm again not sure about how filtering those coefficients would affect stability.

Each of the five biquads must be tunable across 20 Hz to 20 kHz, -18 to +18dB gain (where applicable) and a Q from 0.5 to 10 for the Peaks. The type of filter stays the same though (i.e. a highpass will always stay a bypass, a peak filter will always be a peak filter...).

So my question is: what is the most efficient way to avoid clicks while not risking instability of the filter?


texts:
My objective is to build a noise shape filter from a given transfer function (in one case) and from a given PSD (for another case). Checking my precedent questions you can see that this argument is keeping me busy by long time.
You can check this for my first attempt. First, I have noticed that $SD = \sqrt{PSD}$ and not the amplitude spectral density. Secondly I have tried to use the simple approach of putting a white noise in a filter to get a shaped function. My questions are:
Is what I used as source a real white noise?
- I have used the simulink block: Band-Limited White Noise which asks for: Noise power, sample time, seed. I started filling by trial and error procedure to get an averega power spectral density of 1

   prop = 1;
sample = prop*0.1;
pwr_white = prop;                 % autotuning toperforme pwr == 1 unit^2/Hz
seed = randi([1 23341]);

Applying these parameters to the block in simulink and feeding the output of the Band_imited white noise block into an Averaging Power Spectral Density block, the average PSD value is $\pi$ times smaller than the Noise power pwr_white. This is exactly how should it works. So from this information I am supposing that my noise is actually with PSD = 1.
I have tried a similar procedure in MALAB. Using wgn(t_sim/sample,1,0) I get compleatly different results. But var(wgn(t_sim/sample,1,0)) is 1 and var(band limited white noise) is 10.
I don't know where I am wrong.

-Using $H(s)$ as filter will produce in continuos domain the right noise shape?
Looking the aswer cited before, I think yes.
- taking the simulated noise in continuous domain then I need to obtain the PSD (and then SD) to confront with the noise characteristic that I want to perform:
 I came up to the idea that my initial approach was wrong: I was trying to use a noise shape filter simulated in continous to compute than a welch approximation. In this approach I used the MATLAB function pwelch

for i =[1e-5, 1e-4, 1e-3, 1e-2, 1e-1 1]

    f_range = linspace(i,i*10,1000)
    [pxx_noise_t] = pwelch(noise_t,[],[],f_range);
    loglog(f_range,sqrt(pxx_noise_t),'r','LineWidth',1.2)
    hold on

end

The welch approximation is used here in a hamming window (as default) with 50% of overlapping. The f_range given as a vector enables the approximation to be done on the specified range of frequencies. The for cycle is used to put as many points as possible in the approximation, actually the code computes a number of length(i)-1 = 5 welch approximation each in a decade.
Actually I am not sure of the division by $2\pi$ of the noise_t. I have some doubts about the unit of the PSD. Are as deaful in (unit/rad/s)^2?

However, the approach used seems to be not suitable for the welch approximation, or better, for the implementation of a noise shape filter. One main reason is that I am not working in a discrete domain. So I need to built a filter in discrete. This is something that I have never done before, but off course I have never used the pwelch either.

So I need this digital filter. But even if the basic idea of filter is the same for continuous and discrete domain, I need some help. Even in the definition of white gaussian noise (WGN). Actually, checking this answer and the comment I have noticed some problem in applying the approach used in continuos to the case in discrete.
Now I am starting from understanding FIR and IIR filters but I am not sure which one I have to use and how shoul I convert my H(s) in z domain.


texts: Detecting direction of sound using several microphones

texts:
First of all, I've seen a similar thread, however it's a bit different to what I'm trying to achieve. I am constructing a robot which will follow the person who calls it. My idea is to use 3 or 4 microphones - i.e. in the following arrangement in order to determine the from which direction the robot was called:

Where S is source, A, B and C are microphones. The idea is to calculate phase correlation of signals recorded from pairs AB, AC, BC and based on that construct a vector that will point at the source using a kind of triangulation. The system does not even have to work in real time because it will be voice activated - signals from all the microphones will be recorded simultaneously, voice will be sampled from only one microphone and if it fits the voice signature, phase correlation will be computed from the last fraction of second in order to compute the direction. I am aware that this might not work too well i.e. when the robot is called from another room or when there are multiple reflections.

This is just an idea I had, but I have never attempted anything like this and I have several questions before I construct the actual hardware that will do the job:

Is this a typical way of doing this? (i.e. used in phones for noise cancellation?) What are other possible approaches?
Can phase correlation be calculated between 3 sources simultaneously somehow? (i.e. in order to speed up the computation)
Is 22khz sample rate and 12bit depth sufficient for this system? I am especially concerned about the bit depth.
Should the microphones be placed in separate tubes in order to improve separation?


texts: Sampling frequency lower than half of sampled signal frequency?

texts:
In GNURADIO when RTL-SDR Source block is placed there is a parameter Sample Rate (sps) and Channel frequency.

I can put there for example 32kHz into Sample rate and 107.5MHz into Channel frequency to tune some FM station in my location. FFT block shows a signal of that station in a correct frequency as shown on the picture:

But according to Nyquist–Shannon sampling theorem, Sample rate should be at least 2 times higher than signal frequency, so in this situation I should set at least 215MHz into Sample rate parameter.

Why can I set much lower sampling frequency of a signal and it still looks correctly?


texts: RGB, CMYK, HSV - color models or device-dependent color spaces?

texts:
Some sources say that for example RGB, CMYK and HSV are color models instead of color spaces (like CIEXYZ, CIELAB, sRGB and Adobe RGB). Some sources say that RGB, CMYK and HSV are device-dependent color spaces and spaces like CIEXYZ and CIELAB are device-independent spaces.
So, which distinction is better?