texts: Goertzel filter Poles on unity circle influence of a behavior of a signal at high frequency

texts:
I'm trying to apply a goertzel filter on my sound data of 20500 Hz with sampling rate of 44100[samples/sec]. My bit duration which is passed to the goertzel function is returning a very low magnitude at high freqs. So, I've noticed that as bit duration is lower the magnitude also is lower. The 2 parameters of increasing the freq and decreasing the singal's duration is making the magnitude be closer to 0. For instance: 0.00005.

I saw on wikipedia that the filter has 2 poles on the unity circle of it's Z-Transform - causing it be unstable.

My question is: if I use a very small signal duration, for example 5ms would it make a sense to duplicate this frame's singal several times to make the whole new frame larger so the goertzel have a lot of more samples?


texts:
i have 2 questions:
1-If zeros on the unite circle that will be min or max phase?and this system has low-pass frequency response or not ?
2-Is the min phase has the shortest impulse response?


texts: Pre-Warping is not working on bilinear transform

texts:
New here and I just wanna ask the question of what happens when your pre warping doesnt work on your bilinear transform?
I am assuming my isnt working as part of my research I have read that the higher the frequency is the more non linear it becomes, so I am assuming that my frequency I chose just has to many non linearities that it wouldn't ever get re-aligned to the continuous function critical frequency. It tries  but doesnt get there.
Sampling Frequency = 44.410kHz
Nyquist Frequency = 22.205kHz
Critical Frequency = 20kHz
Ratio = Fs/Nyquist frequency =  20/22.205 = 90.1%
I noticed the closer you are critical frequency is to your Nyquist frequency when applying the bilinear transform it gets ugly. And it looks like the pre warping can only do so much.
The only way I see around this if I just changed my Sampling frequency to like 96kHz which would make my Nyquist frequency = 48kHz would be then the ratio = 42% this would ease the constraint of the non-linearities and have the pre warping affect work much better, but this raises the question is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz? Because the way this is going I dont see it.

CODE:
s = tf('s');

Fs = 96e3;
Ts = 1/Fs;
Fc = 20e3;

GP = -3;
GS = -20;
WP = 2*pi*Fc;
WS = 251327.412;

n = ceil(log10((10^(-GS/10)-1)/(10^(-GP/10)-1))/(2*log10(WS/WP)));

WC_1 = WP/((10^(-GP/10)-1)^(1/(2*n)));
WC_2 = WS/((10^(-GS/10)-1)^(1/(2*n)));

LP = 1/((s/WC_1)^4+2.613*(s/WC_1)^3+3.414*(s/WC_1)^2+2.613*(s/WC_1)+1);
LP_2 = 1/((s/WC_2)^4+2.613*(s/WC_2)^3+3.414*(s/WC_2)^2+2.613*(s/WC_1)+1);

LP_44_410kHz = c2d(LP,1/44.410e3,['Method','tustin','PrewarpFrequency',WP]);
LP_96kHz = c2d(LP,Ts,['Method','tustin','PrewarpFrequency',WP]);

options = bodeoptions;
options.FreqUnits = 'Hz';
options.Title.String = 'With Pre-Warping';
bode(LP_44_410kHz,LP,LP_96kHz,options);


texts: Confused about filter's cutoff frequency in rational sampling rate converter

texts:
A interpolator which increase sampling rate by L times is followed by an anti-imaging filter, the filter is said to have a cutoff frequency of Fs/2L. Is Fs here is the frequency before or after interpolation?

Also during rational sampling rate conversion, the filter in between the interpolator and decimator is said to have cutoff at min(Fs/2L, Fs/2M). But clearly Fs isn't same for both of them. The interpolator changes the sampling rate before the signal passes through the decimator.


texts: Why is there a time-frequency precision tradeoff

texts:
I understand the science - I think. Looking at a spectrogram, the lower frequencies are fuzzier due to a lower ratio of: frequency/time.
My main question is: why can I as a human hear lower frequencies sharply, but the FT cant?
Is this a psychological thing - perhaps my ears suffer from the same issue, but my brain is filling in the missing information?

Notice how the band gets wider, demonstrating the increase in innacuracy.


texts:
I have a signal which looks something like:

The 12 main round peaks is what I'd like to keep. However, the peak detection algorithm I apply is picking up the in-between noisy peaks. Any strategies I could use to get rid of this noise?

I tried smoothing the signal, but the noisy peaks are still there.

Thanks


texts: Is there an intuitive explanation for what complex basebanding a signal is?

texts:
I have an algorithm which takes an input signal, $x[t]$, and (among other things) complex basebands the data and outputs the phase of the complex numbers.
I'm new to the field and the Matlab function which produces the basebanded signal is a kind of black box to me. I've tried reading the code, but it's pretty impenetrable and not commented. I also tried to input different shapes of signals, but couldn't learn much from it. My questions are:

What exactly is the complex baseband of a signal? I sort of understand what "baseband" means, but where does the complex phase come from?

Related, does the phase have any relation to the signal? If I input a sinusoid, what should I expect the phase output to be?


texts: How to recognize hexagonal tiling in boardgame?

texts:
I would like to recognize the boundaries of a hexagonal tiling in a photograph, like in  image below:

It seems to me, that a standard approach at a square grid is to first detect corners (e.g. canny) and then extract the longest lines via a Hough transform or something similar.

This does not look as optimal solution with hex tiling, because the length of outer lines is shorter and it's hard to separate them from other lines.

Is there any algorithm to adress this problem? It would be particulary nice to have a solution in opencv, but I am also interested in general ideas.

update:

With python and opencv I was able to receive this result:

Here is my code:

import cv2
import numpy as np

imgOrig = "test1";
img = cv2.imread(imgOrig+".jpg");
lap = cv2.Laplacian(img, cv2.IPL_DEPTH_32F, ksize = 3)
imgray = cv2.cvtColor(lap,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(imgray,127,255,0)
contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)
size = img.shape
m = np.zeros(size, dtype=np.uint8)
for i, cnt in enumerate(contours):
    if cv2.contourArea(cnt) >= 1:
        color = (255,255,255)
        cv2.drawContours(m, cnt, -1, color, -1)
cv2.imwrite(str(imgOrig)+"contours.jpg", m);

Laplacian of image looks like:

I will try to optimate the parameters of this approach and then try to interpolate the boundaries of the four sections.


texts: Why level of noise can be magnified twice through each numerical differentiation?

texts: I was reading a paper and saw this is mentioned there, but I cannot figure out how this can analytically be proven?


texts: FFT based symbol synchronization for digital demodulation

texts:
Where can I read more about "slow search" methods for symbol synchronization based on FFT for extracting a clock signal from the modulated signal?

I have read here (8.7 Symbol synchronization)

  One method of synchronization is to extract a harmonic of the symbol
  frequency from the received signal. Then a local symbol clock can be
  synchronized by methods that are very similar to the phase-locked
  loops used to recover the carrier phase. If necessary, a start-up
  procedure, such as one that uses a slow search, can be used for
  initialization. Synchronization is then maintained by locking a
  feedback loop to a clock signal that is extracted from the modulated
  waveform.

Then here I found a diagram which seems to be the kind of algorithm I am looking for:

  The timing tone can be extracted by ... computing DFT at the symbol
  frequency (i.e., a single point of the DFT output is needed for each
  data block)

There must be more literature about this technique? I am not interested in high-performance methods for synchronization, but rather recovering the clock signal by any means without conserving memory or processing time.

I figure the basic technique is to

Pass the signal through a matched filter
"Condition" the signal (and make it purely real) by computing the magnitude (as in the diagram)
Passing the signal through a DFT computation
Searching the DFT bins for the highest power frequency
Use the discovered frequency of the DFT bin to generate a local timing signal
Use the discovered phase of the DFT bin (as in the diagram) to adjust the local timing signal offset


texts:
I would appreciate it very much if someone would be able to provide some clarity on plotting phase responses.

For instance, given that the frequency response of a filter can be written as $$H(\exp(j*\Omega))=\cos(\Omega)$$ How would I go about determining the phase response?

The answer key demonstrates a piecewise solution:

This is what I understand so far:

$\cos(\Omega)$ is strictly real, so the phase carries values of $0$, or $\pm\pi$
Between the values of $-\pi/2$ and $\pi/2$, $\cos(\Omega)$ is positive, so the phase is $0$

I understand since the function $\cos(\Omega)$ is strictly real, phases where the function is not positive should either be $\pi$ or $-\pi$ (as opposed to $-\pi/2$ or $\pi/2$, for when it's imaginary), but how do you know for which range it's $\pi$, and for which range it's $-\pi$?

This might be a stupid question, but I've spent a lot of time thinking it through, so I'm thinking I'm missing something, and would appreciate any clarity.

Thank you so much!


texts: Apply butterworth filter (lowpass) to a signal

texts:
I'm working on a DSP lab and I can't figure out how to apply the filter I've created to a sound. The only thing I found online that did something to my sound was signal.filtfilt but that seems to just output silence for some reason.

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

sf = 44100
t = np.arange(0., 0.5+1./sf, 1./sf)
totalTime = 0.5

mysound = totalTime * np.sin(2*np.pi*697*t) + totalTime * np.sin(2*np.pi*1477*t)

#Creation of the filter
N  = 6    # Filter order
fc = 1000/3000 # Cutoff frequency, normalized
b, a = signal.butter(N, fc)

#Apply the filter
tempf = signal.filtfilt(b,a, mysound)
librosa.output.write_wav('3mod.wav',tempf,sf)

I think I'm either making a mistake in the filter creation, or I just can't find the right way to apply it. Since I'm very new to DSP I can't really understand all the documentation material, so I might just be overlooking something simple in one of the documents I've been reading, but I can't figure out what.

EDIT: it seems to work when I change the filter creation code to

#Creation of the filter
cutOff = 1000 # Cutoff frequency
nyq = 0.5 * sf
N  = 6    # Filter order
fc = cutOff / nyq # Cutoff frequency normal
b, a = signal.butter(N, fc)

So it seems that having the wrong calculation for the Nyquist frequency gave me silence. I'm not 100% sure that .filtfilt is the right (or only) way to deploy the filter.

Leaving the question here in case it helps others.


texts: Do lowpass filters affect the integral over the signal?

texts:
For which filters, does the resulting signal have the same integral (/average) as the original? Obviously this is true for a moving average, but is it also mathematically true for e.g. a Butterworth filter?

I am trying to build a lowpass filter for a signal of measurement samples. While there is a lot of noise, the average (over a sizable number of samples) is correct. I need to filter it retaining that property, regardless of the input frequency spectrum.


texts: Significance of outer product in structure tensor / second moment matrix

texts:
How does one interpret the outer product $\mathbf{M} = \nabla I \otimes \nabla I = \begin{bmatrix}\left(\frac{\partial I}{\partial x}(x)\right)^2 & \frac{\partial I}{\partial x}(x) \cdot \frac{\partial I}{\partial y}(y)\\\frac{\partial I}{\partial y}(y) \cdot \frac{\partial I}{\partial x}(x) & \left(\frac{\partial I}{\partial y}(y)\right)^2\end{bmatrix}$ where $\otimes$ is the outer/tensor product and $I : \mathbb{R}^2 \to \mathbb{R}$ is an image? This matrix is called the structure tensor / second moment matrix.
Background:
The structure tensor tells us how much the pixels change between two image patches (infinitesimal small change $\mathbf{v}$).
\begin{align*}
    f(v_1, v_2) &= \sum_{x,y \in N} \left(I(x, y) - I(x + v_1, y + v_2)\right)^2\\
\end{align*}
where $N$ is a neighborhood. The first order Taylor approximation is
\begin{align*}
    I(x + v_1, y + v_2) &= I(x, y) + \mathbf{v}^T\nabla I(x, y)
\end{align*}
Then
\begin{align*}
    f(v_1, v_2) &\approx \sum_{x,y \in N} \left(\mathbf{v}^T\nabla I(x, y)\right)^2\\
              &= \sum_{x,y \in N} \left(v_1\frac{\partial f}{\partial x}(x) + v_2\frac{\partial f}{\partial y}(y)\right)^2\\
              &= \sum_{x,y \in N} \left(v_1\frac{\partial f}{\partial x}(x)\right)^2 + \left(v_2\frac{\partial f}{\partial y}(y)\right)^2 + 2v_1v_2\frac{\partial f}{\partial y}(y)\frac{\partial f}{\partial x}(x)\\
 &= \sum_{x,y \in N} \mathbf{v}^T\mathbf{M}\mathbf{v}
\end{align*}
I know this is an ellipsoid and there is a relationship to PCA and the covariance matrix. But I am trying to figure out the significance of computing the outer product of the same gradient. It looks a bit like the Hessian matrix $\mathbf{H}(I(\mathbf{x})) = \mathbf{J}(\nabla I(\mathbf{x})) = \begin{bmatrix}\frac{\partial^2 I}{\partial x^2} & \frac{\partial^2 I}{\partial x \partial y}\\\frac{\partial^2 I}{\partial y\partial x} & \frac{\partial^2 I}{\partial^2 y}\end{bmatrix}$ but without the 2nd derivative.


texts: Keep the exact frequency content - Multidimensional fft

texts:
In Matlab I've got to generate a Gaussian process embodying a certain frequency content. Not to go multidimensional, I have used reshape function. I explain better with a piece of code:

dk1 = 2*pi*L/L1;
dk2 = 2*pi*L/L2;
dk3 = 2*pi*L/L3;
k1 = dk1*(-N1/2:(N1/2-1));
k2 = dk2*(-N2/2:(N2/2-1));
k2min = min(k2);
k2max = max(k2);
k3 = dk3*(-N3/2:(N3/2-1));
k3min = min(k3);
k3max = max(k3);
k = allcomb(k1,k2,k3);
k1 = reshape(k(:,1),[1,N]);
k2 = reshape(k(:,2),[1,N]);
k3 = reshape(k(:,3),[1,N]);
k = sqrt(k1.^2 + k2.^2 + k3.^2);
idx = (k1~=0);

allcomb may be found here: http://www.mathworks.com/matlabcentral/fileexchange/10064

Afterwards, I proceed as follows:

k30 = k3(idx) + beta.*k1(idx);
k0 = sqrt(k1(idx).^2 + k2(idx).^2 + k30.^2);

C1 = (beta.*k1(idx).^2.*(k1(idx).^2 + k2(idx).^2 - k3(idx).*k30))./(k(idx).^2.*(k1(idx).^2 + k2(idx).^2));
C2 = ((k2(idx).*k0.^2)./((k1(idx).^2 + k2(idx).^2).^(3/2))).*atan2((beta.*k1(idx).*sqrt(k1(idx).^2 + k2(idx).^2)),(k0.^2 - k30.*k1(idx).*beta));

xhsi1 = C1 - (k2(idx)./k1(idx)).*C2;
xhsi2 = (k2(idx)./k1(idx)).*C1 + C2;

Ek0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));

