texts:
I am implementing the method from this paper:
https://dspace.mit.edu/bitstream/handle/1721.1/66243/Picard_Noncontact%20Automated.pdf?sequence=1&isAllowed=y

The main idea is cardiac pulse measurement using a set of frames ($N=300$) from a 10s video so the frame rate equals $30$ fps.

red = [item[:,:,0] for item in imgs]
green = [item[:,:,1] for item in imgs]
blue = [item[:,:,2] for item in imgs]

red_avg = [item.mean() for item in red]
green_avg = [item.mean() for item in green]
blue_avg = [item.mean() for item in blue]

red_mean, red_std = np.array(red_avg).mean(), np.array(red_avg).std()
green_mean, green_std = np.array(green_avg).mean(), np.array(green_avg).std()
blue_mean, blue_std = np.array(blue_avg).mean(), np.array(blue_avg).std()

red_avg = [(item - red_mean)/red_std for item in red_avg]
green_avg = [(item - green_mean)/green_std for item in green_avg]
blue_avg = [(item - blue_mean)/blue_std for item in blue_avg]

data = np.vstack([signal.detrend(red_avg), signal.detrend(green_avg), signal.detrend(blue_avg)]).reshape(300,3)
from sklearn.decomposition import FastICA
transformer = FastICA(n_components=3)
X_transformed = transformer.fit_transform(data)

from scipy.fftpack import fft

first = X_transformed.T[0]
second = X_transformed.T[1]
third = X_transformed.T[2]

ff = np.fft.fft(first)
fs = np.fft.fft(second)
ft = np.fft.fft(third)

imgs - is the initial list of arrays with 300 image pixel values.
As you can see, I split all frames into RGB channels and thus have traces $x_i(t)%$, where $i = 1,2,3$

After standardization, we detrend all the traces and stack them to further apply ICA and then FFT all the three components.

The method then claims that we need to plot power vs frequency (Hz) and select the component that is most likely to be heart pulse.

  Finally, we applied the fast Fourier transform (FFT) on the selected source signal to
  obtain the power spectrum. The pulse frequency was designated as the frequency that
  corresponded to the highest power of the spectrum within an operational frequency band. For
  our experiments, we set the operational range to [0.75, 4] Hz (corresponding to [45, 240]
  bpm) to provide a wide range of heart rate measurements.

Here's how I try to visualize the frequencies:

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

data = ft
print(fs.size)
ps = np.abs(np.fft.fft(data))**2

sampling_rate = 30

freqs = np.fft.fftfreq(data.size, 1/sampling_rate)
idx = np.argsort(freqs)
#print(idx)
plt.plot(freqs[idx], ps[idx])

What I get is totally different since the range of frequencies is from $-15$ to $15$ and I have no idea whether this is in Hz or not.

The three images above are what I get when I execute the code to visualize frequencies and signal power.

I would appreciate any help or suggestions.


texts: Are SOS structures better for linear phase FIR filters?

texts:
I've seen many discussions suggesting that SOS's are often/usually the preferred structure for fixed point hardware filter implementations, since they preserve conjugate pairs, etc., but I've never seen any specific arguments explaining why they are better for FIR filters.

Do these advantages still apply for direct form, linear phase FIR filters which allow the number of multipliers to be reduced by a factor of 2?

Finally, are there any good resources discussing these trade-offs?  I'm trying to implement fairly long FIR filters with small multipliers, i.e. 40+ taps with 8-10b coefficients, and am looking for resources that might shorten this effort.


texts:
I'm a little new to signal processing and I'm trying to wrap my head around convolutions.

I know the definition of convolution for a continuous signal is

$$y(t) = x(t) * h(t) = \int_{-\infty}^{\infty}{x(\tau)h(t-\tau) \, \mathrm{d}\tau}$$

Let's say for example that

$$h(t)=6e^{-t}u(t)$$

and

$$x(t)=e^{-4t}u(t)$$

you end up getting that

$$y(t)=\int_{-\infty}^{\infty}{6e^{-3\tau-t}}d\tau$$

which is divergent. I watched a few YouTube videos about it and they always explain how to do convolutions graphically with two rectangles, does anybody mind explaining how to do it with two exponential functions, or point me to somewhere that does?


texts:
I have multiple, real, audio signals which I am taking the FFTs. Although they don't come from the same recording, those signals represent the same thing (biologically). I was told that I could combine those ffts and average them together to get a good representation of the averaged signal. I am thinking of using the matlab function 'pwelch' for this purpose but I have lots of questions and I could not find answers by googling.

I have already averaged the ffts themselves with the matlab function 'pwelch' to decrease the noise. When combining them together, should I just use the output of a normal FFT rather than the output of pwelch (P)?
My signals are not exactly of the same length, which will make them not usable to combine in a matrix. Should I cut them all to the same length before averaging? Or is there a way round that?
How do I actually combine them? I was thinking something like this:

:

avg = (abs(P1) + abs(P2) + abs(P3)) / 3;

and then use the function pwelch like this:

[Pavg, Favg] = welch(avg,ones(SegmentLength,1),0,NFFTavg,fs,'power');

Any insights or help would be really appreciated. Thank you!


texts: Why is "$1/\mathrm{Hz}$" unchanged while converting$\mathrm{V}^2/\mathrm{Hz}$ to $\mathrm{dB}/\mathrm{Hz}$?

texts:
When computing power spectral density via Fourier transform and Parseval's theorem, one gets units of $V^2/Hz$ (e.g. Wiki). As a popular next step, the following unit conversion is widely used: $V^2/Hz \to dB/Hz$, and it's done by taking the famous $10.\log10(\cdot)$.

Question: while $\mathrm{dB}$ is a log of a ratio of two $\mathrm{V}^2$'s (second $\mathrm{V}^2$ is a fictive one) and I get where $\mathrm{dB}$ comes from, I'm struggling to understand why $1/\mathrm{Hz}$ remains unchanged after taking the logarithm?

Alex

P.S. Some obscure explanation is in the last paragraph here (and I'm afraid the word 'loosely' conceals all the real stuff here):

  Since the argument of the logarithm is in units of $\mathrm{Hz}^{-1}$, this spectral measure can loosely be said to be in units of `dB/Hz'


texts: Signal to Quantization Noise Ratio (SQNR)

texts:
I have this problem in my homework (following picture).

How do find the signal to quantization noise ratio at the output?

I don't want the solution. I just want to know a better approach to tackle it.

Because the way I am trying to do is -  take the inverse fourier transform of H(z) and convolve it with x(n) and find the power. Then divide this recently found power by the noise power. As easy it sounds (probably correct) it is computationally intense.

I end up with 4 convolutions just to find the y(n).


texts: convolution of a real and complex signal

texts:
How do we define convolution of: $$x(t)= \mathrm{sinc}(W t)$$ and $$y(t)= -j\mathrm{sinc}(2W t)$$

Answer: In the frequency domain, both of them are rectangular functions and multiplication of them is equal to the rectangular signal with smaller bandwidth i.e $x(t)$, which means convoluted signal is $\mathrm{sinc}(W t)$ with some scaling factor.

But I do not understand how this could be done, the two signals are orthogonal (because one is real and the other imaginary). So shouldn't convolution be 0? I would appreciate if someone could explain.


texts: Any Ideas for Extracting the lines on a Football Pitch?

texts:
Ok, so my goal is to extract the pitch bounding lines (not the vertical ones down the center or the semicircle) from an image such as:
source
With the goal of speed, id like to avoid an ML approach, and focus on image processing.
So far, i've been using morphological operations to attempt to achieve my result, the algorithm is as follows:
-> greyscale
-> canny algorithm
-> morphological close (square 3x3 structuring element )
-> morphological open (line 1 x 10 structuring element)
-> morphological dilate x4 (line 1 x 30 structuring element)
I then find the contours of the resulting image, fit a line to each set of contours, and chose the two lines of greatest length to be the pitch lines, and generate a mask accordingly to isolate the pitch.
For this example the results are pretty good : RESULT
But the problem is that it works fairly inconsistently, e.g for the other example image I linked above, here is the RESULT + FINAL IMAGE BEFORE FINDING CONTOURS ... completely misses the pitch. I can alter a few parameters to get it to work for this image, but I'm looking for a solution which can work consistently such that it can be used to process a video, extracting only the pitch each frame.
I think the biggest problem is the crow, lots of noise, lines and colours, leading to miss-classifications or objects as pitch lines.
If anyone has any ideas for a way to consistently detect these pitch lines (morphologically or otherwise) i'd love to hear it.


texts: Fourier transform of a Gaussian in Excel returns imaginary coefficients

texts: I am using Fourier Analysis tool in Excel to transform a Gaussian $\exp(-x^2)$ on a uniform and symmetric grid of $x$ values.  I am expecting the result to contain only real numbers (which should furthermore look like a Gaussian), but what I get is a collection of complex coefficients. I am at a loss. A Fourier transform of a Gaussian should be another Gaussian, why complex coefficients? How to get from them expected Gaussian?


texts: change the central frequency of a wavelet

texts:
I am using the discrete meyer wavelet in MATLAB for time series decomposition. I would like to change the central frequency of this wavelet but I do not know if this is possible or if I would need to design a new wavelet.
Thank you very much in advance.


texts: How to index output of convolution product

texts:
Unfortunately indexing the output of convolution product confuses me. Suppose $x[n]$ is of length $15$ (i.e. $0 <= n <= 14$) and $y[n]$ has the length of $9$ ($0 <= n <= 8$). if $z[n]$ be the output of convolution product of $x[n]$ and $y[n]$, that is $$z[n] = (x * y)[n] = \sum_{k=-\infty}^{+\infty} x[k]y[n-k]$$
As a result, $z[n]$ would has the length of $23$, right?
Now the question is what is the range of $n$ for $z[n]$? Is it $0 <= n <= 22$ or $-4 <= n <= 18$?
From formula for convolution I think the latter is true. But the problem arises when I work with seismic signals.
Suppose $acc[t_i]$ is seismic signal of Kobe earthquake where $0 <= t_i <= 41.99 sec$ with sampling period of $T_s = 0.01 sec$ (length of signal $4200$) as depicted below:

Now it is desired to denoise $acc[t_i]$ using moving average filter with impulse response $h[t_i]$:
$$ h[t_i] =
\begin{cases}
\frac{1}{30},  & \text{if $0 <= t_i <= 0.29$} \\
0, & \text{if $0.3 <= t_i <= 41.99$}
\end{cases}$$
where $T_s = 0.01 sec$. After applying $h[t_i]$ to $acc[t_i]$ by convolution, the output is the sequence $ACC[t_i]$ of length 8399, as shown in following figure.

Obviously the results of filtering process is the interval marked by red bullets. However according to convolution formula, I think I should label time axis between approximately $-21sec$ to $63sec$.
Is that right or something is missing?


texts: How to compute Discrete Cosine Transform

texts:
Can anyone please help me understand the 2D equation of DCT?I find it hard to understand it in terms of how to apply the equation.

For example if I have an 8x8 image and I will apply 2D-DCT to every block of that image where will I start the application from left to right, top to bottom, or after left to right, top to bottom next.???

any help would be appreciated


texts:
I am working on a project in Python to detect and classify some bird song, and I have found myself in a position where I need to convert a wave file into frequency vs. time data. This hasn't been too much of a problem, but to be able to classify the different syllables into groups, I need to write something that will detect when the data clusters into a certain shape. To give you an idea of what the data looks like, here is an image of how the data looks when plotted:

I need some way to get each individual syllable (each shape with a separation on either side) and save them either to a variable or to their own files so that I can run Pearson correlation between them using SciPy.

Also, I prefer Python, but I am open to coding in other languages if you have another way to do it.

Thanks!


texts:
I wish to implement the 2D convolution on an FPGA, so Ineed to understand how it is calculated in practice. The main difficulty that I found apparently 2 different ways showcases how to do it.
The first is to center the kernel to a position and calculate the sum of the overlapping products:
http://www.songho.ca/dsp/convolution/convolution2d_example.html
Making sense of indices in 2D convolution operations

And the other one that involves shifting the kernel
https://www.allaboutcircuits.com/technical-articles/two-dimensional-convolution-in-image-processing/

So what is the way?
Or perhaps they are equivalent in way that I don't recognise?


texts: How does Arbib's clustering method work?

texts:
Recently in Image Processing I've been studying image clustering (or segmentation) methods and one of them is Arbib's. I understand the main idea behind it from which you start with one cluster with its cluster representative then pixel values get assigned to this cluster (also involving a learning rate, what is it?) until you reach a certain amount and the cluster gets cloned, the algorithm stops when you reach the desired amount of clusters.

I searched for quite a while on the internet but i didn't find anything about this method, can you please explain to me precisely how the algorithm works, how pixel values get assigned, what is the learning rate (it might be called in the different way because i'm translating from my language), what are the input and output and all the steps of the algorithm?

Thank you and I apologize for my bad english.


texts: Calculation of line of sight system gain

texts:
I'm trying to calculate the overall gain of the transmitter-receiver system for a line-of-sight wireless transmission.

My question is how to calculate the overall gain of the system. Is it as simple as calculating the gain of the transmitter and receiver separately and then adding them together?


texts: How to Calculate Gaussian Kernel for a Small Support Size?

texts: How to calculate the values of Gaussian kernel? I think I understand the principle of it weighting the center pixel as the means, and those around it according to the $\sigma$ but what would each value be if we  should manually calculate a $3\times 3$ kernel?


texts: Preprocessing audio signal for neural network classification

texts:
I need to identify certain features of the audio signal recorded from microphone in stethoscope. These sounds are only samples i've found, but the final signal will be probably a bit noisier (maybe not, i don't know yet). Below are examples of three signals. Normal one and one of the person sick with atrial spatial defect, third one is of person with late aortic stenosis. These two fluctuations are named S1 and S2 respectively and they are my main subject of interest. Depending on a given heart condition S1 may have doubling and there may be some "noisy" sounds (clicks also) between S1 and S2. The length of the whole cycle (S1, S2) may vary depending on the heart rate.

My question is, what kind of processing should i perform to surely extract and recognize those features, i.e. two large fluctuations normally, with slight delay between them and if person is sick it may be either something before S1 or between S1 and S2 or S1 may double a bit. Also normally S1 has a larger amplitude spike than S2 as i understand.

I have very limited knowledge in this field and not sure about what options are available to me. FFT magnitude graphs are showing somewhat recognizable results, but i think there is too much unneeded information there and they are not so distinct.

The samples themselves can be found here: Washington University: Heart sound samples

Normal heart sound:

Atrial spatial defect heart sound:

Late aortic stenosis:


texts: Difficulties acquiring sample period of a signal

texts:
I'm teaching myself some DSP for future projects and my method is mostly learn by doing. I have this time-domain function with its correspondent frequency-domain transform:
$$s_1(t) = (t-2)e^{-t}u(t-2) $$
$$S_1(f)=\frac{e^{-2(1+j2\pi f)}}{(1+j2\pi f)^2} $$
Plotting the frequency-domain function in Matlab, for a frequency vector ranging from [-50,50) in steps of 0.1 (1000 values of frequency).

The question I don't know how to solve, is about the period sample. I know each sample is one of those dots, and that despite it can be arbitrary (choosed by us), there is criteria to follow in order to pick the best sampling rate possible to avoid aliasing.
If I set the frequency vector like I did, that means I already determined the sample rate, am I right? That also means, I can compute the sampling period as:
Ts=1/fs; fs=0.1 being the sample frequency.
Or should I go with another one? Should I try Nyquist (since its a continuous signal)? If so, how could I apply it?
Thanks for any suggestion you may have on this.


texts:
Accidentally asked this question in the general area and was told to ask here, so...

I've been trying to develop a lightweight, relatively-fast-to-decode sound compression format for use in my gaming projects (perfect reproduction is not needed so I only use 16 bit data).

The idea is to split the sound data into 14-sample frames and use linear prediction to reduce the data size by only storing the residuals. To make it even lighter, the residuals are then quantized to 4 bits per sample by reducing their precision so that they are stored as scaled residuals, where the scale is dictated by the frame header. To make the signal less noisy, 16 linear prediction models are generated that best suit the signal in the file.

Each frame ends up being 64 bits (4 bits for the LP model selection, 4 bits for the residual scale (2^n, so only n is stored), and 14*4 bits for the residual data).

What I've been doing is simply solving the Yule-Walker equations for each frame to get the coefficients and then saving them in an array. After calculating the coefficients for all the frames, I then try to create 16 coefficients based on the Euclidean distance from each LP model and averaging them out. That is to say...

Zero-out target buffer

For i = 1, i <= 16, i++
  Create and zero-out a temporary buffer for i LP models

  For each saved LP model
    Find the model in the original target buffer [up to i] with the smallest Euclidean distance between it and this model
    Accumulate the coefficients to the temporary buffer

  Average out the temporary buffer and store it to the final LP model array [again up to i]

As expected, this doesn't yield very good results at all if only because I imagine it's similar to uniform colour quantization - it has very little relation to the actual data.

After Googling for a few months, the only thing I can come up with is 'sparse convolution' but I'm not entirely sure what is meant by that when it comes to LPC or if it's even what I'm after.

Given these parameters (that is: generate 16 LPC models that will minimize the residual error given that the frames are 14 samples long), how would you go about doing it?


texts:
Please recommend books for learning/studying signals and systems
I know about oppenhiem,please write about other books that are easy to understand


texts: Scaling of the PWELCH function in MATLAB

texts:
I am trying to compute the power spectral density of a random signal using the PWELCH function in MATLAB.

Since I think have not understood properly how Pwelch scales the PSD, I wrote a sample program, in which I generate a sum of two sinusoids with given amplitudes - A1 and A2 - and given frequencies - f1 and f2 - superimposed to random noise.

I computed the PSD using the FFT and then using Pwelch. It seems to me that the PSD computed via FFT is correct, because it returns the power corresponding to the sinusoids (peaks at $A1^2$ and $A2^2$).

However Pwelch seems to have a different scaling because the peak values depend on the number of samples I put in the PWELCH window. This should make sense up to a certain extent, but the variation seems too large to me. I think I haven't properly understoon either how to use PWELCH or how to scale the PSD returned by PWLECH.

Can you help me understand this? Thanks. The sample code is reported below.

Thank you
E.

clear all
close all
clc

rng default

Fs = 10000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 50000;             % Length of signal

t = (0:L-1)*T;

A1=2;
A2=3;

f1=5;
f2=20;

X=A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t)+ 2*randn(size(t));

figure(1)
hold on
grid on
plot(t,X)
xlabel('time (s)')
ylabel('signal (Unit)')

%% Power Spectrum via FFT
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L

Y = fft(X);

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

%% via FFT

figure(2)
hold on
grid on
plot(f,P1.^2)
title('Single-Sided Power Spectrum of X(t)')
xlim([0 50])
xlabel('f (Hz)')
ylabel('PSD (Unit^2/Hz)')
set(gca,'xscale','log')
% set(gca,'yscale','log')

%% via PWELCH

nfft= L;
window= rectwin(nfft);

[Pxx_f,f]=pwelch(X,window,[],nfft,Fs);

figure(2)
plot(f,Pxx_f,'r','linewidth',1.5);

legend('via FFT','via Pwelch','location','best')


texts: Breaking down sin wave to frequency domain causing weirdness in phase

texts:
I have tried many many different things and I cannot figure out why my phase is messed up. Here is what I am trying to do. I create a simple signal with a known wavelength, phase, and amplitude like below...
time_vec = np.arange(1, 11, 1)
wavelength = 1/.2
phase = np.pi
amp = 10
created_equation = amp * np.sin((2 * np.pi / wavelength * time_vec) + phase)

I then do an FFT to get individual frequencies. I end up with 4 frequencies from the FFT but only 1 at 2 hz is of any significance, which is exactly as expected. From here I am trying to get the wavelength, phase, and amplitude of this dominent frequency, create a new equation from these, and do inverse FFT to get back to my original, known, signal. Everything works except the phase is acting strange. It is producing mirror images of the expected result. I am using arctan2(Im, Re) for phase but I have also tried arctan(Re / Im) and get the same result. Here is the plot at 0 phase, and 1hz known signal.

and another at 0 phase, and 2 hz....

I can try any combination of phase and frequency and it is always slightly off. I have tried to manually shift the phase and this will solve the issue for any frequency at the same phase, ex. phase=0 and freq=any. But as soon as I change the phase, it breaks again so this is not a valid solution. What is causing this? and how do I fix it?


texts: Suggestions for inhomogeneous particle segmentation

texts:
I need to detect multiple particles of different colors in an image (see example image). In case of very homogeneous particle coating distributions, i.e. all particles look very similar, I can achieve good results by template matching with a synthesized instance of an ideal particle.

However, when I have instances of particles that are white at the same time of particles that are fully coated, this matching approach does not yield good results, especially the dark particles are detected with a smaller probability than the white ones.

After segmenting the particles, a classification is done using a CNN (which works well) - however my accumulated results are negatively affected by the different detection probability of coated and uncoated particles.

I would be very happy to hear some recommendations about what algorithms can be used to segment the particles reliably. I assume those would be based on the particle's shape.

Sample Image

Synthetic Particle Template

Detected Particle Positions (other base image)


texts: What are the most suitable methods for correcting sudden DC offsets or baseline offset? and Why? See the image for the use-case or example

texts:
The signal values should not get affected by the correction. Most methods I found were related to baseline drift or exponential baseline shift/decay. I am unable to find a method suitable for my use case.
In the image it is visible that after sudden dc offset the baseline got shifted. After correction it should stay same.


texts:
I'm trying to understand Wavelet Transform , I've done well enough to understand Continuous Wavelet Transform which was easy enough, where we simply stretch the wavelet and match it with the original signal.
But i couldn't be able to understand the Discrete Wavelet transform in fact I've always had a trouble understanding Discrete of any transform (like DFT , DCT which are still not clear to me).

According to Wavelet transform we Stretch our signal and then match it with the original signal to find out the frequency.
This is the example shown in the book,

  In this hypothetical example the student does fairly well the first
  half of the term then neglects his or her studies for the last half.
  Thus the exam scores for the term were 80%, 80%, 80%, 80%, 0%, 0%, 0%,
  and 0%* We can tell the average of all the scores (40%) and when the
  scores “tanked” after the 4th exam just by looking. Knowing the answer
  in advance, however, is a good way to learn and to verify the wavelet
  transforms. Then we can use them with confidence on real-world data
  where we can’t simply “eyeball” the final values.

  We will now walk through the CWT process step by step using the
  simplest of the wavelet filters on this example. We begin by comparing
  the humble Haar wavelet filter, [1 –1]

Signal -> [80 80 80 80 0 0 0 0]
Filter(or signal that will be stretched later) -> [1 -1]

Comparing the first 2 points with the wavelet filter we obtain 80 – 80 = 0. For
this very simple high-pass filter we can say there was no change in the first 2
exam scores.

This is how we've done it,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 80 80 80 0 0 0 0]

Shifting once to the right and applying filter again we get,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 0 80 80 0 0 0 0]

Shifting once again the right and applying filter again we get,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 0 0 80 0 0 0 0]

Shifting again,

80*1 + 0*(-1) => 80 - 0 = 80

Now We have,

[0 0 0 80 0 0 0 0]

We will eventually end-up with

[0, 0, 0, 80, 0, 0, 0]

This is significant in that this wavelet process of comparison and
shifting has just indicated a large change between the 4th and 5th exam. We
have “found the discontinuity”.

This actually makes sense as we had a large change from 4th to 5th point But how does this discontinuity defines a frequency here ?

Now if we stretch our signal from [1 -1] to [1 0 -1] and repeat the process,

We will eventually end-up with,

= –80 –80 0 0 80 80 0 0 0 0

When we used filter [1 -1] , our output at least makes some sense it clearly shows a rapid change, but what happened here in case of 3 point filter ? what are these weird values (-80 -80) ?

If we keep scaling to 10th times here is what we'll eventually get,

(scale = 10) –320 –160 0 160 320 320 240 160
(scale = 9)  –240 –80 80 240 320 240 160 80
(scale = 8)  –320 –160 0 160 320 240 160 80
(scale = 7)  –240 –80 80 240 240 160 80 0
(scale = 6)  –240 –160 0 160 240 160 80 0
(scale = 5)  –160 –80 80 160 160 80 0 0
(scale = 4)  –160 –80 0 80 160 80 0 0
(scale = 3)  –80 0 0 80 80 0 0 0
(scale = 2)  –80 0 0 0 80 0 0 0
(scale = 1)   0 0 0 0 0 0 0 0

While studying the Continuous Wavelet Transform the scaling made sense but i could't understand the scaling here ? what are these points -80, -160, -240, -320 represents ?

This is the final output,

How can i relate this output from my original signal ?


texts: Image processing based on wavelet transform in python

texts:
I take the dwt2 for an image and saved it's coefficients  (LL,Lh, hl,hh) using pywt.dwt2 (image,'haar'),in my project I have to change them to uint8 but when I change their types and reconstruct the image an error appear prevent the show of the result image .can any body help me with this problem! My code is like:

import cv2
import numpy as np
import pywt
image=cv2.imread('1.tif')
image=cv2.resize(image,(256,256))
[b,g,r]=cv2.split(image)
x=b
coffes=pywt.dwt2(x,'haar')
ca,(ch,cv,cd)=coffes
ca2=np.uint8(ca)
ch2=np.uint8(ch)
cv2=np.uint8(cv)
cd2=np.uint8(cd)
coffes=ca2,(ch2,cv2,cd2)
origin=pywt.idwt2(coffes,'haar')
cv2.imshow('image',origin)
cv2.waitKey(0)
cv2.destroyAllWindows()


texts: Play with a Gaussian Random Set in the Frequency Domain to Obtain Desired Effect in the Time Domain

texts:
Please assume we have a set of 100 random numbers obeying Gaussian PDF in time-domain. Let's index them 1-100.

Within a time-accurate simulation, I apply four operations on this dynamic set every time-step:

Use the entire set for some mathematical operation
Then, remove, say, the last 10 numbers of the set (indices:91-100)
Shift the remaining 90 numbers forward from the indices:1-90 to 10-100
Generate 10 new Gaussian random numbers and put them into the emptied indices:1-10. And the entire set should still obey Gaussian PDF.

I want to obtain these time-domain changes via performing some operations in the frequency domain.

However, I do not know how I should play with the frequency-domain counterpart of the set, so that I obtain the same (or similar) results as if performing time-domain operations above.

The obscure form of the question:

One can generate random samples of Gaussian distribution directly in the frequency domain as explained here, and exemplified with a Python snippet here.

Let assume we generate a one-dimensional set consisting of uniformly-spaced 100 Gaussian random samples, $X(f)$.

I want to manipulate the set as such its time-domain counterpart, i.e. $x(t)$, is:

Index-shifted by N positions in a circular manner,
Then, (only) its first N elements being replaced by new N random samples.

For the first item, this gives the answer. For the second, however, I am perplexed.

Could you please tell us how to insert this new N random samples directly into the existing set in the frequency domain, so that we can avoid possibly(?) redundant FFTs otherwise?


texts:
I'm plotting the FFT power-spectrum of a signal in MATLAB.
I uploaded the 8000 samples time-series signal in a text file here:
http://wikisend.com/download/896484/signal.txt

I'm using the following code:

Fs = 500; % sampling frequency
L = length(y); %number of samples

complex = fft(y)/L; % complex signals
f = 0 : Fs/L : Fs/2; % frequency bins
amplitude = 2*abs(complex(1:L/2+1)); % amplitudes
pow = (amplitude).^2/2; % power

semilogy(f,pow,'-ro');
grid on;

And I'm having the following plot:

The confusion is, where I work they are using their own modified FFT function and if I run that code I obtain the following plot:

If you check the power in above figures the linear power ratio is 16 between my code and their's. And for other measurement I'm getting the ratio 60, 120 or even 1.

I'm confused and need help if one can plot FFT for me and see if my code and plot correct.


texts: Any way to compensate for fractional delay filter's weakness between taps?

texts:
A little background: I have an FPGA (running at 200MHz) that I am creating multipath events using buffers, all of that works.  I take RF, mix down to IF, sample, make adjustments, out the DAC, and back up to RF.  I use a windowed-sinc algorithm (like described here).

To do fractional delays, I compute the coefficients in a Java app and pass them to the FPGA.  Each coefficient is a separate tap in the path on the FPGA, and then everything is summed up before outputting.

Like I mentioned, this all works, but I would like to make it work better.  What I see is that for delays that work out to be closer to the middle of an FPGA clock "bin", I get the worst results.  When the delay is closer to one of the edges of a bin I get better results.

As a test, I have 2462MHz for my RF, 2437MHz for my LO, and get an IF of 25MHz (just using pure tones out a signal generator).  I have the tone going straight through on one tap, and I sum it with a fractional delay of itself.  I looked at the output signal level on a spectrum analyzer for the first 15 completely out of phase points and plotted the results (in an ideal world they would cancel each other out).  In this image the very first plot is actually where the second tap has a delay of 0 (completely in phase) to give a point of reference:

As you can see, the out of phase signals that fall near the beginning and end of a bin perform pretty well (and seem totally acceptable to me), the ones that fall near the center of a bin have the hardest time cancelling out.  It might be acceptable (and as good as I can do), but I am curious if there is a way to perform better.

As an FYI, if I look at points within the first 2 periods of the RF signal (not just the 180 degree out of phase points) at 0.05ns increments, it looks very nice and smooth, the only thing that catches your eye is that there is a slight rise to the output, but that can be attributed to the fact that as I go along, I am getting closer to the center of a bin, and the results are getting worse.


texts: Concept About Estimated Standard Deviation

texts:
I am looking for the concept about how to estimate standard deviation.
Actually I'm not sure how can I get a concept the estimate standard deviation ?
If you know the concept, then would you let me know it here?
Also if you have any reference or examples, it will be help to me.
p.s. Is this different between "estimate standard deviation" and "estimated standard deviation" ? I think it seems to the same.
Updated
What does different between "Method of Least Squares" and "estimate standard deviation"?
Are those the same concept? I don't know estimate standard deviation.


texts:
I'm trying to segment (or extract) different calls from a spectrogram. My current method, calculates a threshold value (variance) and then using RMS can detect voiced/unvoiced signals:

$$ Thres = (1 - c)  * min(x) + c * max(x) $$

Where $c = 10$

This method works ok, for most signals. However, does not produce high results when using different data-sets. I have read this article (http://www.avisoft.com/classifying%20recordings.pdf) which uses Frequency contour cross-correlation in order to classify the different calls.

I'm wondering whether or not someone could provide an insight into this, or at least provide some research based material to which I can read and hopefully gain some kind of understanding.

Thanks


texts:
In wikipedia, there was an example on QPSK modulation where the main signal was 11000110 and I contained 1001 and Q contained 1010, now their output was shown by the Signal in the figure.

Now according to constellation diagram Signal can be represented as

11 -> pi/4
00 -> -(3*pi)/4
01 -> (3*pi)/4
10 -> -pi/4

But my confusion was in the phase vs time diagram shown in another figure given below

For 11 it starts from phase pi/4, but for 00 shouldn't the phase start from -(3*pi)/4 = -135 ?


texts:
MATLAB's icwt docs state inversion to be done by a single integral:
$$
f(t) = 2 \Re e\left\{
   \frac{1}{C_{\psi, \delta}} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \tag{1}
   \right\}
$$
Where does this expression come from, and how is integration over $b$ (translation) eliminated (as in Wiki's formula)? Does it hold if computing discretely over log-distributed $a$ (e.g. $(2^{1/32})^{[1, 2, ...]}$)?


texts: How to determine sampling frequency for an x(t) signal avoiding aliasing?

texts:
The antitransform of the function is given:

$ \hat{x}(f) = \frac{123 + i246\pi f}{246 - 24600 \pi^2 f^2 + i4920\pi f} $

I'm asked to determine which frequency can I sample using the function x(t) avoiding aliasing? (or which range of values)

Update the problem actually suggest five answers:
100, 10, 1, 0.4 or none of these. Turns out the answer is none of these.


texts: Sketch the following Signal: $x(t)=3\sin(\frac{\pi t}{4})$

texts: I feel kind of overwhelmed by it, and would like someone to walk me through the steps of sketching the signal.


texts: What kind of math is introduced in this answer related to DSP?

texts:
What kind of math is introduced in this answer:

https://dsp.stackexchange.com/a/40225/30897

It smells like standard linear algebra course (which I haven't studied, yet), but is it only that?
Thanks!


texts: Adaptive Gaussian Filter for Image Denoising

texts:
I am looking for methods to enhance noisy images, where:

some pixels in the image are very noise,
some other pixels do not contain so much noise.

My first thought is to build an adaptive Gaussian filter.
This means that the Gaussian kernel will depend on the (estimated) noise status of the current pixel (large radius of the Gaussian kernel for noisy pixels and a smaller in the converse case).

Could you help me in calculating the kernel values?


texts: Digital filter simulating hydroacoustic signal distortion

texts:
Preambule:
I'm designing a sound model for my small submarine game. Model is running on the server, and I want to present the client with a mono-channel wav-stream from his hydrophone (20kHz discretization should suffice, I target 20Hz-10kHz band). I want that signal to be relatively-realistic. Specifically, I need to account for the fact that water dissipates high frequencies much faster.

Sound model:
I plan to pre-generate (with help of some RNG and IFFT) raw propeller time-domain sound sample from it's intensity spectrum. Then I plan to amplitude-modulate it with shifted shaft rpm x propeller blade count sine wave to get that propeller beat. At this point i need to materialize that signal on the player's sensor. It will probably involve two filters:

Some filter (wich is the filter in question) wich will account for water dissipating signal energy.
Bandpass filter, wich corresponds to this particular sensor's sensitivity band.
Simple amplitude multiplication to account for range, sensor directivity etc.

According to R.J. Urik "Principles of Underwater Sound", on the distance of r meters from the target, band level (dB) in frequency bin f of such sound can be approximated using following formula:

[1]: ResultBL = SourceBL - 10 * log10(r * r) - r * F(f)
or without wave front expansion term (wich is trivial):
[2]: ResultBL = SourceBL - r * F(f)

where F(f) is some smooth monotonic function (big radical with multiple squared frequencies and constants) of frequency.

The questions (all about this one filter, or lack thereof):

Is it possible to design a time-domain digital filter for given range r, that implements\approximates the signal distortion [2]?
What type of filter would you recommend? Problem is not realtime, but calculations should be fast. I'm not limited to causal filters, so I'm sensing something like non-causal IIR?
What algorithms are applicable to this problem.
If such algorithm involves transfer function, how should I transform\approximate F(f) term in it?
Will it be fast? What factors affect it's speed?
Will such filter work good in the whole band (20Hz-10kHz)?
Will such filter work well on large intensity level variances, e.g. on both very weak and very strong signals?
Is such synthesis computationally-expensive?

Thank you.


texts:
I am running into a problem where I have a control system $S[t]$ that takes a control $C[t]$, so that
$$S[t+1] = H(C[t,t-1,...], S[t,t-1,...])$$ the response of the system is the history of controls and responses.

I am looking for a way to control $C[t]$ to optimize the value of $S[t]$ at any point in time.

What is the literature on this problem called? Where do I start looking? Thanks.


texts: What are advantages of having higher sampling rate of a signal?

texts:
Being a non signal processing science student I have limited understanding of the concepts.

I have a continuous periodic bearing faulty signal (with time amplitudes) which are sampled at $12\textrm{ kHz}$ and $48\textrm{ kHz}$ frequencies. I have utilized some machine learning techniques (Convolutional Neural Network) to classify faulty signals to the non faulty signals.

When I am using $12\textrm{ kHz}$ I am able to achieve a classification accuracy $97 \pm 1.2 \%$ accuracy. Similarly I am able to achieve accuracy of $95\%$ when I applied the same technique on the same signal but sampled at $48\textrm{ kHz}$ despite the recording made at same RPM, load, and recording angle with the sensor.

What could be the reason for this increased rate of misclassification?
Are there any techniques to spot differences in the signal?
Are higher resolution signals prone to higher noise?

Details of the signal can be seen here, in chapter 3.


texts:
I've been trying to work out the logic for this task, and plan to use the KissFFT source package to perform the fast fourier transform.  Please let me know if this looks about right:

Allocate an FFT structure, ie. kiss_fft_alloc(N,0,NULL,NULL)
Where N is the window size I'm using.  The input buffer will be an array of N elements of type kiss_fft_scalar.  The output buffer will be an array of N/2 + 1 elements of type kiss_fft_cpx.
Decode N (window size) number of PCM samples.
For each PCM sample, average each channel's amplitude (unsigned samples) and scale from 0 to 2 (divide by 65536.0), storing the result into the input buffer.
Perform windowing (ie. Hanning) on the input buffer.
Perform fast fourier transform on the input buffer, storing into the output buffer.  Since I am using real values as input, I can use kiss_fftr().
For the N/2 output values, obtain the squared magnitude of the transformed data and convert the values to the dB scale with the following formula:
10 * log10 (re * re + im * im)
Plot the N/2 values from step 6.
Discard the first half of the input buffer, decoding the next (window size / 2) PCM samples and performing scaling and windowing to the data.  This should effectively slide the input window and avoid having to redo math on processed PCM samples.
Loop to step 5, repeating these steps until all samples are processed.
Free the used memory from kiss_fft_alloc().

It was suggested that I subtract a value from the input window before I perform the FFT, so that the resulting DC value has a magnitude of zero.  Should I be subtracting the mean or the average from the input data?

Also, what are the things I need to consider when I choose a window size?  Besides that it has to be an even number as per KissFFT's instructions, is there a benefit to using a small window size, ie. will it provide for a better graph?  I assume that a large window size reduces the number of FFTs that must be performed, is that the only benefit to using a large window size?

Lastly, when I get to the point that the data is ready to plot, how do I go about plotting it?  When I worked on some waveform graph logic in the past, I just plotted 3 values for each pixel along the $x$-axis (min amplitude, max amplitude, RMS amplitude), but I don't know what I'm supposed to do with spectrogram data.

Thank you in advance for any and all guidance you can provide.


texts:
I need to preprocess raw ecg data in R, here is a sample already standardized. I'm not an expert in signal processing nor experienced in working with medical data, so I need concrete answers ideally with code.

Below is a plot of how the data looks like, so perhaps this can give an idea of the problems present.

Question 1)
What strategy is needed in terms of subsampling, filtering to preprocess the data to clean get a clean ECG signal (google how it should look).

Question 2)
Would this strategy work for EMG data also?


texts: What algorithm generates these "bucket coloring" artifacts?

texts:
Context: This is picture (taken with a smartphone*) of a copy of school test. The original answers were in black pencil and the professor made check marks in red.

However, the parts of the answer that touched the red check mark, became red as well.

Looking up close at the copy it self, it seems to have been printed by a laser printer, given the small but visible gaps within the red check mark.

*The artifacts were NOT generated by the phone. They were present in the printed copy itself. Thus, either the scanner or the printer is to blame.


texts: Convolving two complex signals - relationship between phase

texts: Suppose two continuous complex domain signals are convolved, how is the magnitude and phase of the resultant signal related to the magnitudes and phase of the original signals?


texts: Resample an array in real-time (application: pitch-bend)

texts:
I'm working on a sampler software (that plays .WAV files, when notes are played on a MIDI keyboard).

In order to implement pitch-bend feature (i.e. you have a "pitch bend wheel" on the synth and you can turn this wheel to retune the sound in realtime, classical effect used in funk music),

we have to:

Retune: this is possible by resampling the array of the .WAV.
With Python, there is scipy.signal.resample or probably the same in numpy.
But I heard there are algorithms optimized for audio repitching. What are the audio-optimized resampling algorithms?
[Secret Rabbit Code : http://www.mega-nerd.com/SRC/. Why is this important? How to implement it in Python, with numpy for example?]
Retune in realtime: indeed, we can play with the pitchbend-wheel very fast, thus, hundreds of resampling factor could be needed per second! How to do that?
How to resample an array with different resampling factors, in realtime?


texts: Is this a correct method for fixed pattern noise correction?

texts:
I am currently involved in a project involving programming an imaging sensor. Our sensor is giving us noise, so we want to correct for it. Someone else on the project came up with the idea to take a "black" image i.e. put the lens cap on and take an image that is meant to be all black. (Obviously it is not due to noise) At this point on subsequent captures he takes the pixel values from the black image and subtracts them from the regularly captured image.

The image does look better and most of the noise is removed, however I am not convinced that this is the best approach for removing the noise due to the following:

The range of the fixed image is [-172 194] (366 values), versus the standard range of [0 255]. When it is redrawn it gets ranged back to [0 255], and it does look better, however I believe this is incorrect.

I should mention that the new image is taken in low light.

Is this method correct for removing the noise? Why or why not?


texts: Given two waves, 441 Hz and 700 Hz, how many samples long is each wave's period if the sample-rate is 48 kHz?

texts:
The title says it all.

I have two waves, Master and Slave.  Master is at 441 Hz and Slave at 700 Hz.  The sample-rate is 48 kHz.

I want to know how to calculate the period of each wave in terms of samples.  How many samples long is each wave's period?

Please show me how to calculate this

Thanks


texts: What are the known use cases for different wavelet families?

texts:
There are numerous wavelet families that differ by a number of parameters (like the number of vanishing moments, symmetry, etc.).

What are the known use cases where one particular family is known to be more suitable and provide better results?
For example, Daubechies 9/7 are used in the JPEG2000 image format. As far as I understand Daubechies wavelets were designed specially to have a high number of vanishing moments so that makes them more suitable for signal compression tasks. What are the known use cases for other families of wavelets?

And as second part:

How does one decide which one is the most suitable for it's task?
Are there any known heuristics, common knowledge, etc?


texts: Channel weighting of Mono channel in EBU R128

texts:
I am developing a loudness meter conforming to the EBU R128 recommendation.
Currently, it works with stereo files but not with mono files.

In the ITU Recommendation, where the metering algorithms are defined, weightings for the individual audio channels are given (table 3, page 6). From that, I assumed that the weighting for a mono channel is 1.0 (0 dB).

As you can see in figure 1 on page 3, after applying k-filter, mean square and channel weighting to each individual channel, the resulting values of each individual channel are summed to one value. Logic dictates that when having two channels of the same audio signal, the value will be higher than when just having one.

As a consequence, a mono file will result in short-term loudness values, that are 3 LU lower in comparison to a stereo file with same content (on both L and R channels). This is clearly wrong, because established meters like "Klangfreund LUFS Meter" and "MeldaProductions MLoudness Analyzer" will output the same short-term loudness with both files.

How am I interpreting the algorithm definition wrong? How to deal with mono signals?

Side note: It seems that FFmpeg once had a very similar issue, but I could not find the solution to it.


texts: Hi All: I was thinking about buying a "DSP book with matlab". I'm not sure which one is generally recommended  but that's not my question. My question is, when you buy such a book, do they require the DSP toolbox (assuming it exists) or just the regular base matlab. I have the base matlab but I would think the  toolboxes are pretty expensive so I don't want to buy a book only to find out that I need  the dsp toolbox. Or if the answer is that there is no such toolbox, that's fine also. Thanks.


texts: Basics of Normalizing Cross-Correlation with a View to Comparing Signals

texts:
I'm trying to understand how cross-correlation is used determine the similarity of two signals. This tutorial offers a very clear explanation of the basics, but I still don't understand how to use normalization effectively to prevent strong signals from dominating the cross-correlation measure when you have signals with different energy levels. The same tutor, David Dorran, discusses the issue of normalization here, and explains how to normalize the correlation using the dot product, but I still have some questions.

I wrote this python routine to do cross-correlations between every pair of signals from a group of signals:

import numpy as np
import pandas as pd

def mycorrelate2d(df, normalized=False):
    # initialize cross correlation matrix with zeros
    ccm = np.zeros(shape=df.shape, dtype=list)
    for i, row_dict1 in enumerate(
                        df.to_dict(orient='records')):
        outer_row = list(row_dict1.values())
        for j, row_dict2 in enumerate(
                            df.to_dict(orient='records')):
            inner_row = list(row_dict2.values())
            x = np.correlate(inner_row, outer_row)
            if normalized:
                n = np.dot(inner_row, outer_row)
                x = x / n
            ccm[i][j] = x
    return ccm

Suppose I have 3 signals of increasing magnitude:
[1, 2, 3], [4, 5, 6] and [7, 8, 9]

I want to cross-correlate these three signals to see which pairs are similar, but when I pass these 3 signals into the routine I wrote, I don't appear to get a measure of similarity. The size of the cross correlation values is just a function of the energy signal. Period. Even the cross correlation of a signal with itself yields lower values than the cross correlation of that same signal with another signal of higher energy.

df_x3 = pd.DataFrame(
        np.array([[1, 2, 3],
                  [4, 5, 6],
                  [7, 8, 9]]).reshape(3, -1))
mycorrelate2d(df_x3)

This yields:

array([[array([ 3,  8, 14,  8,  3]),
        array([12, 23, 32, 17,  6]),
        array([21, 38, 50, 26,  9])],
       [array([ 6, 17, 32, 23, 12]),
        array([24, 50, 77, 50, 24]),
        array([ 42,  83, 122,  77,  36])],
       [array([ 9, 26, 50, 38, 21]),
        array([ 36,  77, 122,  83,  42]),
        array([ 63, 128, 194, 128,  63])]], dtype=object)

Now, I pass in the same 3 signals, but this time I indicate that I want normalized results:

mycorrelate2d(df_x3, normalized=True)

This yields:

array([[array([ 0.2142, 0.5714,  1., 0.5714, 0.2142]),
        array([ 0.375,  0.71875, 1., 0.5312, 0.1875]),
        array([ 0.42,   0.76,    1., 0.52,   0.18])],
       [array([ 0.1875, 0.5312,  1., 0.7187, 0.375]),
        array([ 0.3116, 0.6493,  1., 0.6493, 0.3116]),
        array([ 0.3442, 0.6803,  1., 0.6311, 0.2950])],
       [array([ 0.18,   0.52,    1., 0.76,   0.42]),
        array([ 0.2950, 0.6311,  1., 0.6803, 0.3442]),
        array([ 0.3247, 0.6597,  1., 0.6597, 0.3247])]],
        dtype=object)

All the max values are now 1!! So we went from having max values that were based on spurious differences to having no difference between the max values at all! I readily confess that I do not understand how cross-correlation is used to detect similarity between signals. What is the analytical work flow of someone comparing signals with cross-correlation?


texts: MATLAB: $\tt fft$ and $\tt ifft$ scaling

texts:
In MATLAB, the outputs of the fft and/or ifft functions often require additional processing before being considered for analysis.

I have heard many differing opinions on what is correct:

Scaling

Mathworks states that fft and ifft functions are based on the following equations:
\begin{align}
X[k] &= \frac{1}{1} \cdot \sum_{n=1}^{N} x[n] \cdot e^{\frac{-j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}}, \quad\textrm{where}\quad 1\leq k\leq N\\
x[n] &= \frac{1}{N} \cdot \sum_{k=1}^{N} X[k] \cdot e^{\frac{+j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad \textrm{where}\quad 1 \leq n\leq N
\end{align}
Scaling by signal length

My peers typically scale the data by $\small \frac{1}{N}$ immediately after the processing the fft.
(We do not consider the raw fft data before scaling.)
%% Perform fft
X_f = fft(x, n_sample, 1) / n_sample;
% fft must be normalized by the number of samples in the data.
% This convention was set by the software developer (Mathworks).

Is this correct?

If so,  why does the MATLAB ifft function expect that we have not scaled by $1/N$ already?
Is there a MATLAB ifft function or function option which does not automatically scale by $1/N$?

Alternatively, is there a better convention which we should be using in placing the $1/N$? For example, placing the $1/N$ in the fft rather than the ifft, or placing an $1/\sqrt{N}$ in both equations, instead of an $1/N$?
Scaling by sampling period

I have heard that the fft and ifft functions assume that the sampling period $T_{\rm sampling} = 1/f_{\rm sampling} = 1$, and that for the functions to be true, the following would need to apply:

\begin{align}
X[k] &= \frac{1}{T_{\rm sampling}} \cdot \sum_{n=1}^{N} x[n] \cdot e^{\frac{-j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad \textrm{where}\quad 1 \leq k \leq N\\
x[n] &= \frac{T_{\rm sampling}}{N} \cdot \sum_{k=1}^{N} X[k] \cdot e^{\frac{+j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad\textrm{where}\quad 1 \leq n \leq N
\end{align}

See links:

Link 1 (see comment to Matt Szelistowski by Dr Seis)
Link 2 (see answer  by Rick Rosson vs that of Dr Seis)
Link 3 (see comment by Matt (Message: 7/16) and comment by Poorya (14/16)
Link 4 (see pg. 10, slide [1,1])
Link 5 (see pg. 8+9) [it seems he is using inverse convention for fft and ifft].

Is this true?

I'm particularly piqued because I cannot find any DFT or DTFT equations on Wikipedia which include the sampling period.


texts: The concept of poles and zeros in filters was introduced to me during our lab (our lab isn't sync with our lecture) through the pole-zero plot generation of filters in MATLAB. I didn't get its importance until it was mentioned in our lecture that poles make the equation infinity (or its denominator zero) while the zeros make the roots in the numerator zero. I still can't quite grasp what poles and zeros are. Can someone give an an intuitive and simple explanation on what poles and zeroes are? Much appreciated.


texts: Difference between freqz vs freqresp in scipy?

texts:
I'm always confused that scipy has

scipy.signal.freqz, but also
scipy.signal.lti.freqresp / scipy.signal.dlti.freqresp

According to the documentation they all compute the frequency response. Surprisingly, the convention for the b coefficients are exactly opposite (as pointed out in this answer as well), and they use a different default of sample points. Other than that, my interpretation of the documentation is that scipy.signal.freqz should be yield the same as scipy.signal.dlti.freqresp. However I can't get the same results trying it out on a small example:

import scipy
import matplotlib.pyplot as plt

filter_coef = [0.2, 0.8]
w1, h1 = scipy.signal.freqz(filter_coef)

dlti = scipy.signal.dlti(filter_coef[::-1], [1]) # accounting for flipped convention
w2, h2 = dlti.freqresp(n=len(w1))

fig, axes = plt.subplots(2, 1)
axes[0].plot(w1, np.abs(h1), label="freqz")
axes[0].plot(w2, np.abs(h2), label="dlti.freqresp")
axes[1].plot(w1, np.angle(h1), label="freqz")
axes[1].plot(w2, np.angle(h2), label="dlti.freqresp")
plt.legend()
plt.show()

Produces:

It looks like they match in terms of magnitude, but for some reason the phase is flipped.

Any ideas what is the difference here? And more generally, when should I use which version of freqz/freqresp?


texts: How does companding compress dynamic range?

texts:
Regarding when you 'compress' the signal, I typically see diagrams like this:

And the companding curve:

I can only see how the companding curve achieves a higher average SQNR by rounding down higher amplitudes more than lower amplitudes when making the quantization, according to a normal PDF, making the quantisation more efficient. I don't see how it reduces the dynamic range of the signal to the 'compressed dynamic range'. The only dynamic range compression I can see is that when the higher amplitudes are expanded again, they'll be expanded to the figure they were rounded down to during quantisation, reducing the upper dynamic range, but not the lower.
The dynamic range compression appears to be referring to the dynamic range of the signal after the compressor and before the expander, but this doesn't make sense. How can PCM samples have a dynamic range that isn't anything you want it to be.. and as the samples are uniform, surely a high rounded down amplitude now has a higher relative amplitude (sample 255/256 for instance whereas it was originally rounded down much more), and of course a lower amplitude is now a much higher amplitude.
What confused me is that there is a correspondence between 1.0 input amplitude and 1.0 output amplitude, which doesn't compress the dynamic range at all.
I hope that the key word in the output is 'relative' i.e. relative to the output range and it's actually mapping it onto a smaller dynamic range, though it looks like it's the same range. It would make sense to me if the dynamic range compression comes from the dynamic range that the samples are being mapped to is smaller than the input, and the companding curve is just a way to do this more efficiently than linear.
Why is the companding curve even called 'compression', surely the compression part is mapping it to a smaller dynamic range, or when using digital compression, the compression part is the reduction of bits in the PCM sample, and you can compress with linear companding, which wouldn't be considered compression at all if you were defining compression to be non uniform quantisation.


texts: Contrast Enhancement for subjects under Flash

texts:
For a machine vision application, we take short duration (50uSec) monochrome pictures with IR flash. The images we get are decent quality however the light color subjects (light blue, light green, yellow) are all lost or become less visible.

We convert these 10 bit images to 8 bit and look at them that way which already contributes to the problem. This is something we can probably solve with a non-linear 10-8 bit companding, so I am not worried about this part. If I look at the images on a display that supports 10 bit I still see them very very weak. My objective is to boost them.

In the raw data when I look at them, assuming surface average is 800 (pretty white), these light color subjects are 760-780 range, so the info is there but not visible. (A side note, if I take a picture of the same object under normal light with very long exposure time, 1msec, I see all details, so flash is also contributing to the problem).

What is the best way to gain these lost info back? We still have access to the 10 bit raw image data, what type of algorithm I can use to gain them back.


texts: How is Point Spread Function (PSF) related to Image Priors in Blind Deconvolution?

texts:
We are researching for our thesis about enhancements in Blind Deconvolution Image Deblurring Algorithm Applied in Image Restoration.

What really is PSF? Is it one kind of image prior? What do these both terms have to do with Blind Deconvolution?

Please enlighten us. Thanks! xo


texts: Extended Kalman filter - why not use a 2nd order Taylor expansion?

texts:
I'm implementing my first EKF and using a Taylor expansion to linearize the model. I'm somewhat new to this. Based on a paper I've read, I'm attempting to linearize the state transition matrix using a Taylor Expansion.

AFAIU, given a state vector X and state transition matrix F, we can write a linearized state transition matrix Φ (probably not the name) in the following manner:

$$
\textbf{F} = \frac{\partial \dot{\textbf{X}}}{\partial \textbf{X}}
$$
$$
\Phi = \textbf{I} + \textbf{F}\Delta t + \textbf{F}^2 \frac{\Delta t^2}{2!} + H.O.T
$$

Such that the state and covariance estimates are:

$$
\hat{\textbf{X}}_k = \Phi \hat{\textbf{X}}_{k-1}
$$

$$
\hat{\textbf{P}}_k = \Phi \hat{\textbf{P}}_{k-1} \Phi^T + \textbf{Q}
$$

Using a 2nd order expansion instead of a 1st, I've found that I get a much more accurate solution. It also allows me to use a larger timestep (minutes instead of seconds). I worry though that in doing so, I'm not actually linearizing, since the 2nd order TSE corresponds to a quadratic approximation.

Assuming the above steps are correct, is there anything wrong with taking the 2nd order TSE?


texts: Why do we say that "zero-padding doesn't really increase frequency resolution"

texts:
Here is a sinusoid of frequency f = 236.4 Hz (it is 10 milliseconds long; it has N=441 points at sampling rate fs=44100Hz) and its DFT, without zero-padding :

The only conclusion we can give by looking at the DFT is: "The frequency is approximatively 200Hz".

Here is the signal and its DFT, with a large zero-padding :

Now we can give a much more precise conclusion : "By looking carefully at the maximum of the spectrum, I can estimate the frequency 236Hz" (I zoomed and found the maximum is near 236).

My question is : why do we say that "zero-padding doesn't increase resolution" ? (I have seen this sentence very often, then they say "it only adds interpolation")

=> With my example, zero-padding helped me to find the right frequency with a more precise resolution !


texts: Is there better than hamming code for correcting errors

texts:
As known, hamming code can be represented by (7,4) where data are 4 bits, and 3 bits for parity,  and (15,11) where 11 bits are data and 4 bits for parity. It can correct up to one bit, or detect two bits without correcting.

Parity bits will cause an issue of spectral efficiency or in data rate, so my question, is there another code better than hamming code which can correct more errors with less parity bits ?


texts: how to get the phasor(frequency domain) representation of a digitalized signal?

texts: I know that we can get the phasor representation(frequency domain) of a signal through Fourier transform, but the picture above gives a method that can get the phasor representation of a signal in a different way, in this method according to the writer, it numerically "demodulate" the signal to create a phasor representation of the signal, however, i can't understand that, anyone can help me? Thanks in advance!


texts: Statistic that is responsive to changes in time series, yet not too volatile

texts:
Let's say I have a fairly volatile time series $X_t$ - it doesn't have any reason to show an upward / downward trend, but it does show drops and spikes from time to time. It can also change level (e.g. for a few weeks it will be at one level, but it will quickly drop down to a new level a few % points down and stay there).
I'd like to have a statistic (basically a time series Y_t) that signals a change in $X_t$ quickly. It should also converge to the same (or close to the same) level as $X_t$ with as little delay as possible. Finally, $Y_t$ should have as little volatility as possible (I'm aware there's a tradeoff between low volatility and responsiveness).
I've tried the rolling mean - if there's a short duration (say 1 period) drop in $X_t$, then the rolling mean reacts quickly to the drop, but shows a delay in reacting to the recovery in $X_t$ (i.e. when $X_t$ is restored to normal values). The delay is the same as the window size minus the drop duration.
At this stage, I'm considering two statistics. The exponential moving average with a window of $5$:
$$E_t=\frac{\sum_{i=0}^t(1-\alpha)^iX_{t-i}}{\sum_{i=0}^t(1-\alpha)^i}$$
where $$\alpha=\frac{2}{\text{window}+1}$$
The other statistic I'm considering is a combination of EMA and median. Let $M_t$ denote the rolling median with a window of $9$. Then
$$C_t=aE_t+(1-a)M_t$$
To test their response to a drop, I did a quick simulation as follows:
a = pd.Series(np.hstack((np.ones(10)*0.81, np.ones(4)*0.77,
              np.ones(4)*0.72, np.ones(5)*0.77, np.ones(10)*0.81)))
x = 0.5
b = a.ewm(span=5).mean()
c = x*a.ewm(span=4).mean() + (1-x)*a.rolling(7).median()

figure(num=None, figsize=(27, 9))
plt.plot(a, marker="o")
plt.plot(b, marker=".")
plt.plot(c, marker=".")

By the time $C_t$ converges to the same values as $X_t$, the difference between $C_t$ and $E_t$ is really small (less than a percentage point).
$C_t$ has the feature that it won't drop as sharply as $E_t$ in response to a drop in $X_t$, but after a few periods it will abruptly drop down to the level of $X_t$. Similarly when $X_t$ recovers, $C_t$ initially shows a slower rise than $E_t$, but then suddenly shoots up, overtaking $E_t$ and catching up with $X_t$.
Now intuitively it seems that overall the EMA $E_t$ would be a better statistic, but I'm not sure. Is there any concrete, rigorous justification for the superiority of any one statistic over the other? For example, is it possible to show that maybe $C_t$ reacts similarly to two different scenarios, and is therefore unsuitable to differentiate between those scenarios?
Finally, given the requirements of quick response, quick convergence and low volatility, is there any other statistic I may have missed that's better than these two?


texts:
How can I know exactly what frequencies are not influenced by the filter?

I mean, I have a signal in the range of frequencies $5-100 Hz$ and I apply a high-pass, $4^{th}$ order Butterworth filter with cut-off frequency of $8 Hz$. The filter will not cut sharply at $8 Hz$, instead it will have a slope related to the filter order: the higher the order, the steeper the slope right?.

So, how can I compute the frequency from which the signal is left unchanged (so no attenuation due to the filter)?

Is there a way to obatin this info directly in matlab?


texts: Python: How to calculate image perimeter orientation histogram?

texts:
I would like to find the perimeter orientation histogram of an object inside an image.
The object is returned into an array but I don`t know how to find the orientation of its perimeter.
I think the first step to find the orientation histogram of the perimeter is to find the perimeter itself. That is to return the perimeter pixels coordinates as a 2-D array of X and Y.
What I found in Python is to just return the total number of pixels in the object perimeter but no information about the coordinates of perimeter.

Can anyone help me find the perimeter coordinates of an object and then calculate the orientation of the perimeter pixels?

More information about perimeter orientation histogram can be found in this paper in section 3.3 Feature Extraction.


texts: Estimation of the amplitude of a frequency component in a signal

texts:
I'm working on a device that produces a signal with sampling rate of 100 kHz. Every 40 ms, I receive 16000 points, from which I have to extract a frequency close to 3.55 kHz. I have developed a method that extracts the frequency of the signal very close to the Cramér-Rao lower bound, but I need to extract the amplitude of the signal at that frequency very precisely.

The first idea that occured to me was using a digital implementation of a lock-in amplifier, but the problem is that a lock-in amplifier needs a low-pass filter. And I'm trying to avoid using filters as much as I could, because filters introduce a systematic shift in my readings, while I need extremely high precission and accuracy.

What's the best amplitude estimation algorithm I should use?


texts:
I am very new to Fourier analysis, but I understand that through the use of the Fourier transform a signal in the time domain is displayed in the frequency domain, where frequency values are normally displayed along the x-axis, and amplitude is displayed along the y-axis.  However, at one point in the textbook I am using, the following is stated:

Let us assume that we have the function $f(t) = \cos(\omega_0 t)$.  The spectrum then consists of two delta-functions

$$F(\omega) = \pi \delta(\omega - \omega_0) + \pi \delta(\omega + \omega_0)$$

__

This confuses me.  When we have $f(t) = \cos(\omega_0 t)$, then I would assume that the Fourier transform should yield an amplitude of $1$ at $\omega = \omega_0$ and $0$ elsewhere. But the delta function is defined as:

$$\delta(\omega - \omega_0) = \left\{ \begin{array}{1 1} \infty & \quad \omega = \omega_0 \\ 0 & \quad \omega \neq \omega_0 \end{array} \right.$$

So wouldn't this give an infinite value at $\omega = \omega_0$?

If anyone can explain the intuition behind the statement in my textbook, then I would be very grateful!


texts: Why is there a mismatch between normalized cross-correlation and perception for basic waveforms?

texts: To my ears, sawtooth and square waves are more similar to each other than either one is to a triangle wave. However, when I compare their magnitude spectra against each other using normalized cross-correlation, the square and triangle waves are the most similar to each other, followed by the sawtooth and triangle waves. Does anyone know what the source of this apparent mismatch is?


texts: How is "Joint class histogram" algorithm deploy using Matlab?

texts:
Image registration and image retrieval depend on image similarity measure to find mutual information between 2 images. Using histogram, joint histogram, and joint class histogram we can find the degree of similarity between those images. Hopefully, I understand and deploy both first techniques.Here, I found joint histogram algorithm and deploy using Matlab.

For example:
this paper was mention this histogram and used as image similarity measurement, but I couldn't found how the authors calculate it.

However, I search for many days to find good paper, book or article regarding joint class histogram but their authors only mention this technique without enough detail. I think this technique still infant and doesn't have enough study yet.

I think joint class histogram algorithm will the best answer for my question.


texts: How to predict the cramped frequency of a digital filter based on an analogue frequency?

texts: I have a second order analogue high pass transfer function (unity gain at infinity). It's magnitude response hits the -20 decibel line at a frequency of 5706 Hz (the corner frequency is 18000 and sample rate is 44100). When I convert this analogue filter to a digital IIR via the BLT method, the digital filter's magnitude response hits the -20 line at 11149 Hz. This will happen with the standard RBJ audio EQ cookbook formula i.e. it's not a problem with my maths. I know that digital filters suffer from cramping and that's what I'm interested in. Is there a formula to predict what an analogue frequency will be cramped to when the digital filter is realized given the analogue filter's frequency? i.e. something that will tell me 11149 given 5706. Just to be super clear I'm not asking how to do pre-warping using: $ tan(\frac{w_c}{2}) $.


texts: Is the sum of all FFT bins from Audacity equal to RMS amplitude?

texts:
Say I open an audio segment with Audacity.  Then I perform a size 4096+ FFT with a rectangular window.  The program outputs bin magnitude values in dB.  If I add them all together like this:

sum = 10*log10(10(bin1/10) + 10(bin2/10) + 10(bin3/10)...).

I get a value very close to the file's calculated "total RMS amplitude".  (Or in this case, 3dB off since Audacity considers a full-scale square, instead of sine, to be 0dBFS).

Is this expected, a coincidence, or just an approximation?


texts: What's the advantage of adaptive IIR filter against FIR?

texts:
Adaptive IIR filters is not straightforward, and may be unstable. Many people say that adaptive IIR filters use less coefficients than FIR filters. What I'm curious about is how many coefficients can IIR save?

I tried to use adaptive IIR filters to estimate transfer function of a 32-order FIR filter. Assume the IIR filter has $M+N+1$ coefficients: $a_1, a_2, ..., a_M, b_0, b_1, ...b_N$. I found the estimation result is acceptable only when $M+N+1 \ge 30$, i.e. only 2 coefficients can be saved.

In actual projects, for example, a 50 MHz FPGA, a 32-order FIR will produce about $(32 / 50 ~{M}) / 2 = 0.32 ~{\mu s}$ delay, so

What will happen for IIR?
Can adaptive IIR filters really reduce number of coefficients and reduce signal processing time delay?


texts: find a line go through the center mass and that line must be inside object area

texts:
I have an object (orange color object) and I know its centroid/center mass C. I want program to automatically find another point ( can be point 1, point 2, point 3 or point 4 and so on) such that the line connect point C and one of those points (so as green line in figure 2) must stay inside the object. Moreover, the green line length must be as long as possible, or at least 3/4 of the maximum length. The yellow point in figure bottom is not qualified because the violet line which connect it to the center mass C has a part which is not inside the object.

Any idea? I am doing it on Matlab (image processing)


texts: Symmetry and periodicity of frequency-shifted discrete Fourier transform

texts:
The unitary discrete Fourier transform (DFT) of a sequence of numbers $x_n$ to $X_k,$ with integer $0 \le n < N$ and $0 \le k < N,$ can be defined as:
$$X_k = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}\tag{1}$$
and the inverse discrete Fourier transform (IDFT) as:
$$x_n = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k e^{2\pi ikn/N}\tag{2}$$
If $x_n$ is modulated (multiplied) by a unit-magnitude zero-phase complex sinusoid $e^{-2\pi ibn/N}$ before DFT, and the IDFT output is demodulated (divided) by the same, then we get another transform pair from the family of generalized discrete Fourier transforms, parameterized by the constant $b:$
$$
\begin{aligned}
X_k(b) &= \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_ne^{-2\pi ikn/N}e^{-2\pi ibn/N}\\
&= \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_ne^{-2\pi i(k+b)n/N}
\end{aligned}
\tag{3}
$$
$$
\begin{aligned}
x_n &= e^{2\pi ibn/N}\frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k(b)\cdot e^{2\pi ikn/N} \\
&= \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k(b)\cdot e^{2\pi i(k+b)n/N}
\end{aligned}
\tag{4}
$$
Whereas DFT samples frequencies $2\pi k/N,$ the frequency-shifted transform samples frequencies $2\pi (k+b)/N.$ This can be visualized on the Z-plane:
Figure 1. Z-plane representation, showing the unit circle, of the frequencies sampled by A) DFT and B) the frequency-shifted DFT, with $N = 12$ and $b = 1/3.$ The shift $2\pi b/N$ appears as the angle between the real axis and the vector from origin to one of the sampled frequencies.
Question: What are the time-domain symmetry and periodicity properties of such frequency-shifted Fourier transforms when extending $n$ beyond $0\le n<N$ in the inverse transform, and how does this depend on the parameter $b?$
For DFT (or with $b = 0$) the time-domain extension has period $N$ with no time-domain symmetry imposed by the transform.