B = sigma_iso*sqrt((2*pi.^2*L.^3.*Ek0)./(V.*k0.^4));

%% C-matrix Calculation for all those k satisfying the condition |k| > 3

C = zeros(3,3,N);
% In this case, the approximation proposed in Mann[94],Mann[98] and IEC
% 61400-1 3rd Ed. will be used
C(1,1,idx) = B.*(k2(idx).*xhsi1);
C(1,2,idx) = B.*(k3(idx) - k1(idx).*xhsi1 + beta.*k1(idx));
C(1,3,idx) = B.*(-k2(idx));

C(2,1,idx) = B.*(k2(idx).*xhsi2 - k3(idx) - beta.*k1(idx));
C(2,2,idx) = B.*(-k1(idx).*xhsi2);
C(2,3,idx) = B.*(k1(idx));

C(3,1,idx) = B.*(k2(idx).*((k0./k(idx)).^2));
C(3,2,idx) = B.*(-k1(idx).*((k0./k(idx)).^2));

nn = normrnd(0,1,[3,N]);

u = zeros(3,N);
for l = 1:size(C,3)
    u(:,l) = C(:,:,l)*nn(:,l);
end

Then, I need a 3D fft of u; thus, I use reshape once more as

u = reshape(u,[3,N1,N2,N3]);
U = fftn(u);

But I believe that this flow is somehow affecting and corrupting the right frequency content.

Generally, I would go for:

[k1,k2,k3] = ndgrid(k1,k2,k3);
k = sqrt(k1.^2 + k2.^2 + k3.^2);

further to u calculation. The core of the code would switch to 5D instead of 3D.

Do you have any idea how to improve this code? Or may you spot the flaw?

Any support or criticism are welcome.


texts: how does this equation correspond to smoothing?

texts:
Please help me understand smoothing of data. This is a follow up to my previous question posted here. Especially the top answer by Junuxx where he says a way of smoothing a function $f(x)$ is:

$$
f'[t] = 0.1 f[t-1] + 0.8 f[t] + 0.1 f[t+1]
$$

here we can see that for every point in $f[x]$, we are taking a weighted average of that point and its two adjacent points, to get a smoothed version of $f[t]$ called $f'[t]$.

A paper on speech enhancement explains that an equation of the form

$$
y[i] = a[i]y[i-1] + (1 - [i]) x[i]
$$

helps us get the value of y as a recursive smoothing of x. Here $a[i]$ acts as a smoothing parameter and it is itself calculated as

$$
a[i] = \alpha + (1 - \alpha)p[i]
$$

where $p[i]$ is calculated elsewhere and alpha is a constant. $y[i]$, $a[i]$, and $x[i]$ are all arrays with $i$ elements.

How can I relate this equation of $y[i]$ with the equation of $f'[t]$? Both of them are for smoothing data, however equation for $f'[t]$ contains weighted average of consecutive points in the array for $f[x]$ itself while the equation for $y[i]$ does not contain consecutive data points for $x[i]$. How can we comprehend this equation as a smoothing of data in $x$?

If this question is not relevant when the equations are taken out of context then I will be happy to provide more details.


texts: Is it an excessive pre-processing step in image processing?

texts:
I wanted to consult with some experts whether or not this pre-processing step is excessive for what I'm trying to achieve. I am in a similar scenario as the author of the article, where I have many many colonies that I'm trying to segment and distinguish from the background. So I found this article:

A System for Automatic Cell Segmentation of Bacterial Microscopy Images

Preprocessing ->

Power-Law Transformation,
Homomorphic Processing,
Log of Image,
Subtracting Original image by Gaussian Kenrel of Image,
Vincent Dome Transform (never heard of this, so I just used a different de-noising algorithm)
Morphological Operator,
3x3 Median Filter

Then we continue with the main part of the recipe:

II. Background Segmentation

III. Cell Segmentation

IV. Convex Hull Segment Splitting

V. Post Processing

I have no problem going through with this I've already Implemented step 2-6 (not sure how to apply 1). My only concern that it's time consuming to iterate one test for every time I change a value of a controllable variable.

Update

Finally got to ask my professor (an expert of Computer Vision and Machine Learning) and she goes "no this is not excessive at all"


texts:
I would like a simplified measure to quantify spikes in my time series, I have one series that has many and I would like to compare it to another that has very few. I have thought through several options for example fitting a sine wave and looking at the R values as a measure of smoothness etc but I thought that a measure of the Coefficient of variation of the differences in the series may be suitable. I was wondering if this was a good approach and also as the mean values of the differences are very small I was wondering if theres a common alternative for the CV which of course would become meaningless.

Also just to mention that there is no proportionality between the means and standard deviations of the data.

Another way to ask this question would be if anyone had suggestions for the inverse, a successful measurement of smoothness? or a way to summarise jerking in the data


texts: Calculate phase shift between two signals

texts:
Consider the measurements of two modulated AC signals (voltage and current) which may have noise, harmonics and a DC bias. Each measured signal consists of a list for the time and a list for the corresponding values.

My goal is to find out the phase shift in degrees between these two signals in a robust way.

I would like to realize this in the frequency domain perhaps by calculating and evaluating the fundamental waves of the signals.

Any help or sample code (MATLAB) would be very much appreciated.


texts: Lattice structure realization of an all-pole system

texts:
Given an IIR system with the system function,

$$ H(z) = \frac{1}{2 + 1.4z^{-1} + 1.8z^{-2}} $$

I need to draw the lattice structure realization.

I have drawn the lattice structure for $\frac{1}{1 + 0.7z^{-1} + 0.9z^{-2}}$, which is shown in the figure below.

Now since $H(z) = \frac{1}{2(1 + 0.7z^{-1} + 0.9z^{-2})} $, should I multiply the signal $x(n)$ by $2$ or $1/2$ before feeding to the system? My intuition is that I should multiply with $\frac 1 2$, but then since $2$ appears as the multiple in the denominator polynomial, I got confused.


texts: differences in fundamental frequency/pitch between music and speech