texts: About the orthogonality of columns in Zadoff-Chu matrix

texts:
As we know that Zadoff-Chu matrix is similar to Walsh-Hadamard matrix where every columns in those matrices is orthogonal with the any other column.

For Walsh-Hadmard matrix, the orthogonality is clear and can be demonstrated easily. But regarding the Zadoff-Chu matrix, I don't see that orthogonality is right as shown below:

a = (sqrt(2) + 1j*sqrt(2))/2;
W = [a 1 a -1; -1 a 1 a ; a -1 a 1; 1 a -1 a];  %%Zadoff-Chu matrix

C1 = W(:,1);     %% Take the first column

orth_results = [];
for j = 1 : size(W,2)
    results = (abs(sum(W(:,j).*C1)).^2);  %Project the first column to every column in matrix W
    orth_results = [orth_results results];
end

Following the condition of orthogonality between columns, orth_results should be equals to 0 except in the fist value, but what I get is orth_results = [8.0000         0    8.0000      0];

What's the problem of that? Is there any issue in the matrix itself ?


texts: DDS - Supersampling above the FPGA clock frequency

texts:
I'm using a Xilinx Virtex 6 FPGA with 200 MHz clock connected to a 1.6 GSa/s DAC to generate sine waves of up to 800 MHz using 8 parallel Xilinx DDS Compiler IP cores.

Initially, I calculate the 32 bit wide phase increment using the formula provided in the manual with a 32-bit wide phase accumulator, 200MHz clock frequency and desired output frequency.

Then 8 DDS modules are instantiated each with a phase offset of $0,\frac 18 \Delta_\phi,\frac28 \Delta_\phi,\ldots, \frac78 \Delta_\phi$. These output 8 samples of a sine wave on every 200 MHz rising clock edge, matching the higher sample rate of the DAC.

The DDS produces correctly sampled output waves up to a frequency of 200MHz. However for f > 200 MHz, the output waves are the wrong frequency. For example, if I input $f =300\,\text{MHz}$, the generated wave is 100 MHz; if I input $f = 400\,\text{MHz}$, the generated wave is 200 MHz.

My suspicions are there is an overflow at some point in the calculation of the phase increment and phase offset values due to their 32-bit width. My question is how would I get around this to correctly produce waves of up to 800 MHz?

I am not concerned with effects of undersampling the higher frequency waves, i just want to see the correct amount of samples per cycle i.e. 2 in the case of $f =800\,\text{MHz}$.


texts: Why do we use the term "cepstral" in MFCCs, and how does it related to the cepstral domain?

texts:
As I know, we can get the cepstrum by taking the Inverse Discrete Fourier Transform (IDFT) of the log of the power spectrum.