texts: I am developping a program, that extracts very expressive frames from a video (meaning: frames, which could summarize a video) for a university seminar. One of my ideas was, to analyze the audio data (e.g.: if it's very loud, there must be sth. happening and so on). I'm using Praat for this analysis. Now, I stumbled upon the pitch, which is supposed to be the fundamental frequency (formant f0). Can anyone tell me what usually happens to the pitch, when somebody talks or when you can hear only music? And do you have any other ideas what the audio data of a film might be able to tell me?


texts: Fast implementation of LOESS for equally spaced samples (possibly using convolution)

texts:
Is it possible to do LOESS smoothing on equally spaced data using convolution?

The Wikipedia article on LOESS mentions a possibility of using a FIR filter for LOESS. I could not find anything more about it, though.

I'm currently using statsmodels.nonparametric.smoothers_lowess.lowess
in Python but it's slow. I tried switching to a Savitzky-Golay smoothing filter but I could not get comparable results on my data.

I mostly use this for detrending.


texts: Kalman for 3D position and 1D orientation

texts:
Is the following state correct for a Constant Acceleration (CA) model KF applied for tracking an object which moves in (x, y, z) and have the freedom of rotation $\phi$ around the $z$ axis?

\begin{equation}
 x=
 \begin{bmatrix}
 x, \dot x, \ddot x, y, \dot y, \ddot y, z, \dot z, \ddot z, \phi, \dot \phi, \ddot \phi
 \end{bmatrix}
\end{equation}

UPDATE:

I can only measure $(x, y, z)$, directly. After, I computed $\phi$, rotation around the $z$ axis, using $\arctan(\frac{\delta y}{\delta x})$, indirectly. The measurements that I input to KF are: $x, y, z, \phi$.

Assuming:
\begin{equation}
\pmb{\alpha}=
\begin{bmatrix}
1  & dt & 0.5 \times dt^2 \\
0  & 1  & dt \\
0  & 0  & 1
\end{bmatrix}_{3 \times 3}
\end{equation}

\begin{equation}
\pmb{\beta}=
\begin{bmatrix}
1  & 0 & 0
\end{bmatrix}_{1 \times 3}
\end{equation}

$A$, the state transition matrix which applies the effect of each state parameter at time $t-1$ on the state at time $t$, is built using:
\begin{equation}
x_{k+1}=
\begin{bmatrix}
\pmb{\alpha}_{3 \times 3}  & 0_{3 \times 3} & 0_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & \pmb{\alpha}_{3 \times 3} & 0_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & 0_{3 \times 3} & \pmb{\alpha}_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & 0_{3 \times 3} & 0_{3 \times 3} & \pmb{\alpha}_{3 \times 3} \\
\end{bmatrix}_{12 \times 12}
\cdot
\begin{bmatrix}
x\\ \dot x\\ \ddot x\\ y\\ \dot y\\ \ddot y\\ z\\ \dot z\\ \ddot z\\ \phi\\ \dot \phi\\ \ddot \phi
\end{bmatrix}_{k}
\end{equation}

In the same manner, the transformation matrix $H$, which map state vector parameters into the measurements domain, is built using:
\begin{equation}
H=
\begin{bmatrix}
\pmb{\beta}_{1 \times 3}  & 0_{1 \times 3} & 0_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & \pmb{\beta}_{1 \times 3} & 0_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & 0_{1 \times 3} & \pmb{\beta}_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & 0_{1 \times 3} & 0_{1 \times 3} & \pmb{\beta}_{1 \times 3} \\
\end{bmatrix}_{4 \times 12}
\end{equation}


texts: Filter difference equation and impulse response

texts:
I am having trouble specifying the transfer function of the following system as seen below:

From the initial view I work out the difference equation to be:

$$
y(n) = x(n) - p^Nx(n-N) +py(n-1)
$$

However should the constant before y(n-1) be negative p according to the following standardised equation for a difference equation?

\begin{align}
y(n)=b_0 x(n)+&b_1 x(n-1)+\cdots+b_M x(n-M)\\
-&a_1 y(n-1)-\cdots-a_N y(n-N)
\end{align}

Furthermore I am trying to work out the impulse response and it is required that the first $N$ values ($n=0,1,2,\ldots,N-1$) of the impulse response of the filter are defined as:
$$
h(n) = p^2 $$

I cannot get this to work from my current workings but feel it may be related to the first problem.

Thanks for taking the time to read my problem and any help would be greatly appreciated!


texts: what is continuous wavelet (cwt) ,wavelet packet (wpt) and stockwell (S-T) computational complexity?

texts:
What is computational complexity of continuous wavelet transform
cwt ?
What is computational complexity of wavelet packet transform wpt ?
What is computational complexity of Stockwell transform transform wpt ?
for the cwt : imagine I have a signal with length of $N$ and I
scale vector is 0:1:s
what will happen if I discretize scale vector to 0:dj:s for better accuracy?
for the WPT : Imagine signal size is $N$ and I am decomposing it in
$L$ levels.

note:

for Continuous case I guess its $O(S\times N^{2})$ and for 0:dj:S, It must be $O(\frac{S}{dj}\times N^{2})$ , Am I right?
for WPT transform it must be $O(2^{L}\times NLogN )$ . Is that right?
For S transfrom I think it is $O(N^{2}Log N)$ . Is it true?
I have calculated the running time of Stockwell and Continuous wavelet transforms. but the time results are strange.
N=3000;            $O(Stockwell)=1.039\times 10^{8}$  it takes 2.5 sec
N=3000, scale=300; $O(cwt)=2.7\times 10^{9}$ it takes 0.58 sec
why this is happening here ? [P.S I have run them 100 times and got the time average]


texts:
How to do this entirely in the time-domain, without using frequency domain?
Let's say I have a continuous-time system that is a differentiator, and X(t) is WSS:
$$x(t) \rightarrow \boxed{\frac{d}{dt}} \rightarrow y(t)$$
I can change the differentiator operator into convolution as follows:
$$X(t) \rightarrow \boxed{\frac{d}{dt} \delta(t)} \rightarrow Y(t)$$
$$h(t) = \frac{d}{dt}\delta(t)$$
$$Y(t) = h(t) * X(t)=\frac{d}{dt}x(t)$$
$$\begin{matrix}R_{XX}(\tau)\end{matrix}\longrightarrow \boxed{\begin{matrix}h(\tau)\end{matrix}}\longrightarrow \begin{matrix}R_{XY}(\tau)\end{matrix}\longrightarrow \boxed{\begin{matrix}h(-\tau)\end{matrix}}\longrightarrow \begin{matrix}R_{YY}(\tau)\end{matrix}$$

cross-correlation of X(t) and Y(t):
$$R_{XY}(\tau) = h(\tau)*R_{XX}(\tau)$$
$$R_{XY}(\tau) = \Big(\frac{d}{d\tau}\delta(\tau)\Big) * R_{XX}(\tau)$$
$$R_{XY}(\tau) = \frac{d}{d\tau} R_{XX}(\tau)$$

autocorrelation of Y(t):
$$R_{YY}(\tau) = h(-\tau)*R_{XY}(\tau)$$
$$R_{YY}(\tau) = h(-\tau)* \frac{d}{d\tau} R_{XX}(\tau)$$
$$h(\tau) = \frac{d}{dt}\delta(\tau)$$
$$h(-\tau) = \frac{d}{dt}\delta(-\tau)$$
$$\delta(\tau) = \delta(-\tau)$$
$$h(-\tau) = \frac{d}{dt}\delta(\tau)$$
$$R_{YY}(\tau) = \Big(\frac{d}{d\tau}\delta(\tau)\Big) * \Big(\frac{d}{d\tau} R_{XX}(\tau)\Big)$$
$$R_{YY}(\tau) = \frac{d^2}{d\tau^2} R_{XX}(\tau)$$
However, the book says the answer should be:
$$R_{YY}(\tau) = -\frac{d^2}{d\tau^2} R_{XX}(\tau)$$
Somewhere I lost a sign...
and the reason why this sign is important because if I take the Fourier of $R_{YY}(\tau)$, it needs this sign to make $S_{YY}(\omega)$ positive as is required for PSD.


texts: My Power Law (Gamma) Transformation Code Doesn't Compatible With Gamma Value

texts:
I'm trying to apply Power Law (Gamma) Transformation. Formula is simple: s = c*rγ.
(c=1 , r= intensity , s= outputintensity). Here is detailed information: https://theailearner.com/2019/01/26/power-law-gamma-transformations/.
I tried code (witch C#) below. But it doesn't work properly.

 public void gammatransform(Bitmap CikisResmi)
    {
        int ResimGenisligi = Img.Width;
        int ResimYuksekligi = Img.Height;
        int Greyscale;
        double y =0.4;

        for (int i = 0; i < ResimGenisligi; i++)
        {
            for (int j = 0; j < ResimYuksekligi; j++)
            {
                OkunanRenk = Img.GetPixel(i, j);
                Greyscale = Convert.ToInt16(OkunanRenk.R * 0.3 + OkunanRenk.G * 0.58 + OkunanRenk.B * 0.12);
                Greyscale = Convert.ToInt32(Math.Pow(Greyscale, y));
                if (Greyscale > 255)
                {
                    Greyscale = 255;
                }
                CikisResmi.SetPixel(i, j, Color.FromArgb(Greyscale, Greyscale, Greyscale));

            }
        }
        return;

}


texts: While finding the ESD/PSD of a signal why we always prefer to find it via Auto-correlation function then the square of the FT of the signal?

texts:
In a video i saw that while calculating ESD or PSD of a signal time auto correlation function was used when it can be also done by the square of FT of the signal.Why we followed that approach even when the later process is more straight-forward and easy.

How does it helps us in later stage?


texts:
For best tree decomposition using wavelet packets, i have tried both Shannon and log energy entropy functions. I noticed that log energy leads to images with higher PSNR values. Could you please help me explaining the reason?

The Matlab expression for the two cases is as following: [Arec,bitstream,PSNR,MSE,D,Drec,s,p_stream]=spiht_wpackets('lena256.bmp',1,'CDF_9x7',5,4,'greedy','shannon'); [Arec,bitstream,PSNR,MSE,D,Drec,s,p_stream]=spiht_wpackets('lena256.bmp',1,'CDF_9x7',5,4,'greedy','logenergy');

Lena256 (image with size 256*256)


texts:
Is there a public database where I can find or anyone can provide me real signals.

I'm looking specifically for pairs of natural signals (e.g. from Sensor networks) to test some cross-correlation and cross-spectrum methods.


texts: Confusion implementing chebyshev type 2 band pass filter in matlab

texts:
I want to design a band pass Chebyshev Type II filter using analog prototyping. The order of filter is 20 with a value of 60 dB stop band attenuation and 0.75 dB pass band ripple where,
Pass band edge = 800 Hz
Stop band edge = 2000 Hz
Sampling frequency = 6000
My Matlab code is below
clc;  close all;   clear all;
Rp = 0.75;    Rs = 60;
fp = 800;
fs = 2000;
Fs = 6000;

fn = Fs/2;     % nyquist frequency
wp = fp/fn;     % normalized Pass band corner frequency
ws = fs/fn;     % normalized Stop band corner frequency

[n,Wn] = cheb2ord(wp,ws,Rp,Rs);
[num,den] = cheby2(20,Rs,Wn);
freqz(num,den,512,Fs);

I have also attached snapshot of output which shows that output  frequency plot looks like low pass,but i think ,why output plot is not band pass?How can it look band pass?


texts: How can I upsample 22 kHz speech audio recording to 44 kHz, maybe using AI?

texts:
I have a 22 kHz mono audio recording, which is mainly speech, a reading. I would like to upsample somehow to 44 kHz, to improve the audible quality. I have read about there are AI methods for upsampling pictures, even videos to a higher resolution. Maybe there is some similar methods for audio too.
The recording comes from radio broadcast, so the recording was excellent studio quality, streamed in 128 kbps, 44 kHz, stereo, but somebody, who did the recording, downsampled it from the original 44 kHz broadcast to 22 kHz, mono. Hence the sound quality is far from optimal, sounds like an old telephone, I think, due to the missing high frequencies.
The sound sample https://drive.google.com/file/d/1AqdQomHjpKqYM4Wq8SdHGOPDR8Irzjj3/view?usp=sharing


texts: Channel estimation using bit reversed order of FFT matrix

texts:
Is it possible to use the bit-reverse order of FFT instead of FFT to estimate the channel?
Based on my reading, I found that bit reversed order of FFT matrix can diagonalize the channel, but unfortunately it can't estimate the channel based on pilot data. It means, that when having $y = h*x$, where $h$ is a channel, $*$ is convolution operation, and  $x$ is a signal. Therefore, in order to get the frequency signal based on bit-reversed order FFT matrix,  $Y = F^b y = F^b (h*x)$. Therefore, the signal $Y = HX$ where the $H$ and $X$ are the frequency-domain channel response and signal gotten using $F^b$, respectively.
Based on my analysis, I got that $H = Y/X$, but when estimating $H$ based on pilot data taken from $X$. For example, $X_p = X_{1:4:end}$, the channel response $H$ cannot be estimated using $H = Y/X_p$. Why ? and is it possible to estimate $H$ based on $F^b$ and $X_p$ ?


texts: When sinusoidal input starts at n=0, why are transient response associated with z-transform poles of digital filter?

texts: In http://www.eng.ucy.ac.cy/cpitris/courses/ECE623/presentations/DSP-LECT-10-11-12.pdf, it says that when sinusoidal input $X(z)$ starts at n=0 (with n<0 having zero input) and the input passes through digitl filter $H(z)$, resulting in output z-transform $Y(z) = H(z)X(z)$, transient response is associated with poles of the filter while steady-state response is associated with the poles of the input. Is this true? If so, how does one prove this?


texts: Downsampling oversampled image (smallest size which retains image 'detail' within acceptable limits)

texts:
I have a large collection of oversampled images (family photos of varying 'sharpness' scanned at a uniformly high dpi) and would like to downsample them to optimise storage space whilst retaining detail within some reasonable threshold.

I am hoping to find (or build) a tool that can help me identify what I (a confessed layman) would call the 'natural resolution' of the scanned image.

For example a 70 year old black & white 6x4" print scanned at 600dpi might be downsampled to 253dpi and look 'sharp' (for some definition of sharp). The downsampled image might then be upsampled to 600dpi again with a deviation from the original that falls within an acceptable threshold. Downsampling the original to 252dpi then upsampling might result in too much deviation, therefore the tool may decide that 253dpi is the sweet spot for that image. Ideally the tool could determine the sweet spot without having to iteratively downsample and upsample.

There are many complications (noise, dust, sharp edges that represent the edges of the photo instead of the photo itself) so the tool may require user supervision. At the moment I'm just looking for pointers as to what techniques and terms I should be familiarising myself with. I have put many hours into googling this issue but all results are either overly simple pc mag stuff (e.g. 'always scan at X dpi') or impenetrable research papers.

Edit (clarification, 2014-01-11)

The paragraph above about downsampling then upsampling is just to illustrate my point about determining loss of detail. In my actual workflow I only intend to downsample (no subsequent upsample step). The goal is to find an algorthm that on a per-image basis will determine the lowest resolution I can downsample to without sacrificing image detail (or sacrificing within an acceptable limit).


texts: Series Sum Operator Properties and Identities

texts:
As I was reading my book, I came across the following expression:
\begin{equation}
\sum _{n=-N}^N\sum _{m=-N}^N \phi \left(n-m\right)
\end{equation}
Then the book said, if we let $k = n-m$, then:
\begin{equation}
\sum _{k=-2N}^{2N}\left(1-\frac{\left|k\right|}{2N+1}\right)\phi \left(k\right)
\end{equation}
Can somebody please tell me if this is a property of summation operator and how they got there. Also please refer me to a reference that carries such properties of summation. Thank you in advance.


texts: How to do Finger Tracking which is occluded by an object?

texts:
Basically I want the hand gesture recognition even when I have an object in hand by using a camera (MS Kinect) in front.

My fingers are not visible directly to the camera due to the object. For example if I could project some electromagnetic rays which can pass through the object but not through the hand so that I get some trace of hand in the image. Even using markers on the fingers can be an option. The material of the object can be plastic or something of this sort.

Any help or reference material would be highly appreciated. Thanks


texts: Blending artifacts in OpenCV image stitching

texts:
I am using OpenCV to blend a set of pre-warped images. As input I have some 4-channel images (*.png or *.tif) from where I can extract a bgr image and an alpha mask with the region related to the image (white) and the background (black). Both image and mask are the inputs of the Blender module cv::detail::Blender::blend.

When I use feather (alpha) blending the result is ok, however, I would like to avoid ghosting effects. When I use multi-band, some artifacts are appearing on the edges of the images:

The problem is similar to the one raised here, and solved here. The thing is, if the solution is creating a binary mask (that I already extract from the alpha channel), it does not work for me. If I add padding to the ovelap between both images, it takes pixels from the background and messes up even more the result.

I guess that probably it has to do with the functions pyrUp and pyrDown, because maybe the blurring to create the Gaussian and Laplacian pyramids is applied to the whole image, and not only to the positive alpha region. In any case, I don't know how to fix the problem using these functions, and I cannot find another efficient solution.

When I use another implementation of multiresolution blending, it works perfectly, however, I am very interested in integrating the multi-band implementation of OpenCV. Any idea of how to fix this issue?


texts: What is the z-domain transfer function of variable delay?

texts:
Is it possible to obtain the transfer function of a signal that has a variable delay?

For instance, I have a signal consisting of 10 samples:
0,1,1,1,1,0,0,0,0,1

If I would apply a unit delay to the signal, the transfer function would be z^-1 and the signal is now:

0,0,1,1,1,1,0,0,0,0

What would be the effect e.g. transfer function if I delay the rising edge of the pulse by 1 delay as above, but the falling edge would be delayed by 2 samples? The signal then would be:

0,0,1,1,1,1,1,0,0,0

Could this be derrived as a combination of z^(-1) and z^(-2) or this transfer function is non-linear?

Thank you for your answers.


texts:
My knowledge of wavelets is less than epsilon. Bear with me.
If I have a signal of two well separated sinusoids (15 and 48 Hz) plus some random noise, I can clearly make out the two in a spectrogram (the two stripes in my picture);

t=0:0.002:1;   % fs = 500 Hz
x=4*sin(2*pi*15*t)+2*cos(2*pi*48*t);
xn = x + randn(size(x));
figure(1);
plot(xn);
figure(2);
spectrogram(xn, 64, 60, [], 500);

but, using 'wscalogram', I cant say I can tell the components apart:

    coefs = cwt(xn,1:64,'db8','scalCNT');
    wscalogram('image',coefs,'scales',1:64,'ydata',xn);

Is there a way to read from the scalogram that there are 2 distinct sinusoids and if so, how can I separate them using wavelet decomposition and filtering? I thought that perhaps my choice of wavelet ('db8') is not optimal but I can't see much difference using other types ad larger (or smaller) number of scales. Clearly I am missing something about where and when and on what to apply wavelets.

Thanks


texts: transfer function of a sampler in the s domain

texts:
I would like to modelize my whole system into the S-domain. This is a mixed system, there a numerical part (corrector, ADC, DAC) and an analogic part (plant transfer function, sensors, etc...). I know that it is possible to transform a DAC (Digital to Analog converter) or a ZOH into the s domain thanks to the Padé's approximation. Nevertheless I do not know how to do for the ADC as it is non linear like the DAC. Is it possible to do the same with the Padé approximation. Does anyone know the transfer function in the S domain or the approximation of the transfer function of the ADC ?
Thank you and have a nice day !


texts:
In an attempt to get myself up to speed with PLLs I'm trying to implement Example 1 from this tutorial, in Python:
https://wirelesspi.com/phase-locked-loop-pll-in-a-software-defined-radio-sdr/
(search for "Example 1")

I have output that is close to what the article shows, except that my cos_out[n]  seems to be "leaning" to the left:

My phase-difference signale_D[n] is flipped vertically from what the tutorial shows:

I must be missing a - sign or have a sin() where I need cos() but I've checked and triple checked my code.

My code:

import numpy as np
k = 1
N = 15
K_p = 0.2667
K_i = 0.0178
K_0 = 1

samples = np.linspace(0, 100, 100, endpoint=False)
input_signal = np.cos(2*np.pi*(k/N)*samples + np.pi)

integrator_out = 0
phase_estimate = []
e_D = [] #phase-error output
e_F = [] #loop filter output
sin_out = []
cos_out = []

for n, sample_in in enumerate(input_signal):

    # phase detector
    try:
        e_D.append(sample_in * sin_out[n-1])
    except IndexError:
        e_D.append(0)

    #loop filter
    integrator_out += K_i * e_D[n]
    e_F.append(K_p * e_D[n] + integrator_out)

    #NCO
    try:
        phase_estimate.append(phase_estimate[n-1] + K_0 * e_F[n])
    except IndexError:
        phase_estimate.append(K_0 * e_F[n])

    sin_out.append(-np.sin(2*np.pi*(k/N)*n + phase_estimate[n]))
    cos_out.append(np.cos(2*np.pi*(k/N)*n + phase_estimate[n]))


texts: $\tt ifft()$ function - absolute vs real form

texts:
I have multiple files from an experiment in frequency-domain that I would like to use ifft() function to convert to the time domain in R to apply signal processing techniques.

My question is about the principle of this conversion:
at the moment my frequency values are non-complex numbers (I believe they are in their absolute form).

So, to convert these columns using ifft() function, should I use the absolute values or there is a way break these values down to the complex form and then do ifft() on those?
My other question is, which form (absolute or real) is the best form to save my data in the time domain for further signal processing techniques?

If this helps, my data looks like this:

| Time (s)| 0.0000 |0.000164|...
|:--------|--------|:-------|:-------
| 1.52kHz |  2747  |  350   |...
| 3.05kHz |  2996  |  420   |...
| 4.57kHz |  4078  |  300   |...
|   ...   |  ...   |  ...   |...

Total of $513$ rows $\times$ $122041$ cols, sampling frequency of $1562.5\textrm{ kHz}$.


texts: Signal Acquisition in Compressed Sensing

texts:
I'm trying to wrap my head around compressed sensing, so I've been reading this intro to the topic.

I'm completely keeping up when they discuss exploiting the sparsity of a signal in some domain to compress it.  I also think I get why you'd want to choose a sensing matrix that is incoherent with the basis in which your signal is sparse.

The point where my understanding breaks down is the "Undersampling and Sparse Signal Recovery" section.

Concretely, let's say I have a signal f of length n.  I have a sensing matrix which is incoherent with the basis in which f is sparse that is m x n, and m < n.

If I understand correctly, I'm supposed to measure the signal by taking its inner product with each row of my sensing matrix, so I'll end up with samples of length m.

As I understand it, compressed sensing is useful in cases when sampling at the Nyquist/Shannon frequency is prohibitively expensive or slow.  I must be missing something really fundamental, but didn't I have to do the work to acquire all samples of f so that it could be correlated with my sensing matrix?  As I understand it, I've compressed the data now, which is great, but I didn't save any energy in the acquisition process.


texts:
I am fairly new to DSP have been reading up on digital filters and have a few questions. I have a signal with very high frequencies whose DC component i want to estimate. The magnitudes of these frequencies can be very high too. I can average out the signal, and this works, but the high frequencies mean high sample rate. I want to do this a lot faster. If i can filter the high frequencies, i should be able to average out the signal with a lower sample rate. To do this i need a low pass filter. If i can design a low pass filter that can attenuate all frequencies above the cutoff frequency (Fcutoff) by a factor of K or more decibels (K can be any number, -1000db, -100000db etc;) then i should be able to estimate the DC component (DC component estimate will be a function of K) much faster.

a) Is it possible, in theory, to design such a low pass filter, where i specify not just Fcutoff but also K. And if i were to implement such a filter, which filter would be the best one. FIR or IIR? Butterworth, chebychef, elliptic etc?

b)I was reading up on the butterworth filter on wikipedia and other online sources. Wiki states that the butterworth filter for high frequencies has an attenuation of 20n db/decade, where n is the filter order. Consider a filter of order 5. If i want all frequencies >= 1Hz to be attenuated by >= -100db, then i can move the cutoff to 0.1hz. And if >= 1Hz is to be attenuated by >= -1000db, i can move the cutoff to 0.01hz. Doing this would still allow the DC component to pass through and filter everything else out. Theoretically this should be possible, right?