And we can get the Mel-frequency cepstral coefficients (MFCCs) by taking the Discrete Cosine Transform (DCT) of the log filterbank energies.

But I didn't understand the following two points:

Why do we use the term "cepstral" in MFCCs when we use DCT instead of IDFT?
How does MFCCs relate to the cepstral domain according to the reference here?


texts: How to calibrate three orthogonal cameras？

texts:
I am trying to find out intrinsic parameters(focal length ,principle point,and lens distortion coefficients)as well as the extrinsic parameters(rotation and translation of each camera relative to the reference frame) of three cameras in the flowing picture.

They are placed at nearly orthogonal angles around a 12x12x12 inch working space. There are four steps:

capture images of chess board by two cameras at the same time
use toolbox calib_gui(matlab) to calibrate two cameras separately and get Calib_Results_left.mat and Calib_Results_right.mat.
put these two results into stereo_calib(matlab) and run Calibration
repeat step 1~3 on another pair of cameras

However, I was stuck at step 3! None of pair images was found inconsistent. The errors were shown in the following form:

Disabling view 2 - Reason: the left and right images are found inconsistent (try help calib_stereo for more information)

I checked the order of images but didn't find problems. I doubt that my application is not binocular calibration, because intersection angle between two cameras is too big
for a stereo-camera!