c) when you have the magnitude response of a filter in decibles do you use the 10log rule or the 20log rule? I am assuming the 10log rule. Is this correct?

Thanks for your answers!


texts: What is the least JPG-compressible pattern? (camera shooting piece of cloth, scale/angle/lighting may vary)

texts:
I am trying to design a cloth that, from the point of view of a camera, is very difficult to compress with JPG, resulting in big-size files (or leading to low image quality if file size is fixed).

It must work even if the cloth is far away from the camera, or rotated (let's say the scale can vary from 1x to 10x).

Noise is quite good (hard to compress), but it becomes grey when looking from far, becoming easy to compress. A good pattern would be kind of fractal, looking similar at all scales.
Foliage is better (leaves, tiny branches, small branches, big branches), but it uses too few colors.

Here is a first try:

I am sure there are more optimum patterns.
Maybe hexagon or triangle tessellations would perform better.

JPG uses the Y′ Cb Cr color space, I think Cb Cr can be generated in a similar way, but I guess it's better to not use uniformly the full scope of Y' (brightness) since camera will saturate the bright or dark areas (lighting is never perfect).

QUESTION: What is the optimum cloth pattern for this problem?


texts: Difference between $\mathbb{E}[\mathbf{x} \mathbf{x}^{\rm{H}}]$ and $\mathbb{E}[(\mathbf{x}-\boldsymbol{\mu}) (\mathbf{x}-\boldsymbol{\mu})^{\rm{H}}]$

texts:
Let us have a random vector $\mathbf{x} \sim \mathcal{CN} (\boldsymbol{\mu}, \boldsymbol{\Sigma})$ with $\boldsymbol{\mu} \neq \mathbf{0}$. What can we say about the relationship between the elements of $\mathbf{x}$ in the following two separate cases?

1) If $\boldsymbol{\Sigma} \triangleq \mathbb{E}[(\mathbf{x}-\boldsymbol{\mu}) (\mathbf{x}-\boldsymbol{\mu})^{\rm{H}}]$ is a diagonal matrix;