I just want to find out rotation and translation of each camera relative to the reference frame. Does anybody have any ideas?


texts: How would sigma-delta modulation work in software?

texts:
I've recently come across this YouTube video that shows a PCM stream being played back via the 1-bit internal "PC Speaker", and after hours of googling I've come to think this is done via Sigma-Delta modulation.

So far, what I gather from the process is that it takes the original digital signal, oversamples and approximates it via PWM while applying a low-pass filter to the 1-bit output. I'm not quite sure if this is right.

Not having a background in DSP, I'm wondering:
How could this process be implemented in software?


texts: OFDM sensitivity to sampling time offset (with half of time interval between samples)

texts:
In simple QAM modulation symbol synchronization is crucial: the samples need to be taken (or interpolated and resampled to ) exactly corresponding to the time when the symbol was emitted.

This is why symbol syncronization methods like Muller-Muller was developed.  I call this - for the sake of contrasting it with coarse OFDM symbol syncronization (achievable without interpolation and resampling) - fine-symbol-syncronization or smaller-than-time-interval-between-samples syncronizaiton.

Now, the question is, if such fine syncronization is important or not for OFDM. In OFDM there is the cyclic prefix which is much longer than the time interval between samples and there is also the FFT. It is not clear if in the case of OFDM the fine-symbol-syncronization is really neccessary. It might be that the guard interval+FFT makes it unneccessary.

One reason why fine-syncronization should not matter is that the purpose of the guard interval (cyclic prefix) is to make demodulation insensitive to channel response (multi path due to reflections).

Is it important with OFDM to carry out fine, i.e. smaller-than-time-interval-between-samples syncronization ?

EDIT: this book seems to discuss this question in Section "5.2.4 Synchronization Errors", it says that timing offset causes rotation - increasing in strength for sub-carriers further away from the center - need to read it more carefully.


texts: Convolving with the frequency response of a filter

texts:
With scipy I can use signal.butter and signal.lfilter on a time signal.
Performing a Fourier transform then shows me that the filter worked.

Using signal.freqz I can get the frequency response of the filter. Is there anything stopping me from simply taking the FT of my time signal and convolving it with my frequency response?


texts: Fourier transform is an isomorphism...but we don’t get when each frequency appears?

texts:
Statistician here who wants to get some DSP knowledge for time series analysis.

I’ve known for years that if we hit a function with a Fourier transform, we have an inverse Fourier transform that will recover the original function. However, doesn’t the interpretation of the Fourier transform in the frequency domain lack the time component? In other words, we can say that 100Hz appears in the signal with some intensity, but we can’t say if that appears at the beginning of the signal or the end. Those are very different signals to me, yet I’m supposed to be able to recover each by applying the inverse Fourier transform?

There seems to be an inconsistency: we either lose the time information and can’t invert, or we retain the time information and can invert.

What is the resolution to this apparent inconsistency?

(I have a hunch that it has to do with the imaginary part of the Fourier transform, though I’m not sure how.)


texts:
I have read that non-uniform quantization boosts the smaller amplitude signals by a large amount. However the larger amplitude signals receive a small gain. As shown in the below diagram (Compressor Input and Output).

The Input - Output Characteristic diagram of compressor is given as below.

My query is as seen from the Input - Characteristic diagram, the step size is very less for smaller amplitude signal and hence the gain has to be less. On the other hand the step size is more for the large amplitude signal. Hence the gain has to be more. But the compressor input and output are in reverse way. Can you please help where I was going wrong?


texts: Image processing free open source libraries