OR

2) If $\mathbf{R} \triangleq \mathbb{E}[\mathbf{x} \mathbf{x}^{\rm{H}}] = \boldsymbol{\mu} \boldsymbol{\mu}^{\rm{H}} + \boldsymbol{\Sigma}$ is a diagonal matrix.

$ $

In 1), the elements of $\mathbf{x}$ are all independent. But what is the interpretation of 2)?


texts: Representing a positive 0-limiting continuous time domain function, in terms of spectral samples, outside of a window

texts:
Bit of an odd question

Assume you have a known function $x(t)$ which is continuous, but approaches but is not equal to, 0, on the interval $[-\infty,\infty]$ . You know also that it approaches 0 exponentially fast on some interval $[T_a, T_b]$ e.g. the Normal or Cauchy Distribution, or some pulse with characteristics you already know with the properties as described. You also know this time domain function's closed form Fourier Transform $X(f)$, which is also continuous:

$$ X(f) \triangleq \mathscr{F}\{ x(t) \} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ dt $$

Say now, you wanted to reconstruct an approximate $\hat{x}(t) \approx x(t)$ from only spectral samples $f_k$ coming from the Fourier domain, and a given $t$. This is perhaps in the form of an IFT of an FT. It is acceptable that this reconstruction of $\hat{x}(t)$ be equal to $0$ outside of the given time domain range $[T_a, T_b]$ at a specified error which is determined beforehand.

I have only dabbled a little in signal processing, but my current understanding suggest that one can only represent this reconstruction $\hat{x}(t)$signal by one of two ways

Since you know the closed form of $x(t)$, you can take a Fourier series with coefficients calculated on the given interval. This will enable one to take the Fourier Transform of this FS, and then apply the Inverse Fourier Transform to gain the desired representation. By using the FS however, it means that your reconstructed signal is inherently assumed to be periodic outside of your given interval, such that if one were to provide a $t$ such that $t<T_a$ and $t>T_b$  you would not get 0, but instead a value, in the period you're in, that corresponds to a value inside $[T_a, T_b]$
Since you know the close from of $X(f)$, you can directly sample in the frequency domain for a given spectral sampling rate. This however would still give the same result as 1. because the sampling impulse train is periodic, while $X(f)$ is continuous, meaning one would need an infinite number of samples, which is definitely not desirable!

The possible idea I had for a solution is to say that $x(t)$ is multiplied in the time domain with $w(t)=\operatorname{rect}(t/T)$ for some period $T$. Assume for simplicity now that our $[T_a, T_b]$ is symmetric about 0. This would mean we would require convolution with a $\operatorname{sinc}$ function in the frequency domain. But since the $\operatorname{sinc}$ function is continuous, our convolved $W(f) \star X(f)$ would then be continuous which then means we would need again infinitely many samples!

I have a feeling something is missing in my understanding but I've searched everywhere for something that hints at an answer, but everything talks about band-limited signals, not time-limited signals, and any that regard time-limited signals only applies to a solution on a specified interval.

So specifically I am curious how to use the additional information (that we know either 1. $x(t)$ is 0 outside a given interval, or 2. we have the closed form of $x(t)$). Also note that this assumed a linear shift invariance for $t$ in $x(t)$. Perhaps an additional approximate integration step between each sample in the frequency domain is required? Or some extension with Poisson Summation?

Any advice is greatly appreciated!

edit:
- thanks Robert for the suggested clarifications


texts:
Forgive me if this is the wrong site for this question. Here's my problem: I am creating a game where players have to search for a hidden detail within the overall image. I downloaded an 11MB compressed photo from the internet which is 85MP in size (14000 x 6000 pixels). I had a very small image within this (a pixel-sized number). The color and brightness of this number closely matches the surrounding pixels so that it blends in. Now, when I went to save this image, the compression destroyed the hidden number, even at JPEG 12/Maximum. This maximum compressed JPEG increased the file size to 50MB even though the original was just 11MB. An uncompressed save file (TIF/PNG) reaches up to 250MB.

Is there anything I can do to compress the image while still preserving the integrity of the number? The number only encompasses a few pixels, but the addition of this small data has increased the file size 25x. I have access to Photoshop and Gimp. Thanks in advance, if you need to see the file then I can provide a link to it (at 250MB size or 50MB size).


texts:
Let $H(e^{j\omega})$ represent the frequency domain response of a minimum phase filter such that $H(e^{j\omega})^{-1}$ exists. Is there a closed-form expression or a bound for the following expression?
$$
\left\lVert H \right\rVert_2^2 \left\lVert \frac{1}{H} \right\rVert_2 ^2
$$
EDIT:
$L^2$-norm is defined as a summation over all frequencies. Obviously, the lower bound will be 1 (if averaging), and equality holds if and only if $H(e^{j\omega})=\text{constant}$. But is there a better bound for a general filter response?
Thanks.


texts: Implementing a Number Theoretic Transform (Fermat number transform--257). Should I consider alternatives to a radix-4 cooley-tukey?

texts:
This is kind of a comp sci question, but I figured I could use some input from FFT experts.

I've already got a radix-4 cooley-tukey implementation of the NTT briefly described on page 9-10 of https://who.rocq.inria.fr/Gaetan.Leurent/files/SIMD.pdf

which was pretty convenient when it came to putting together an Sse2 implementation (parallel computations of 4 32-bit unsigned integers instead of 1 at a time).