texts:
This is the first time that I post a question  here so excuse me for asking such a beginner question.
what free open source 3rd party libraries which have good documentation(because I'm new to using 3rd party libraries) would you suggest me for fourier transform,inverse fourier transform , meshgrids , dynamic arrays and etc to have a better performance and higher speed in visualC++ when implementing image processing filters?
Note that I'm not permitted to use openCV. In fact I'm going to compare the results of my code with results of openCV


texts: What are the advantages of using fft2 over fft

texts:
Assume we have a matrix x of size (8,8), . As known,  FFT(x) performs 1D-FFT transformation, column wise. However,  FFT2(x), performs 2D-FFT transformation.

In that case, what's the advantage of using 2D-FFT over 1D-FFT ?


texts:
I'm currently working on a piece of DSP/Synthesis software and I have run into the issue of generating band limited waveforms. I have been doing some research to try and find a good solution but I haven't had much luck. I have found some things about a BLIT to generate complex waveforms but I have not found a good example of an implementation. I am not very good with Math or Mathematical notation but I am great with dealing with c/c++ code. If anyone has a good explanation of how to generate a BLIT using c/c++ or resources about that I would much appreciate any help that you could give.

Thanks

P.S. I know this could also be a Stack Overflow question because of the programming nature but I thought I would try the DSP community first. Thanks again.


texts:
double** table_new_pulse(ULONG length, ULONG nharms, double dcyc)
{
    ...

    ULONG harm = 1;
    double step = TWOPI / length;
    double dcyclerad = M_PI * dcyc;
    for (ULONG i = 0; i < nharms; i++)
    {
        double amp = sin(harm * dcyclerad) / harm;
        for (int j = 0; j < length; j++)
        {
            table[j] += amp * cos(step * harm * j);
        }
        harm += 1;
    }

    // normalize
    ...

    return table;
}

The preceding C function tries to approximate a pulse wave of duty-cycle dcyc by summing cosines of the first nharms each with an amplitude of sin(sin(harm * dcyc) / harm) and storing the values in a lookup table. If length is 1024 and nharms is 8, when dcyc is 50% the resulting waveform  has equal positive and negative peak amplitudes:

But as dcyc decreases from 50%, the negative peak amplitude also decreases. If dcyc is 33.3%, the negative peak amplitude is -0.540:

If dcyc is 10%, the negative peak amplitude decreases even further to -0.167, rendering the resulting waveform almost unipolar:

Can anyone explain why this happens in the simplest way possible? How can one generate a bandlimited pulse wave with varying duty-cycle and with equal positive and negative peak amplitudes?


texts:
The context is a supervised learning problem. I have a dataset with features X and a target y. Among that dataset there is the date of the day encoded as an integer starting at 1. I have around 10k instances by days. I have found out that encoding the time on a circle (creating $cos(2*\pi t / T), sin(2*\pi t / T))$ for  $T = 7$ or $30$ as supplementary variables work pretty well to predict the outcome y.
Now I have two problems : in the test set the time also starts at 1 , so there will be some unknown phase and the test set will come in a streaming way, that is instance by instance.
Are there any streaming algorithm that would help me to rapidly get the phase for a given period and allow me to realign the circle encoding ?


texts: How to represent the nonlinear model as a state space in Unscented Kalman Filter

texts:
There is an Autoregressive model of order 1 (AR(1)) that is excited by a non-linear signal as the input:
$$x_t = \rho x_{t-1} + u_t \tag{1}$$
The time series $u_t$ is generated from a nonlinear map,
$$u_t = f(u_{t-1},\mathbf{w}) \tag{2}$$ where $f$ is the nonlinear function.
The observations are $$y_t = x_t + v_t \tag{3}$$ where $v_t$ is the measurement noise that is an Additive White Gaussian Noise.

Q1: Can I re-write the model (1) as a state space in the following way:

$$x_t = Ax_{t-1} + f(u_{t-1},\mathbf{w}) $$

$$y_t =  x_t+ v_t \tag{4}$$

Is the above representation correct? If not then I shall be grateful for the correct technique to represent it.


texts:
I am trying to understand the SSIM index to compare two images for similarity, but I am having difficulties understanding how it is done.

Therefore, I would like to ask, if someone can explain it step by step, possibly also using some images to visualize the process ?

That would be of great help !

What I have undertstood so far:

I take two pictures
I split picture 1 in N squared subpictures
I split picture 2 in N squared subpictures at the same location as in picture 1
I compute the SSIM index for each local subpicture


texts:
I'm programming in c# using opencv (emguCV). I use motion detection approach to detect some changes in a paper wall. For example it's may be colored by a laser pointer. My problem is that when I use motion detection, some unwanted changes are occurred by wind, sun and etc. Also I want to save the time of event so I can't use thresholding methods because they find point of laser totally, not their times.
Is there any approach to filter mentioned effects?


texts: Digital Signal Quantized Processor Complexity

texts: In “traditional” digital signal processing, the complexity is computed as a number of multiplication the operation requires, e.g. the computation of the N-point DFT via the decimation-in-frequency FFT requires (N/2)log_2(N) complex multiplications (and the same number of additions, but in the DSP the additions are assumed to require less time and their contribution is traditionally omitted, unless specifically required). However, we typically do not define which numbers we multiply. My question is how we define the complexity in the DSP sense, if we apply some quantization rule for the numbers to be multiplied: e.g. we proceed from multiplying two real numbers in all operations of the processing devise to multiplying the pairs of int8, float 16, etc. The latter operations are evidently “simpler” that the multiplication of two reals, but how can we quantify and express this complexity reduction for the DSP element?


texts: How do I perform a spectral averaging? Amplitude spectral density

texts:
I haven't been doing much signal processing over the years and I'm wondering if my logic is sound. I'd appreciate any input.

So, I have a data set with independent variable ($A$) in Frequency ($1\textrm{ kHz}$-$101\textrm{ kHz}$) and the dependent variable ($B$) with units $\textrm{ Pa}$ (Pascals). These are experimental data.

I currently have a graph of $B$ ($\textrm{ Pa}$) versus $A$($\textrm{ Hz}$).
However, I'm required to represent my data in graph of $\frac{\textrm{ Pa}}{\sqrt{\textrm{Hz}}}$ versus $\textrm{ Hz}$. ($C$ being the new variable with units of $\frac{\textrm{ Pa}}{\sqrt{\textrm{Hz}}}$, Amplitude Spectral Density, or the square root of the Power Spectral Density)

So my question is:

Would it be grossly inaccurate to just average my dependent variable $A$ (with Units of $\textrm{ Pa}$) over the square-root bandwidth ($100\textrm{ kHz}$) of data $B$ to derive $C$?
Or should I simply (not so simple,but) ifft(A). Then use PSD = [fft(x).*conj(fft(x))], and take the square root of PSD?


texts: Which Noise Reduction Algorithms Are Used in Commercial RAW Image Processors?

texts:
I'm trying to guess what noise reduction algorithms are used in commercial processors for raw images from digital cameras. I find this fairly easy to do for the sharpening algorithms (most use unsharp mask as at least one option, for example, and this is clear from the behavior and adjustable parameters), but I know next to nothing about noise reduction. I would like to

identify as well as possible the basic algorithms being used (of course there will be some proprietary optimizations), and
identify the corresponding academic papers that derive them, to learn their inner workings and use as a springboard to the contemporary literature (via tracking citations).

The parameters for the noise reduction algorithms in two common processors, Capture One and Lightroom, are shown below. Here are some further clues.

They're fast. Each works in about a second or less on my modest laptop on 24 MP raw images.
They're old. Lightroom was released in 2007, and I don't think the core options have changed since then.
They operate separately on luminance and color noise.
I don't think they're wavelet-based, since (as far as I know) these applications do not use wavelets anywhere else.

Question 1: Is there a well-known, industry-standard algorithm that meets this description?

Among all raw processors available, the best noise reduction seems to come from DxO's Photolab. (This is of course debatable, but it consistently places at or near the top of comparisons I've seen.) I'd like to try to understand how it operates. A screenshot of its parameters is attached.

DxO describes it as a local algorithm that works pixel-by-pixel:

  “the denoising algorithms (in PRIME) analyze the structure of Raw
  images in depth: more than a thousand neighboring pixels are surveyed
  for each pixel. This extensive exploration identifies similar data
  (for) use (in reconstructing) image information.”

Besides that, we know a few more things:

It's relatively new (released 2014).
It's slow. Users report it takes a minute or more per photo.
It's recommended for very noisy, high ISO photos. DxO and users report no advantage over faster noise reduction techniques for photos with mild noise.

Question 2: What algorithm could DxO plausibly be using?

I realize this is a somewhat unusual question, so thanks in advance for letting me access your expertise.


texts: area under the phase encoding gradient in MRI

texts:
I have a question about how I could compute the area (duration and amplitude) for the phase encoding gradient in a typical MRI gradient echo sequence. So, let us assume that my field of view in the phase-encode direction is 20 cms and I have a k-space matrix of size 256. From the tutorial here, I can see that the k-space spacing is $1/20e-3$ per meter. Let us also assume that the maximum gradient strength I can achieve is 10mT/m.

I am now stuck and not sure how I can compute the area that I need under the phase encode area. I am aware that the area will change at each phase encoding step and at least from the pulse sequence diagrams that I see for a GE sequence, it seems that the amplitude is adjusted (going for negative max to positive max) and the gradient time is kept constant during each phase encode step.

However, I am at a loss as to how to compute this area?


texts: Show that $S[n]= \frac1M \sum e^{\frac{2\pi l n}M}$ for $l = 0$ to $M - 1$

texts:
I am learning about image processing. Currently, I am reading about upsampling and down sampling.
Here is some background what I am reading:

I was wondering if someone show me a proof of how equation (70) is true. It simply says "It can be shown that".

Thanks!


texts: Variation of Pulse width with center frequency in Pulse Compression technique

texts:
Assume a chirp with BW = 1 Hz (frequency range 0-1 Hz) for a period of 700 sec. The pulse compression technique (obtained through auto correlation function) applied on the chirp  generates a pulse with the main lobe width of 1.0 sec.
However, on changing the frequency range to 1-2 Hz ( the BW is same) , the width of the main lobe of the generated pulse changes to 0.6. The image shows

I would like to know the reason behind this. Theoretically, the width of the main pulse is dependent only on the BW. (=1/2BW)


texts: Signal Processing using Fourier Transform

texts:
So I'm trying to understand how MRI machines work. I understand all the concepts of it, the parts, what they do, how the machine works, etc. The part I'm having trouble with is the fourier transform business.

In the MRI machine, when the hydrogen atoms release the radio frequencies previously absorbed, the transceiver absorbs it and uses fourier transform to convert those signals from the time domain to the frequency domain.

The part I don't understand is, how does the fourier transform end up producing the image of the body. How do you go from math to a detailed picture of the brain??? Can anyone explain that?


texts: Masking the mains hum to improve anonymity

texts:
Today, a news blurb appeared that claimed that attackers could identify the location in which an audio record was taken by matching the mains hum to known data. Sadly, I fear that this might be feasible. The mains hum has quite a tight tolerance, but any aberrations (that cannot be done away with) could be spotted easily. In fact, power suppliers have to observe the hum and adjust their stations accordingly.

How could you manipulate a record such that this information is destroyed completely, so that it is not reconstructible by any means? Filtering alone might leave enough remains to still perform that identification. After all, all filters have a finite number of poles.


texts: Applying a Power on a Sum of Gaussians to Enhance Resolution

texts:
There is so-called Power Transform technique in signal processing where the unresolved signal is sharpened by raising each data to a constant positive power. The example shown below is a sum of two Gaussians, $G_1$ and $G_2$ plotted against time. In order to resolve the peaks one can apply a power of 5 to each data point. The result is the resolved pair as shown in the bottom figure.

The above operating is equivalent to raising the sum of Gaussians $G_1$ and $G_2$ to the power 5, i.e., $(G_1 + G_2)^5$ then according to the binomial theorem, we should get 6 peaks, because not only we have the first and the second peak raised to power 5, we will have a "mixture" of the two peaks as well- the other terms of the binomial theorem. No matter how well we zoom, we cannot see the missing $G_1^4$$G_2$, $G_1^3$$G_2^2$ peaks? If power transform of 5 is applied to two peaks, we are left with two peaks not six.

Where do the peaks corresponding to the other binomial terms go?

Thanks.