In the interest of squeezing out every last bit of efficiency for repeated runs, however, should I consider an alternate FFT implementation (Bruun's algorithm, for instance)?  If you need more context to answer this, I'll be happy to provide more, this is the first time I've ever had to implement an FFT and I'm not sure what information is important.  I'd also appreciate any further explanation you can provide; my understanding is still pretty shallow.


texts:
I would like to make an acoustic measure of one room, so I can compare it with a acoustic measure of some other room. Then compare these data. So when I take an acoustic measure of the first room again, I could recognize that i am in the first room again. Next downside is that I would try to use my iPhone's microphone and the build speakers for that task.
Is this even remotely possible?

Can please someone point me in the right direction.


texts: Why is the number of frequencies decomposed in scipy.signal.stft() equal to the hop size?

texts:
I asked this over on StackOverflow first but figured you fine people might know more about signal processing (duh).

I'm trying to calculate the spectrogram of an audio file in python. There's a nice library function for doing this in scipy.signal.stft(), but it's not behaving as I expect it to.

Since a spectrogram is just a Short-time Fourier Transform, I expect the library function to chop my audio file up into segments ('time frames'), and then do a Fourier transform on each time frame.

This is what the function does as far as I can tell, but the size of the output matrix is behaving strangely. I'll explain:

The parameters you can specify for the function is nperseg and noverlap, which are the number of audio samples in each time frame and how much each time frame should overlap with the previous respectively. From this, we can caluclate the hop size, i.e. how many samples should be between each time frame: H = nperseg - noverlap.

This is where the strangeness comes in. The number of frequencies 'measured' in the Fourier transform of each time frame is exactly equal to the hop size, and it's always evenly spaced between 0 Hz and half the audio sampling frequency.

Say I have an audio file sampled at 22 kHz and I've set the function parameters so that the hop size is, say, 5 audio samples between each time frame. Then the output will show 5 different frequencies evenly spaced between 0-11 kHz, like this: [0kHz  2.75kHz  5.5kHz  8.75kHz  11kHz].

Why only these 5? What if I want to have a hop size of 5 but also see what the Fourier transform looks like for more frequencies? What if I want to have a hop size of 706 samples and 149 measured frequencies between 50 kHz and 65 kHz? Is this impossible?


texts: Why would someone need a Spatial Light Modulator with a maximum phase shift $> 2\pi$

texts:
I am working with Spatial Light Modulators for a project and I was wondering why does the producer think someone would need a maximum phase shift $>2 \pi$ since the device works as a monitor and every phase-element monitored on the device could be reduced between $[0,2\pi]$ using a modulo?

Marko


texts: Compressive Sensing - Location of Non Zero Elements

texts: I am new in compressive sensing and I would like to know if there already exists a deterministic or probabilistic approach to estimate the locations of the $k$ non-zero entries in the $k$-sparse vectore $x$.


texts: What is a $BT$ (Bandwidth-Time) product with reference to modulation?

texts:
I would like to know what does Bandwidth-Time product mean. I understand that Bandwidth ($B$) = 1/Symbol Time ($T$), hence $BT = 1$.

But how can it vary?
What is its significance?
e.g. When we say GFSK is GMSK when $BT = 0.5$, what does that mean exactly?


texts: I am stuck on a problem that only has a center frequency of $\pi/2$ for an IIR filter pole-zero plot of this using 2 poles and 4 zeros

texts: Using a pole zero placement method, design an infinite impulse response band-pass filter with center frequency of $\pi/2$ radians per second. Hint: Use two poles and four zeros.


texts: NAN coming after Band Pass Filtering in MATLAB

texts: I have applied a Band Pass Filter having cut-off frequencies [100 200] Hz and order 4 on a 10 sec signal (done in MATLAB). It is working fine. When I have increased the order to 6, though the filter has been designed but applying it to the same 10 sec signal, the filtered output is showing NAN. Can anyone help me out on this issue? And what is the reason for such behavior in filtering?


texts: Phase distortion of frequency warping by using Praat

texts:
I am trying to do speaker conversion by frequency warping method. I modify the function Spectrum_shiftFrequencies() in Praat by scaling the FFT spectrum of a voice. When one scaling factor is multiplied to spectrum, the timbre of source voice is changed and the voice quality is good. However, using multiple scaling induces several frequency bands of echos into the synthesis audio. I guess that the multiple scaling method breaks the phase information of a spectrum.
Fig.1 single scaling and multiple scaling

Here is the code that I modify in Spectrum_shiftFrequencies() in Praat.

for (long i = 1; i <= thy nx; i++) {  //scan target frequency component
    double thyf = thy x1 + (i - 1) * thy dx;  // target freq
    double myf =  0;                 //  from source freq

    if(thyf>0 && thyf <=target_formant[0]) // 2 scaling
        myf=thyf*source_formant[0]/target_formant[0];
    else
        myf=source_formant[0] + (thyf-target_formant[0])*(my xmax-source_formant[0])/(my xmax-target_formant[0]);

        if (myf >= my xmin && myf <= my xmax) {
            double index = Sampled_xToIndex (me, myf); // index:  source freq index

            thy z[1][i] = NUM_interpolate_sinc (my z[1], my nx, index, interpolationDepth);
            thy z[2][i] = NUM_interpolate_sinc (my z[2], my nx, index, interpolationDepth);
      }

}

The real and image parts of target spectrum are interpolated from source spectrum. After checking the amplitude spectrum of synthesis audio, the contour of spectrum is changed as the expected scaling function. The phase spectrum is also changed, but that induces echo.

Does the phase spectrum require special process in this multiple scaling (nonlinear) conversion? or How to eliminate the echo ? or What is the correct phase conversion of this frequency mapping ?

Any help is appreciated. Thank you.


texts: What technique can I use to remove useless information from an image?

texts:
Is there any way to remove useless information (unimportant details) from an image?

If, for example, we have a specific image, so this contains both useful and useless information. If this image can be recognized directly without the useless information, the unnecessary information can be removed (e.g., illumination, some noises, variations, etc.).


texts: Phase of a signal with multi delayed windowed (rectangular window) signal with different amplitudes

texts:
I have a signal $x(t)$ which is composed of multi delayed linear chirps with different amplitudes and phases but they share the same $w$ and $\alpha$, the signal can be expressed as follows:

$$\DeclareMathOperator{\rect}{\mathrm{rect}}
 x(t) =\sum_{i} a_i \cos\left(\omega t+\frac{\alpha}{2} t^2 +\phi_i\right) \times \rect\left(\frac{t-T_i}{\Delta T_i}\right)
$$

where $a_i$ can be any real number. The signal $x(t)$

differentiated w.r.t   time $t$

My goal is not to estimate the different parameters $a_i$,$\phi_i$,  $T_i$ or $\Delta T_i$ but to estimate the time varying phase term $\omega t+\frac{\alpha}{2} t^2 $.

Any ideas/tips?


texts: Does phase margin exists only for linear systems?

texts: If possible, how to determine the phase margin for nonlinear systems?


texts:
I'm analyzing guitar string plucks and sustains.  I'm having good success with auto correlation using FFT's.  Now I'd like to detect plucks while the string is still vibrating.  Since I already am periodically taking FFT's to find the pitch, I thought maybe I can take advantage of the fft results and look for changes that might indicate a repluck.  Do I simply add up the bins for some kind of power measurement?

Any ideas?


texts:
I am writing an app to work with synthetic time series data from a physics experiment. In our experiments we always have $1/f$ noise in our time series, but I haven't been able to find code/packages to generate $1/f$ noise in synthetic time series data. Can anyone tell me how to do this? I have a feeling it should be pretty simple and that I'm missing something obvious, but I don't know what it is. Any advice?

Also, I looked through the list of questions and don't see anything too on point apart from this: making pink noise (1/f) using list of frequencies


texts:
I was reading "Optical Recognition Systems" book and in chapter 3.3.2.3 where methods for full genetic programming are described, a paper is describing an ADF (automatically defined functions) like this:

  The ADFs have the terminal set Tf ={X, N, S, W, E, NE, NW, SE, SW} and the
  function set Ff = {AND, OR, NOT}. The main program has the terminal set Tc = {I,
  L, NIL} and a set of functions (Fc) containing movement functions, logical operators,
  four ADFs and the HOMING operator.

I am confused about the meaning of HOMING operator. An explanation or some links to point me in the right direction would be great.


texts: Calculating N-point DFT of a signal based on another signal's DFT

texts:
Recently I have asked this answer.

Now I would like to know a little bit more about expressing N-point DFT's of signals in terms of one another.

Having N-point DFT X(k) of a certain signal x(n), how can I calculate N-point DFT of a signal $x_{s}=x(n)+(-1)^n \cdot{}  x(n)$ . Assuming $N$ is even.

Having thought about it a little bit, I came to a conclusion that we cancel half of the samples out and multiply the value of the rest of samples by a factor of 2. The spectrum will now have greater amplitude of peak but I am not sure how it will look like in terms of the location of frequency peak.


texts: Fourier transformed acceleration from Fourier Transform of velocity and position using Omega arithmetic

texts:
Let's say the position of an object is given by simple sine function. By elementary calculus, I can calculate the acceleration in the time domain and find its Fourier transform.
I can also calculate the Fourier transform of velocity and multiply it by frequency or multiply the Fourier transform of position by frequency square. All these three methods should yield the same Fourier transformed accelerations(within numerical accuracy).

Source: http://prosig.com/wp-content/uploads/pdf/blogArticles/OmegaArithmetic.pdf

However, as you can see in the attached plot the magnitude of Fourier transformed acceleration is off by orders of magnitude. Can you please point out what I am doing wrong? I understand that the peak frequency is the same but the magnitudes and overall shape of plots is worrying me.

Here is my code

N =2**8
a =1
t=np.arange(N)
t = np.linspace(-100,100,N)
fs= 2
freq_normalized = rfftfreq(N)*fs

y = np.sin(a*t)
vel = np.cos(a*t)
acc = -np.sin(a*t)

fft_y = np.abs(rfft(y))
acc_y = [a*b**2 for a,b in zip(fft_y,freq_normalized)]

fft_v = np.abs(rfft(vel))
acc_v = [a*b for a,b in zip(fft_v,freq_normalized)]

acc_a = np.abs(rfft(acc))

acc_ifft_v = irfft(acc_v)
acc_ifft_y = irfft(acc_y)
acc_ifft = irfft(acc_a)

fig, ax = plt.subplots(figsize=(12, 8))
ax.spines['left'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_color('none')
ax.spines['left'].set_smart_bounds(True)
ax.spines['bottom'].set_smart_bounds(True)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
plt.plot(t,y)
plt.ylabel(r'$ y $', fontsize=16)
plt.xlabel(r'$ t $', fontsize=16)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from position')
plt.plot(freq_normalized,acc_y)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from velocity')
plt.plot(freq_normalized,acc_v)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from acceleration')
plt.plot(freq_normalized,acc_a)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()


texts: Why is `fftfilt` (i.e. `fft` of both inputs, then element-wise multiplication, then `ifft`) faster than direct convolution?

texts:
I have the following matlab code

N = 500;
M = 32;
K = N + M - 1; % Length of full convolution

x = randn(N);
y = ones(M);

% first one: convolution
r = conv(x, y, 'same');

% second one: multiplication
z = ifft(fft(x, K, K).*fft(y, K, K));
idx = round((K-N)/2+1):round((K+N)/2); % Indices of the central portion
z = z(idx); % Obtaining the central part of the result

May I know which opertation is faster? getting $r$ or $z$? Why is this so? Thanks in advance


texts: How to optimally normalize filter coefficients?

texts:
I have designed an IIR filter in Matlab that comprises the following biquads
0.244   0.002   0.244   1   -1.432  0.543
0.244   -0.345  0.244   1   -1.555  0.777
0.244   -0.388  0.244   1   -1.652  0.943

I want to implement this filter in a 16-bit DSP. Before converting the coefficients to short int, I therefore normalize them by dividing all coefficients by the largest magnitude present in the matrix (in this case 1.652). This results in
b0      b1      b2      a0       a1     a2
0.147   0.001   0.147   0.605   -0.866  0.328
0.147   -0.208  0.147   0.605   -0.941  0.470
0.147   -0.235  0.147   0.605   -1      0.570

I then convert these coefficients to 16-bit values ranging from +32767 to -32768. I implemented the filter using the Transposed Direct Form II.

The I realized that the gain is very low although the filter shows 0 dB in the passband when testing it in Matlab.
The algorithm is as follows. It iterates through all biquads and reintroduces their ouputs to the next one as necessary.
for (k = 0; k < N_sos; k++) {

    y_n = (short) ((s1[k][1] + b[k][0] * x_n) >> 15);

    s1[k][0] = s2[k][1] + b[k][1] * x_n - a[k][1] * y_n;

    s2[k][0] = b[k][2] * x_n - a[k][2] * y_n;

    s1[k][1] = s1[k][0];

    s2[k][1] = s2[k][0];

    x_n = y_n;
}

IIR_out = y_n;

Notice that I have not taken the a0 coefficient into account. Since it usually is 1, it does not matter. But I wonder if the missing gain is caused by this missing coefficient. I tried multiplying y_n with it by exchanging
y_n = (short) ((s1[k][1] + b[k][0] * x_n) >> 15);

with
y_n = (short) a[k][0] * ((s1[k][1] + b[k][0] * x_n) >> 15);

but this heavily distorts the signal.


texts:
I am interested doing the consistency check of the Kalman filter. There are several measures like Posterior Cramer Rao Lower Bound (PCRLB), NIS (Normalised Innovation Squared), NEES(Normalized Estimation Error Squared) autocorrelation, and MSE (Mean Square Error). I had following doubts.

In PCRLB equation, how can I initialize the J matrix i.e. J0 because its a recursive equation which requires an initialization. In a paper Posterior Cramer–Rao Bounds for Discrete-Time Nonlinear, Filtering eq.(29) gives an idea on how to initialize. But I could not understand how they initialized the J0 matrix.

If I do not know the true states and the true covariance matrices (P, Q and R) because I do not do any simulation, is it still possible to calculate PCRLB because one thing is clear for me is that NEES (Normalized Estimation Error Squared) and MSE (Mean Square Error) can not be done for real data. I can never get the exact true states and covariances.

NIS and autocorrelation can be done for real data. Can Monte Carlo runs be done for real data as well because I need to change R and Q at each run, and the rest will be same.

Any help in this context will be highly appreciated.


texts:
Assume I have only access to only one period of a signal that there is a transient effect on the signal as well.

I am pretty sure that there is a particular frequency in the signal.
Is there a way to extract the exact value of DC and magnitude of that frequency by FFT (somehow erasing the effect of transient from the signal).

To make it more clear, check this figure:

I divided my signal into seven sections.

For each section, I want to know the DC value; and magnitude of the particular frequency that I always know exist in my signal.

My problem is that there is a transient effect on all sections that affect the next section for a short time. This phenomenon is much more clear in T6 in the figure.

Is there a way to get rid of the transient in signal processing?


texts:
What does the output signal look like when you insert $2000 Hz$ square pulses to  chebyshev filter with parameters: (Matlab syntax). I tried to do it like this but I don't understand where it goes wrong.

%Making the square pulses

Fs = 1000000;
t = [0:1/Fs:50e-4];
f = 2000;

%square pulse with 50% duty cycle

xt = 1*square(2*pi*f*t, 50);

%Making the filter

Wp = [1000 3000];
Ws = [100 3900];
[n,Wn] = cheb1ord(Wp,Ws,3,20,'s');
[b,a] = cheby1(n,3,Wn,'s');

%Filtering with the coefficient

xt_filtered = filter(b,a,xt);

%Plotting the filtered signal

plot(t,xt_filtered);

But I get something like this which doesn't make sense..?


texts: What's the physical meanning and principle of applying anisotropic diffusion in transform domain？

texts:
I am reading several papers applying anisotropic diffusion in the transform domain like Shearlet, curvelet. All of them have got excellent results.

The physical meaning and principle of applying anisotropic diffusion in space domain is quiet straightforward -- the gradient.

I am confused as to the physical meaning and principle of applying anisotropic diffusion in the transform domain.

Papers, for example:

Shearlet-Based Total Variation Diffusion for Denoising
Glenn R. Easley, Member, IEEE, Demetrio Labate, and Flavia Colonna
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 18, NO. 2, FEBRUARY 2009


texts: How to check if the poles of transfer function are stable?

texts:
$$H(z)=\frac{az}{[z-(b+ja)][z-(b-ja)]}$$

The two poles are looking like this:
$z_{\inf, 1}=b+ja, z_{\inf, 2}=b-ja$

I know they must be inside the unit circle: $|z_{\inf,i}|<1$

So I replace the term: $|b+ja|<1$ by the first pole, but how I must continue?


texts: Switching camera angle based on several images

texts:
Given a number (3 or 4) of different images of a given subject from different angles, how do I interpolate to an angle that isn't available on any of the photos?

What is this technique called and is there any software that does this available? (Preferably Open Source.)


texts: Term for signal waveform correction or demodulation

texts:
Let say a periodic signal waveform $v(t)$ can be modeled as:

$v(t) = f(\phi) A \sin(2\pi ft + \psi)$

where $A$ is a constant amplitude, $f$ is signal frequency, $t$ is time and $\psi$ is a constant phase.

While $f(\phi)$ is $f = k \cos(\phi)$, where $k$ is a constant and $\phi$ is a function of time $t$.

Now the signal $v(t)$ is corrected by removing the dependency to $
\phi$ factor i.e. $v'(t) = v(t)/f(\phi)$.

My question: is the process getting $v'(t)$ from $v(t)$ considered signal waveform demodulation ? or undistortion ?

The generic term I can think safely would be signal waveform correction but it would be better if there is a precise terminology.


texts: Zero Phase Filter: Determining Initial Conditions for Forward Backward Filtering

texts:
Is anybody familiar with Gustafson's algorithm for minimizing transients in forward backward filtering [1]? I'm trying to implement it and my first guess was to check Matlab's filtfilt.m, since they are referencing the paper. In the Matlab function also a linear equation system is solved in order to find initial conditions zi that minimize the startup transients, but the relationship between reference and code is not obvious to me. The only lines of code regarding the minimization are (nfilt is the length of the coefficient vectors):

zi = ( eye(nfilt-1) - [-a(2:nfilt), [eye(nfilt-2); zeros(1,nfilt-2)]] ) \...
 ( b(2:nfilt) - b(1)*a(2:nfilt) );

Can anybody point me in the right direction on how those lines relate to the algorithm described in Gustafson's article?

[1]  Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992.


texts: MVDR Beamformer, minimize but do not change pattern in main beam

texts:
Is there a way to minimize all signals that are not coming from the mainbeam?

I would like to keep the main beam, or at least +-3° around the main peak stable.

When I estimate the covariance from the data itself, it always also minimizes my wanted signal in the main beam. I can not estimate the covariance in the absence of the signal.

Any advice?


texts: How to remove the modulation before doing frequency offset estimation?

texts: To use DFT/FFT (or maximum likelihood) method to estimate the frequency offset introduced by the channel, we need to remove the modulation on the received data samples in the front. If the unknown data samples belong to a QPSK constellation, how should this modulation removal be done?


texts: Optimal windowing strategy for filtering live signal

texts:
I have a stream of 250 samples per second arriving at a pipeline stage which looks something like this simplified example.

class PipeLineStage(...):
    historical_buffer = []
    def _process(self, sample):
        historical_buffer.append(sample)
        altered_sample = ...
        return altered_sample

_process gets called indirectly by the pipeline manager whenever a new sample arrived from the device and the returned value it passed to the next stage.

I want the effect of the function to apply a 50hz notch filter.

So my question is: how should I ideally window/buffer the signal of the incoming samples?

Specifically, How should I:

decide how many samples are optimal to filter out a 50hz signal when the sample rate is 250hz? (with minimal delay for live viewing)
manage the windowing of a continuous signal; worth applying an lfilter over overlapping buffers e.g. samples 1-100, then 50-150, then 100-200 etc? or is it just as good to clear the buffer after each application of the lfilter and have separate/discrete windows, e.g. samples 0-100, 100-200 etc?

In other words, my concern it whether applying a butterfilter with lfilter to discrete windows will introduce discontinuities or other artefacts in the filtered signal significant enough to justify the performance cost of a rolling window.

Does the way I'm thinking about this make sense?


texts: Will using window in Welch's method introduce artificial peak in estimated PSD?

texts:
I'd like to estimate the power spectral density of the signal attached here, the sampling rate is 100 Hz.

If I don't do windowing (i.e., use boxcar window), the result seems to be a simple $1/f^2$ noise:

import json
from scipy import signal
from pylab import *
x = json.load(open('x.json'))
f, p = signal.welch(x, fs=100, nperseg=256, window=signal.get_window('boxcar', 256))
plot(f[1:-1], log(p[1:-1]))

However, if I apply, e.g., hanning window, there will be an unexpected additional peak at 25 Hz:

f, p = signal.welch(x, fs=100, nperseg=256, window=signal.get_window('hanning', 256))
plot(f[1:-1], log(p[1:-1]))

Is this some artifact?


texts: How can i record 1 MHz signal and do a real-time processing in a PC?

texts:
I have a 1 MHz signal from transducer that needs to be analyzed in real-time condition, so the sampling needs to be more than 2 MHz. It consists of digital filtering, turn into spectogram and feed into a model for prediction. I am still stuck at what the device that sufficient for that task. i have tried to use high speed ADC and arduino, but the arduino cannot send the data to the PC rapidly. is there aother method that i could use?

Thanks


texts: How to get the area(Pixels count) inside ring-like shape?

texts:
I have the following image, and I do like to count the pixels inside the ring to get the area.

I did some morphological operations as a kind of post-processing to make the image as much as possible with clear smooth edges.
I tried to do that in different ways as you can see in the code down, but none of them was optimal.
Can you please advice me how to count the pixels inner area of the circle?
Note: some pixels inside are not totally black, they are with low intensity, that's why I was trying to do Otsu thresholding.
Thanks in advance

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread, imsave
# import scipy.ndimage as ndi
from skimage import morphology, filters, feature

seg = np.squeeze(imread('prediction.png')[...,:1])
# meijering alpha=None,
# rem2 = morphology.remove_small_objects(seg, 4)
resf = filters.meijering(seg, sigmas=range(1, 3, 1),  black_ridges=False)

sobel = filters.sobel(resf)
# diam = morphology.diameter_closing(sobel, 64, connectivity=2)
gaussian = filters.gaussian(sobel, sigma= 1)
val = filters.threshold_otsu(gaussian)
resth = gaussian < val

# Morphology
SE = morphology.diamond(2)
# SE = np.ones((3,3))
# SE = morphology.disk(2)
# SE = square(7)
# SE = rectangle(3,3)
# SE = octagon(3, 3)

erosion  = morphology.binary_erosion( resth, SE).astype(np.uint8)
dilation = morphology.binary_dilation(resth, SE).astype(np.uint8)
opening  = morphology.binary_opening( resth, SE).astype(np.uint8)
closing  = morphology.binary_closing( resth, SE).astype(np.uint8)
#thinner = morphology.thin(erosion, max_iter=4)

rem  = morphology.remove_small_holes(resth, 2)

# entropy  = filters.rank.entropy(resth, SE)
# print(seg.shape)

plt.figure(num='PProc')
# 1
plt.subplot('335')
plt.imshow(rem,cmap='gray')
plt.title('rem')
plt.axis('off')
# 2
plt.subplot('336')
plt.imshow(dilation,cmap='gray')
plt.title('dilation')
plt.axis('off')
# 3
plt.subplot('337')
plt.imshow(opening,cmap='gray')
plt.title('opening')
plt.axis('off')
# 4
plt.subplot('338')
plt.imshow(closing,cmap='gray')
plt.title('closing')
plt.axis('off')
# 5
plt.subplot('332')
plt.imshow(seg,cmap='gray')
plt.title('segmented')
plt.axis('off')
# 6
plt.subplot('333')
plt.imshow(resf,cmap='gray')
plt.title('meijering')
plt.axis('off')
# 7
# 8
plt.subplot('334')
plt.imshow(resth,cmap='gray')
plt.title('threshold_otsu')
plt.axis('off')
# 9
plt.subplot('339')
plt.imshow(erosion,cmap='gray')
plt.title('erosion')
plt.axis('off')
#
plt.show()


texts: What are the differences between Gaussian down-sampling and bicubic down-sampling in Matlab? Which is more accurate for simulating low resolution?

texts:
I see in some of the technical papers, good practice for downsampling is to first pass the image through a Gaussian Filter and then take sampling to avoid problems like aliasing and so on.

However, I also see a lot of people simply use the bicubic smoothing function with a scale size below 1 (e.g. 1/2, 1/4) and so on.

What are the differences in characteristics between using these two methods? Which one is more preferable in the specific case of simulating a low res image from a high-res image?

Thanks


texts: Estimation / Reconstruction of an Image from Its Missing Data 2D DFT

texts:
Given the 2D DFT of an image i.e. a NxM matrix of complex numbers, with some missing lines (or even partial lines), considering we have zeros in the missing positions.

Any suggestions for an algorithm to "interpolate"/reconstruct the missing values?


texts: Running an FFT filter on a large data set

texts: I have a non realtime application where I need to run a bandpass FFT filter on a data array of between 5k and 10k data points. Do I break it up into (say) 256 point chunks, run the FFT on that and simply move onto the next 256 items, or is that too naive?


texts:
I have an SNR measure, which is a ratio of two linear functions, and I need to maximize the sum rate of a cellular system given by

$$R = \sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}\mathrm{log}_{2}\left(1 + \frac{f(x_{i,j})}{g(x_{i,j})}\right)$$

subject to the constraints

$$0 \leq x_{i,j} \leq 1$$
$$\sum_{j=1}^{N_{2}}x_{i,j} = 1 \quad i = 1, \dots, N_{1}$$

where, $f$ and $g$ are linear functions in $x_{i,j}$

Can anyone suggest a way (rough algorithm) to solve such a problem? I am unable to figure out how to account for the second constraint.


texts: Computing White Speckle Metric in a Document Image

texts:
The paper 'Metric-based No-reference Quality Assessment of Heterogeneous Document Images', discusses about measuring the quality of characters in a document image. I'm having difficulty to understand the white speckle metric in page 7.

  Small white speckle measures how much fattened character strokes have shrunken
  existing white connected components inside characters, or have created
  new ones by connecting strokes of characters. A histogram of white
  connected components in a document image is computed, and we have
  already found the most frequent font size.Then the white speckle is computed by summing up the histogram bins
  between 1 pixel and 1% of font size squared. The sum is then
  normalized by dividing by the area under the histogram between 1 and
  font size squared.

My questions are:

How a histogram of white connected components in a document image is computed?
How a white speckle is computed by summing up the histogram bins between 1 pixel and 1% of font
size squared? Lets say for example, the most frequent font size is 32, so I have to sum up the frequencies from histogram bin 1 to 32^2 (1024)? Is that right?
Honestly, I dont see any relation of computing or summing up the histogram bins between 1 pixel and 1% of font size squared to the small white speckle measure. Can you help me see the relation?

Thanks.


texts: What's the use case and explanation of frequency-domain subsampling?

texts: In my solution to a certain problem I use subsampling in the frequency domain. Amazingly, it produces good results. Is it really that awkward to see subsampling in frequency domain? What the consequences really are?


texts:
As a non maths-pro, I'm looking for some pointers.

I am rewriting the audio core in my emulator to improve accuracy, and am getting a bit stuck on the specifics of the technique.

I am trying to resample a sample buffer of 4,466 samples, at a sample rate of 223,807.5 Hz (from a simulated Yamaha SN76489), to a sample buffer of 880 samples at a sample rate of 44,100 Hz, which can then be output to a sound card.

Now from what I understand, as the ratio between 44,100 and 223,807.5 Hz is 5.075, I need to turn this into a fraction, which works out as 203/40, and then interpolate by the denominator (leaving 0 samples inbetween), run a low-pass filter, and then decimate by the numerator, giving me my 880 samples of 44,100 Hz output.

(1) Am I along the right lines here? Before I was just downsampling by a straight 5 (i.e. missing out 4 of every 5 samples) and whilst it worked well enough, some notes are slightly out of key in some games, and there is aliasing. Is this new approach correct?

(2) Also, I am not sure what would do the trick in terms of a low-pass filter. I found the following pseudocode on wikipedia that seems to fit the bill:

    // Return RC low-pass filter output samples, given input samples,
    // time interval dt, and time constant RC
    function lowpass(real[0..n] x, real dt, real RC)
       var real[0..n] y
       var real α := dt / (RC + dt)
       y[0] := x[0]
       for i from 1 to n
           y[i] := α * x[i] + (1-α) * y[i-1]
       return y

But I have no idea how to get RC and dt - what are these in relation to the filter I wish to achieve? As I don't know how to select suitable values here, the output signal is heavily garbled at the moment. If a better/more understandable formula is available then by all means suggest it. Sorry for the newbish question and I hope this is the right place for it. Thanks :-)


texts: Cannot understand the attenuation in signal as the audio bit depth of adc changes

texts:
I am working on Active Noise Cancellation Project. I came across this strange problem, where the signal is getting attenuated as the audio bit depth changes and I cannot understand it!
My circutit :I connected Audio Output( Head phone terminal) of PC to Audio IN of NI myRIO 1900 device using AUX cable (two way 3.5mm jack cable). I have a pre-recorded signal(.wav) of sinusoidal waveform with Amplitude 1, Frequency 2kHz, 40000 Sample Rate, 80 Number of Samples, Bit depth( bits per sample) = 16.
Now, I just played .wav file using windows media player and tried to record the waveform on NI myRIO module. Surprisingly, I am getting an attenuation of the signal corresponding to the volume reduction in PC(not linear!!). Only I could see the waveform close to amplitude 1, when I increased the volume of the system to 100%. Please explain this !!

The ADC of NI myRIO has its Resolution as 12 bits but the audio Output(.wav file) is 16 bit data. I am expecting it does not have to matter except in quantization errors, as at the PC Audio Out, the audio data should be converted to analog and and again at the 'NI myRIO Audio In', this analog signal will be sampled! (please correct me if I am wrong !)
Another question is, I have tried the same with sinusoidal sound of Amplitude 2.5 and all other specifications are same as given above, now the audio signal is clipped off at 1V above and below -1V.

 ADC(NI myRIO) nominal range is 2.499V to -2.5V! Is it bacause of the bit depth conversion?
Please explain these things !TQ


texts: What features capture the semantics of a sentence from a gien audio?

texts:
I have been thinking lately and looking up on what could be the possible features that capture the semantics of an audio and came across MFCC: MFCC features represent phonemes (distinct units of sound) as the shape of the vocal tract (which is responsible for sound generation) is manifest in them(layman language).

So is MFCC the only feature that can capture the semantics in an audio?
How well do spectrograms capture the actual semantics of a sentence(audio)? Or are scalograms a better alternative?


texts: How to find and fix the delay of a modulation

texts:
I am creating a model in MATLAB to simulate power line communication, but I found a problem. I'm very new to the DSP area, I get my code below from a MATLAB example.

    %mod/demod GMSK
    hMod = comm.GMSKModulator('BitInput', true, 'InitialPhaseOffset', pi/4);
    hAWGN = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (SNR)','SNR',0);
    hDemod = comm.GMSKDemodulator('BitOutput', true,'InitialPhaseOffset', pi/4);
    % Error calculator by viterbi algortim
    hError = comm.ErrorRate('ReceiveDelay', hDemod.TracebackDepth);
    for counter = 1:10
      % mod 10 words of 8-bits
      data = randi([0 1],8,1);
      modSignal = step(hMod, data);
      noisySignal = step(hAWGN, modSignal);
      receivedData = step(hDemod, noisySignal);
      errorStats = step(hError, data, receivedData);
    end

My question:

How can I discover the delay of my mod/demod process and fix that? Basically what I saw debugging the code was: my output (receivedData) only start to receive the data after the third iterarion of the for-loop.
And another doubt, If I want to store the values of the input data likeFirst data send -> 11110000
First output data in the vector -> 11110000

How should I handle the delay to store the values like above? (Basically is, store the first value send in the first position of my matrix/vector, store the second value send in the second position of my matrix/vector. All of this into the for-loop;


texts: Doubt about wide sense stationary random process

texts:
I have white Gaussian noise $F[n]$ with zero mean and autocorrelation $R_F[n_1,n_2]=\delta[n_1-n_2]$.

If now I consider the random process defined as

$$X[n]=u[n]e^{-kn}F[n]$$ Is $X[n]$ a wide-ense stationary random process for all values of $k$?

I tried to solve the autocorrelation function to check if it depended only on the lag but my result was that $X$ is WSS only for $k=0$, but I think that since it's a sort of a linear filtering it should be that it is WSS for all $k$.


texts: Nyquist Theorem adding two same frequency near to Nyquist Frequency with phase shift

texts:
This is my first question on this platform. Sorry if I made mistakes.

What happens if we add two or more same frequency signals near to Nyquist Frequency with phase shift and sample them?

For example, assume that we have a signal of 18000Hz. And add another signal with same frequency, but having 180 degree phase with our first signal. Then, sample this total signal at 36000Hz. What happens if we reconstruct this signal? Does this situation result in aliasing?

Note: I'am trying to figure out how exactly sound records are sampled when musicians playing their instruments. For example, assume that two keyboard players playing a note at 20000Hz which will be sampled at 44100Hz. Since they are human-being they are not supposed to hold the rhythm like in electronic music. Therefore, there will be a phase shift in their playing and wouldn't this cause aliasing in reconstructed signal?

Thanks for the answers.


texts: Transfer function coefficients and impulse response - FIR filter

texts:
I am currently working on a FIR filter (varying in time). I have the real coefficients of the filter, but I would like to know if they correspond to the coefficients from the impulse response (as $b[k]$) and how I could get the coefficients of the transfer function from it.

The example is that I'm trying to plot the frequency response of the filter in Matlab, but the coefficients doesn't seem the match the definition of filter(b,a,x) with a = 1 and b my coefficients.

It doesn't get me the right result (as I tested them with another FIR function but with only one tap) : I should get after the filtering, values included in $[-20,20]$, but with the basic function filter, I only get results in this range $[-0.4,0.5]$. I'm plotting the results in dB and the shape is exactly the same despite getting the wrong values.

In the definition of the filter function, it says the coefficients should be the ones used in $H(z)$, which I think is not the case.

I'm sorry if I'm missing a big part but I have some comprehension problem about how filters work.


texts:
I'm working on a simple web app that allows the user to tune his/her guitar. I'm a real beginner in signal processing, so please don't judge me too harshly if my question is inappropriate.
So, I managed to get the fundamental frequency using an FFT algorithm and at this point the application is somehow functional. However, there is room for improvement. Right now, I send raw PCM to the FFT algorithm, but I was thinking that maybe there are some pre/post algorithms/filters that may improve the detection. Can you suggest any?
My main problem is that when it detects a certain frequency, it shows that frequency for 1-2sec and then jumps to other random frequencies and comes back again and so on, even if the sound is continuous.
I'm also interested in any other type of optimization if anyone has experience with such things.


texts: finding phase and group delay from frequency response

texts:
(From Schaum's Outlines, DSP, second edition, problem 5.25, second part of problem)

What's the procedure to find the Phase and Group Delay of:

$$
H(e^{j\omega}) = e^{-j\theta}\left( \frac{e^{-j(\omega-\theta)}- |\alpha|
}{1-|\alpha|\ e^{-j(\omega - \theta)}} \right)
$$

Book says that group delay should be:

$$
\tau(\omega) = \tau_0(\omega - \theta) = \frac{1-\alpha^2}{|1-\alpha\ e^{-j(\omega-\theta)}|^2}
$$

I want to say the phase is $\phi(\omega) = (\omega - \theta)$ and group delay is: $\tau(\omega) = - \phi ' (\omega)= 1$... but i know that's not correct because the book gives a more complicated group delay.  it looks like book version has something to do with "$\tau_0$", which i'm guessing is the initial group delay at time zero??  actually, i have no idea... its not explained in the book. (Any ideas how to do this same problem using the method in second part of problem?)

I feel there's a 50% chance i need to do this the complicated algebra and calculus way using this expresss as start:

$$
\phi(\omega) = \arg\Big\{ H_R(e^{j\omega})+jH_I(e^{j\omega}) \Big\}
$$

but then, i'm not really sure if there's a shortcut to just convert it to polar form instead of slugging it out in complex Cartesian form alegbra. like it seems the phase is $\phi(\omega) = (\omega - \theta)$...

with a lot of algebra I obtained this result:

$$
\tau(\omega) = \frac{(-\alpha^2_I - \alpha^2_R + 1)}{1 + \alpha_{R}^{2} + \alpha_{I}^{2} -2 \alpha_{R} \cos(\omega) - 2 \alpha_{I} \sin(\omega)}
$$

still would be nice if there were a simpler way..like the way the second part of the problem suggests...but doesn't really explain...

Fancy factoring to match book result:

$$1 + \alpha_{R}^{2} + \alpha_{I}^{2} -2 \alpha_{R} \cos(\omega) - 2 \alpha_{I} \sin(\omega)$$

$$= (1-\alpha e^{-j\omega}) (1- \alpha^{*} e^{j\omega})$$

$$=|1-\alpha\ e^{-j\omega}|^2$$

other factor:

$$|z|^2 = z_I^2 + z_R^2 $$

$$(-\alpha^2_I - \alpha^2_R + 1)$$

$$=(1 -(\alpha^2_I + \alpha^2_R) )$$

$$=(1 - |\alpha|^2)$$

can remove absolute value as long as $\alpha$ remains squared.

$$=(1 - \alpha^2)$$


texts: Why is the green colour not represented as black in Cr component of YCbCr?

texts:
I am testing conversion from RGB to YCbCr in Python's pillow imaging module (formerly PIL), and am getting predictable results overall, except for green, as represented in the Cr channel component of the YCbCr version of image.

For blue, I get white (RGB (255,255,255); V in HSV of 100%) in the Cb channel, and black (RGB (0,0,0); V in HSV of 0%) for yellow in the same channel. In Cr, I get white for the red colour, but (RGB (21,21,21); V in HSV of 8%) for green. Shouldn't green be black in Cr?

I thought I would check here before claiming there is a bug in pillow. Perhaps I'm missing something about YCbCr, and the green of it is not the same as the RGB green. But I don't think so, as blue, yellow and red do match RGB.


texts: How to infer mathematically if the transmission error is correctable after applying Longitudinal and Vertical Redundancy Check LRC,VRC?

texts:
I have the following message in binary:
Transmitted Message:
1 | 0 | 1 | 1 | 0 | 1
1 | 1 | 0 | 1 | 1 | 0
1 | 1 | 0 | 0 | 1 | 1
1 | 0 | 1 | 0 | 0 | 0
Received Message:
1 | 1 | 1 | 0 | 0 | 1
1 | 0 | 1 | 1 | 1 | 0
0| 0 | 1 | 1 | 1 | 1
0 | 0 | 1 | 0 | 0 | 0
The bits in Bold are the Even Parity bits for both LRC and VRC, how can I infer mathematically if the code is correctable or not? can you please explain to me?
Thanks in advance.


texts: How to make ground-truth image for document image binarization task?

texts:
The two images are from DIBCO (Document Image Binarization Contest). The objective is to find a solution to binarize the top image (original image) to be similar to the bottom one (ground-truth image).

But if I want to make my own ground-truth image for my document image, how do I do it? Any idea, method, or technique is appreciated.


texts: The way to measure chromatic tuners' precision

texts:
I have a question about measuring precision of chromatic tuners.

I want to divide 2 cases.

CASE1:

link; Korg AW-2G's precision

If you check the link or this picture, it says that
Korg AW-2G has the precision of +/- 1 cent.

What's the way or measurement to get this result (+/- 1 cent)?

CASE2:

There are people using Arduino (ATMEL328 MCU) with microphone or piezoelectric diaphragm for the ADC input.

One blog introduced the autocorrelation method
for detecting the instrument's pitch.
After calculating autocorrelation, using parabolic interpolation gives
"accurate" results. (This is not sarcasm, don't get me wrong. :D)

Suppose I made a chromatic tuner using ATMEL328 and a piezoelectric diaphragm.
Also, suppose I chose sampling rate 9240 Hz, 512 samples, the autocorrelation method (with parabolic interpolation).

In this case, how should I define the precision of this chromatic tuner?
I can't buy tuning forks to compare the results.
Let's suppose I'm testing with electrical guitar, classical guitar, and violin.
So any ideas?

Best Regards~