Highpass filter design that is better than the Moving Difference filter The question is actually from the problem set 1.3 of the "Wavelets and Filter Banks" book (authors Nguyen/Strang).
Problem 4: Invent a Highpass filter K with 3 or 4 taps (coefficients) that is better than the moving difference H1. The goal is:
$$
 \left| K(w) \right| < \left| H_1(w) \right| \;\;\; for \;\;\; 0 < \left| w \right|  < \pi/2
$$
And,
$$
 \left| H_1(w) \right| < \left| K(w) \right| < 1 \;\;\; for \;\;\; \pi/2 < \left| w \right|  < \pi
$$
We know that
$$
H_1(ω) = \frac{1}{2} (1-e^{(-jω)} ) \\
\left| H_1(ω) \right| = \left|sin(ω/2)\right|
$$
I'm not sure if the best way to tackle this is by going for an anti-symmetrical filter (4 taps, k1 = -k3, and k2 = -k4). Should I just assume coefficients until I actually meet the condition 1 stated? Or is there another approach? Perhaps setting an equation system to find out about h1 and h2 for K(w)?
Also, does condition 2 mean that for this interval of |W|, the filter needs to be invertible too?
Any help is dearly appreciated.

[ 0.25      -0.6035         0.25       0.1035] seems to work.
Here is how I did that: Since the "breakeven point" is specified at $pi/2$, we know that $|K(0)|^2 = 0$, $|K(\pi/2)|^2 = 0.5$ , and $|K(\pi)|^2 = 1$ . The magnitude of the spectrum is thus given at 4 different points. DC and Nyquist are real and we have a conjugate complex pair at $\omega = pi/2$. We can use this to populate a 4 point frequency domain vector and perform an inverse FFT to get a 4 tab impulse response.
The only thing we don't know is the phase at $\omega = pi/2$. A highpass filter generally has a falling phase and the phase should be zero at Nyquist. Hence we should interpret the phase at DC to be 180 degrees and so it's reasonable to guess the phase at $\omega = pi/2$ half way in between and, $\phi = -pi/2$ seems like a reasonable first guess.
y = ifft([0 .707i 1 -.707i].');

You can play around with the phase and see if it can be made better (whatever that means), but given the amplitude constraints, any 4-tab solution must more or less look this way in the frequency domain.
EDIT
Turns out a phase of 75 degrees gives optimum performance in a sense that it both maximizes the gain at $\omega =3pi/4$ and minimizes the gain at $\omega =\pi/4$


Why use both a high pass and a low pass filter in a Butterworth implementation for noise reduction? I am a new to audio processing .  I want to develop an application that reduces audio noise using a Butterworth filter. I found some existing code doing this, but I still do not understand the use of 2 filters (low and high pass filter) in a reversed order.  What I know from my experience, the high pass filter will eliminates the low frequency samples, so there is no need to apply the low pass filter.
The code begins with :
b,a = scipy.signal.butter(5, 1100/(Frequency/2), btype='highpass') # ButterWorth filter highpass
filteredSignal = scipy.signal.lfilter(b,a,NewSound) # Applying the filter to the signal

and ends with:
c,d = scipy.signal.butter(5, 100/(Frequency/2), btype='lowpass') # ButterWorth low-filter
newFilteredSignal = scipy.signal.lfilter(c,d,filteredSignal) # Applying the filter to the signal

You probably meant to reverse them:  The LPF should be at 1100 Hz, and the HPF should be at 100 Hz.  Then you're keeping everything between 100 Hz and 1100 Hz, and throwing away lower and higher frequencies.

will eliminates the low frequency samples

Also remember filters don't eliminate everything in the stopband, they drop off with frequency, so some of the content outside the bands will still exist.
It will eliminate 0 Hz, but other frequencies will only be attenuated, not eliminated.


why does the spectral envelope of human speech not change w.r.t. pitch when taking a Fourier transform? In the context of speech recognition (recognizing individual speech sounds), the pitch of a certain person can change at different times.

Excerpt from Statistical Signal Processing by Steven Kay:

  This is a natural variability due to the nature of human speech. The
  spectral envelope will not change with pitch since the Fourier
  transform of a periodic signal is a sampled version of the Fourier
  transform one one period of the signal.

What does the part in bold above mean?

The spectral envelope, which determines where the formant frequencies are, is determined by the shape of the mouth, tongue, lips, and nasal coupling.  That is independent of the pitch, which is dependent on the tension or stiffness of the vocal cords.


Laplacian of Gaussian Approximation and Gaussian Blur as the Solution of Heat Equation While I was reading SIFT paper(Lowe, 2004), I came across the method that he apply "heat diffusion equation" to Gaussian function to derive that
$$ \frac{∂G}{∂σ}
= σ∇^2G
$$
I searched Wikipedia and found that original heat equation is

$$ \frac{∂u}{∂t}
= σ∇^2u
$$
I wonder how can time derivative equation be applied to sigma derivative equation. Is there anybody to understand me?

I'm not sure I fully understood what's the issue you're having.
Yet I will show a simple property of the Gaussian filter which might make things clearer.

For simplicity, I will use 1D Signal.
Yet it is easy to extend it into 2D.

The Problem

For $ u \in \mathbb{R} $ (1D, unbounded domain), show that a Gaussian convolution with the initial condition solves the linear diffusion equation:

$$ {u}_{t} = {u}_{xx}, \quad u \left( t = 0; x \right) = f \left( x \right) $$

That is:

$$ u \left( t; x \right) = f \left( x \right) \ast {g}_{ \sigma \left( t \right)} \left( x \right) $$

Where $ \ast $ denotes convolution.

Find the function $ F \left( \cdot \right) $ relating the
evolution time $ t $ to the Gaussian standard deviation $ \sigma \: : \: \sigma = F \left( t \right) $.

Solution

Looking at:

$$ \frac{\partial u \left( x, t \right ) }{\partial t } = \frac{\partial^2 u \left( x, t \right ) }{\partial x^2} , \; u \left( x, t = 0 \right) = f \left( x \right) $$

Defining $ U \left( s, t \right) $ as the Fourier Transform of $ u \left( x, t \right) $ with respect to the parameter $ x $.
This leads to:

$$ \begin{align}
\mathfrak{F} \left( \frac{\partial u }{\partial t} \right) & =  \frac{\partial U \left( s, t \right ) }{\partial t}  \\
\mathfrak{F} \left( \frac{\partial^2 u}{\partial x^2} \right) & = -4 {\pi}^{2} {s}^{2} U \left( s, t \right ) \\
\mathfrak{F} \left( f \left( x \right) \right) & =  F \left( s \right)
\end{align} $$

Hence the equation in the Fourier domain is given by:

$$ \frac{\partial U \left( s, t \right ) }{\partial t} = -4 {\pi}^{2} {s}^{2} U \left( s, t \right ) , \; U \left( s, t = 0 \right)  = F \left( s \right) $$

The solution of this Linear Ordinary Differential Equation is given by:

$$ U \left( s, t \right) = F \left( s \right) {e}^{-4 {\pi}^{2} {S}^{2} t } $$

Hence the solution in the time domain is given by (Using the identity $ \mathfrak{F} \left( {e}^{- \alpha {x}^{2}} \right) = \sqrt{\frac{\pi}{\alpha}} {e}^{- \frac{{\left( \pi f \right)}^{2}}{\alpha}} $):

$$ u \left( x, t \right) = \frac{1}{\sqrt{2 \pi 2t}} {e}^{- \frac{{x}^{2}}{4 t}} = \frac{1}{\sqrt{2 \pi {\sigma \left( t \right)}^{2}}} {e}^{- \frac{{x}^{2}}{2 {\sigma \left( t \right)}^{2}}} $$

Comparing it to the classic Gaussian Function yields that $ \sigma \left( t \right) = \sqrt{\left( 2 t \right)} $.

Remark
This is a solution of mine to a problem I was given in a course named Variational Methods in Image Processing.


Measuring time delay of audio signals Before anyone yells at me I fully realize that this question has been asked numerous times. I assure you that I've read through the existing questions and answers but I am still confused about part of the problem.

I have a sound source that plays music (A) in a closed environment. I have a mic  that I'm using to record  A. I am left with two wav files which share the same characteristics and length (number of samples).

My goal is calculate the time it took for A to reach the mic.

I am trying to perform the calculation using cross-correlation (numpy):

# Delay estimation
corr = numpy.convolve(original_audio, recorded_audio, 'full')
delay = int(len(corr)/2) - numpy.argmax(corr)
distance = delay / sample_rate * 343 # sample_rate == 22050, m/s = speed of sound
print("Distance full: %.2f cm" % (distance * 100))

I consistently obtain values in the 300,000 cm range. The distance between the speaker and mic is approximately 2 feet.

This is all pretty new to me so I'm sure I'm missing something obvious.

Thanks in advance.

Are you sure you shouldn't be using numpy.correlate instead of numpy.convolve?  To estimate delay, you want to cross-correlate your signals, not convolve them.  You'll possibly end up with a much larger delay by convolving.

Trying something simple:

x = [1, 0, 0, 0, 0 ];
y = [0, 0, 0, 0, 1 ];
conv = numpy.convolve(x,y);
conv
array([0, 0, 0, 0, 1, 0, 0, 0, 0])
corr = numpy.correlate(x,y,"full");
corr
array([1, 0, 0, 0, 0, 0, 0, 0, 0])


Estimates sub-pixel shift directly in phase region I try to estimates shift estimation directly in phase region, by following the proposed method in this Sub-pixel Shift Estimation of Image based on
the Least Squares Approximation in Phase Region by Fujimoto, Fujisawa and Ikehara (Proceedings of 26th European Signal Processing Conference, EUSIPCO '16, pp. 91–95. IEEE, 2016 PDF). Here is the flow chart described in the paper

If you're still with me, my issue raised when the author try to proceed with the subtraction of the integer shift with the phase difference $θ(k1,k2)$ (Equation (14) to smooth the phase difference. The authors however did mention that they obtained the slope $(a′,b′)$ of the integer shift using conventional phase-only correlation(POC) (Section III, paragraph above Equation (14)).

How can this step can be done? Since the phase difference is of a $60\times60$ matrix (assume the image dimension is so), while the integer shift consists only TWO values. How exactly were the slopes obtained?

Full matlab code I implemented:

function [ output_args ] = phasecorrlsa( refIm, shifIm )

F=(fft2(double(refIm)));
G=(fft2(double(shifIm)));
[m, n]=size(refIm);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];
R=(F.*conj(G))./abs(F.*conj(G));
r=(ifft2(R));
[ap, bp, rhat]=lsa(angle(r));   %this is my slope a' and b'
[~,w] = max(r(:));
[del_hat2p, del_hat1p] = ind2sub(size(r),w);
del_hat2p=del_hat2p-1;
del_hat1p=del_hat1p-1;
theta=angle(F)-angle(G);
R=exp(1j*theta); %E6
theta=atan2(imag(R),real(R)); %E9
[a, b, thetahat]=lsa(theta);
del_hat1_ts=(m/(2*pi))*a
del_hat2_ts=(n/(2*pi))*b
thetapp=theta-rhat;   %equation 14, I guess something amiss here
%figure;
%surf(theta);
%figure;
%surf(thetapp);
[app, bpp, thetapphat]=lsa(thetapp);
del_hat1=(m/(2*pi))*app+del_hat1p;
del_hat2=(n/(2*pi))*bpp+del_hat2p;
if del_hat1>n/2, del_hat1=del_hat1p-m; end
if del_hat2>m/2, del_hat2=del_hat1p-n; end
output_args=struct('a',del_hat1, 'b', del_hat2);

end

function [ a, b, hat ] = lsa( theta )

[m, n]=size(theta);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];
B=regress(theta(:), X);
a=B(1);
b=B(2);
hat=reshape(X*B,m,n);

end

Please lighten me on this issue. Thank you so much!

I am unable to comment on your Matlab code, but eq. (14) seems straightforward to me.  You have a phase difference field $\theta$ which depends on the two spatial wavenumber components $k_1$ and $k_2$, and which formes a wrapped plane.  The shift you seek is the slope of this plane, expressed as two scalar components $a$ and $b$.  Since phase wraps, the slopes are decomposed into integer and fractional parts, $a = a' + a''$ etc.  In eq. (14) you obtain the residual phase field by subtracting the phase due to the integer part $(a', b')$.  You need the actual $k$ values when you compute $(a'k_1 + b'k_2)$, not just the $k$ indices.  Computing the $k$ values is similar to making a frequency axis for a 1D Fourier spectrum analysis.


MATLAB Plot of FT(Cos) Displays Weird Impulse It is known that:
$$
\mathcal{F}\{\cos(2\pi t)\}=\frac{\delta(f-1)+\delta(f+1)}{2}
$$
However, on MATLAB, I used F=fftshift(fft(x))/N;  to obtain the FT of $\cos(2\pi t)$ and I obtained the following:

My question is why does each $\delta$ appear to be made up of two thin close parallel lines. Note that I have used $1024$ sample points and constructed $f$ a vector from $-51.2$Hz to $51.2$Hz and $t$ a vector from $-5$ to $5$.

That's because Matlab preforms FFT, the formula that you referred to is the FT of your signal, the difference between the two can be derived by convolving the DTFT (which has the same form of the FT) with the sinc function, which in your case is close to a delta function.  For that reason you get a slight different frequancy domain from the theory you would expect


Amplitude loss during resampling with IIR I have a problem in that I am resampling a signal and loosing amplitude. The steps I do are the following, in c++:

Upsample/interpolate signal by 16x by padding with zeros.
Run an IIR Bessel 2nd order low pass filter with cutoff frequency of downsampler.
Downsample/decimate the signal by 13x (I take every 208th data point to achieve this).

The problem I am having is that I get a large amplitude loss after running the Bessel filter and lots of ringing. If I change the cutoff frequency of the filter to match the new sampling rate of the upsampler then I get a clean signal, as I wanted, but I still get the large loss in amplitude. I am loosing around 13x to 23x in intensity. I figure the zero padding is reducing the signal a lot. I'm really stuck on this and haven't done this type of processing before so any help is appreciated.

Try adding a gain of 16 before your low pass filter, or equivalently using a low pass filter with a passband gain of 16 instead of 1.


Determine if an impulse is minimum phase and compute a zero-phase equivalent The impulse is given as a sequence: $x: {10; -3; -1}$. To determine if it is minimum-phase I compute the $Z$-transform (the geophysical case, positive powers of $z$) of it: $X(z)=10-3z-z^2$. Then I factor $X(z)$ into dipoles $(a_0+a_1 z)$: $X(z)=10-3z-z^2=(5+z)(2-z)$. Then I observe that all the dipoles satisfy the $|a_0|>|a_1|$ thus being minimum phase, so the whole sequence $x$ is minimum phase.

Here I'm stuck about the zero-phase equivalent. I can evaluate the spectrum and it's magnitude by evaluating $X(z_k)$ at $z_k=e^{2 \pi i \frac{k}{N}}$, where $k=\{0,1,\dots,N-1\}$, and $N$ - length of the impulse.

From this point on how can I find the respective zero-phase impulse that has the same spectrum magnitude $|X(z_k)|$?

To find the zero-phase impulse response, simply take the Fourier transform of your signal, find magnitude and find the inverse Fourier transform of that. Mathematically speaking,

$$\mathscr{F}^{-1}\{\text{abs}\left(\mathscr{F}\{x[n]\}\right)\}$$

This essentially says that the entire Fourier transform of your new sequence is real-valued. All real-valued numbers have zero phase.


How to deconvolve dependent part of signal from independent part? I have a problem of the following form.  There are two signals, x(t) and y(t).  The model for the system is such that:
$$x(t) = x'(t) + f(y(t))$$
where $f(y(t))$ is a variable offset introduced by interaction of $x$ with $y$.  The function $f$ is unknown.  As an additional constraint, $x'$ is expected to vary slowly and $y$ varies rapidly (and slowly as well).  $x$ has both slow and fast components as a result.  $f$ can probably be approximated pretty well by a linear function: it is monotonic increasing and fairly smooth.

The question is how to estimate $f$ and $x'$ given $x$ and $y$?

Sample data:

t = linspace(0, 100, 100)
y = sin(pi * t / 5) + t/50 + 1
x = sin(pi * t/100) + 0.05 * y * y

x in blue, y in green, xprime in red

First of all note that there is a certain ambiguity in the problem formulation because the approximation of $x(t)$ has to be divided into the part represented by $x^{\prime}(t)$ and the part represented by $f(y(t))$, and this division is not unique. The approach I chose is to start with a parametrization of $f(y(t))$, and approximate the given $x(t)$ as well as possible by this function. The remaining approximation error is then taken care of by $x^{\prime}(t)$. I use a least squares criterion for approximating $x(t)$ by $f(y(t))$.

I use the following parametrization for $f(y(t))$:

$$f(y(t))=a_1y^2(t)+a_2y(t)+a_3\tag{1}$$

with unknown parameters $a_i$, which need to be determined given the known functions $x(t)$ and $y(t)$.
Of course other choices for $f(t)$ are possible, but let's for illustration purposes use the above parametrization.

First you need to define a grid of values $t_i$ and set up an overdetermined system of equations:

$$x(t_i)=x^{\prime}(t_i)+a_1y^2(t_i)+a_2y(t_i)+a_3,\quad i=0,1,\ldots N-1\tag{2}$$

with some $N$ (which must of course be much greater than 3, the number of unknowns). Define now a $N\times 3$ matrix $\mathbf{Y}$ with rows $Y_i=[y^2(t_i), y(t_i),1]$, a $N\times 1$ vector $\mathbf{x}$ with elements $\mathbf{x}_i=x(t_i)$, and an unknown $3\times 1$ vector $\mathbf{a}$ with elements $a_i$. The optimal least squares solution for $\mathbf{a}$ is given by

$$\mathbf{a}=(\mathbf{Y}^T\mathbf{Y})^{-1}\mathbf{Y}^T\mathbf{x}\tag{3}$$

which in Matlab or Octave is simply implemented by the command

a = Y\x;

The values $x^{\prime}(t_i)$ are then simply given by

$$x^{\prime}(t_i)=x(t_i)-[a_1y^2(t)+a_2y(t_i)+a_3]\tag{4}$$

with the parameters $a_i$ as determined from (3). Note that $x^{\prime}(t)$ as defined by (4) is zero-mean. This is no restriction of generality because the parameter $a_3$ defining the constant part of $f(y(t))$ can as well be added to $x^{\prime}(t)$. This is part of the ambiguity inherent to the problem.

Using your example, the corresponding Matlab/Octave code would be something like

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
a = Y\x;
f = a(1)*y.*y+a(2)*y+a(3)*ones(100,1);
xprime = x - f;
plot(t,x,t,xprime,t,f)

From which you get the following plot ('x' blue, 'xprime' green, 'f' red):

As you can see, 'xprime' is zero mean, but his could easily be changed by adding the constant term determined by $a_3$ to 'xprime' (and subtracting it from 'f').

UPDATE:

So far no constraints have been imposed on $x^{\prime}$. If smoothness of $x^{\prime}$ is desired, this can be accomplished by adding more equations to the overdetermined system (2). These equations can be of the form

$$w\left[a_1(y^2(t_i)-y^2(t_{i-1}))+a_2(y(t_i)-y(t_{i-1}))\right]=w[x(t_i)-x(t_{i-1})],\\\quad i=1,\ldots,N-1\tag{5}$$

with some positive weight factor $w$ determining the relative importance of these additional smoothness constraints. If the matrix $\mathbf{Y}$ and the vector $\mathbf{x}$ are augmented by the values given in (5), the overdetermined system can be solved as before by solving equation (3). I computed the solution for $w=10$. Furthermore I added the constant term $a_3$ to $x^{\prime}$, which is an arbitrary choice (you might as well attribute the constant term to $f(y(t))$). The result is shown below (same functions and colors as in the above plot):

The plot was created in Octave by the following code:

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
w = 10;
y2 = y(2:100); y1 = y(1:99);
x2 = x(2:100) - x(1:99);
Y2 = w * [ y2.*y2 - y1.*y1, y2 - y1, zeros(99,1)];
Y = [Y;Y2];
X = [x;w*x2];
a = Y\X
f = a(1)*y.*y+a(2)*y;
xprime = x - f;
plot(t,x,t,xprime,t,f+a(3))


Convolution in Continious time Dear all,

Please find attached image,
I am not able to understand why is the underlined term is to be included.
Shouldn't it be simply the exponential term?

The $u(t-2)$ term needs to be included because the upper limit of the definite integral is $t-2$, not $\infty$:

\begin{align*}
h(t) &= e^{-(t-2)}u(t-2)\\
y(t) &= x(t) * h(t)= \int_{-\infty}^{\infty} x(t')h(t-t') dt'\\
 &= \int_{-\infty}^{\infty} x(t')e^{-(t-t'-2)}u(t-t'-2) dt'\\
u(t-t'-2)&= \begin{cases} 1 & t'<t-2 \\ 0 & \text{otherwise} \end{cases}\\
y(t) &= \int_{-\infty}^{t-2} x(t')e^{-(t-t'-2)}(1) dt' + \int_{t-2}^{\infty} x(t')e^{-(t-t'-2)} (0) dt'\\
&= \int_{-\infty}^{t-2} x(t')e^{-(t-t'-2)} dt'.
\end{align*}


Why doesn't JPEG use 1D DCT for image compression? I know that JPEG uses 2D dct and splits the images in 8x8 blocks. Why doesn't it simply split the image in one-dimensional vectors in $\mathbb{R}^{64}$? Wouldn't it simplify the math? My guess is that this way doesn't take advantage of pixel correlation. Is this right or is there any other reason?

If JPEG used a 1d DCT it would only be eliminating the redundancy in a single direction. You could compress the image with a 1d DCT then turn it sideways and compress it again. Since the DCT is separable doing this double 1d compression procedure would be (roughly) the same as just using a 2d DCT transform to begin with.


Some subtleties in the PSD of digitally generated noise with $S_X(f) = A(f)$ I understand that the title is not the most descriptive, but I didn't know how to describe my problem in a more literal way. The situation is as follows:

I wanted to digitally generate noise with an arbitrary power spectral density $S(f) = A(f)$. To do this, I generated a white noise series $X$ with $S_X(f) = 1$ by simply taking $n$ independent samples from a Gaussian distribution with $\sigma = 1$ and then digitally filter them with filter amplitudes $F(f) = \sqrt{A(f)}$. One can easily show that this leads to $S_Y(f) = A(f)$, where $Y$ is the series resulting from filtering $X$ with $F(f)$.

I have done this in practice using Mathematica's FrequencySamplingFilterKernel and ListConvolve. For this to work one needs to know the desired $A(f)$ of course, which in my example I wanted to be a Lorentzian of center frequency $f_0$ and full-width-at-full-maximum b given by
\begin{equation}
A(f) = \frac{1}{1+\left(\frac{f-f_0}{b}\right)^2}\end{equation}

So I took this equation, sampled it, created filter amplitudes and filtered the white noise. The resulting PSD $S_Y(f)$ looked as desired when analysed with Welch's Method.

But now I am interested in a subtlety in the above. As has been pointed out the me, the above Lorentzian is not a proper power spectral density of a real signal; $A(f) \neq A(-f)$. But the method still worked, in the sense that I did not get any complex valued data, and that the periodogram showed the correct spectrum at positive frequencies. So my question is, using this method, what happened at the negative frequencies? Did I essentially create a PSD for which
\begin{equation}
A(f) = \frac{1}{1+\left(\frac{f-f_0}{b}\right)^2} + \frac{1}{1+\left(\frac{f+f_0}{b}\right)^2}\end{equation}
or perhaps instead it is more like
\begin{equation}
A(f) = \frac{1}{1+\left(\frac{f-\textrm{sgn}[f]f_0}{b}\right)^2}\end{equation}
where $\mathrm{sgn}[\dots]$ is the sign function.

Those are the only reasonable scenario's I can think of, resulting from the above. They are very similar in behaviour, except around $f=0$, Here the second one has a discontinuous derivative at $f=0$; is that even physically possible?

My question is thus what PSD I ended up with; one of the two I mentioned above or even a third one?

You are only looking at the positive frequencies. The negative frequencies are also there. So yes, the spectrum is indeed (with some abuse of notation)

\begin{equation}
A(f) \approx \frac{1}{1+\left(\frac{f-\mathrm{sgn}(f)f_0}{b}\right)^2}\end{equation}

such that the PSD is mirrored with respect to $y$-axis.

I don't speak Mathematica, but I can tell you how it works. Your FrequencySamplingFilterKernel is in fact a simple frequency sampling approach to build FIR filters with arbitrary-magnitude frequency response (similar to Matlab's fdesign.arbmag). Given you have $N$ points of the frequency response, it essentially calculates an impulse response (by taking the inverse  DFT) whose Fourier transform is equal to your desired spectrum when evaluated at the corresponding frequencies. That is:

$$h[n]=\frac{1}{N}\sum_{k=0}^{N-1}H(k)e^{j\frac{2\pi}{N}k}$$

where $H(k)$ is the desired DFT (That is, $N$ samples from $\sqrt{A(f)}$).

Notice that the default impulse response $h[n]$ in the documentation of FrequencySamplingFilterKernel (under details and options) is an odd-length FIR filter with even symmetry. From this point everything is just regular FIR filtering with $h[n]$, and when a real process is filtered it will generate a PSD with even symmetry:
$$S_Y(f)=|H(f)|^2S_X(f)$$
where $|H(f)|^2$ is equal to $|A(f)|$ at the given $N$ desired points.

If you want to do it in Matlab:

N = 100;
f = 0:0.01:1;
b = 1;
f0 = 0.75;
A =sqrt(1./(1+((f-f0)./b).^2));
d = fdesign.arbmag('N,F,A',N,f,A);
hd = design(d,'freqsamp');
w = -pi:.001:pi;
freqz(hd,w);
x=randn(1000000,1);
y=filter(hd,x);
psd(y);

with low-order FIR approximation (N=8):


how to find breakpoints of a signal based on the frequency of the signal There are 4 cos signals with 10,20,30 and 100 Hz each.
The main signal is the sum of those 4 and I have to create a low pass filter, on the main signal for the frequency's > 50 Hz.

% design of low pass filter
f = [0 0.6 0.6 1];
m = [1 1 0 0 ];
b = fir2(30, f, m);
[h, w] = freqz(b,1,128);

hold on;
plot(f,m,'b')
plot(w/pi, abs(h), 'r')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude')
legend('Ideal', 'fir2 designed')
legend boxoff
title('comparison of frequency Response Magnitytes');
hold off;

That's the example code that our professor gave us.On maltab using help fir2 there is a similar example

% Example 1:
%   Design a 30th-order lowpass filter and overplot the desired
%   frequency response with the actual frequency response.

f = [0 0.6 0.6 1];      % Frequency breakpoints
m = [1 1 0 0];          % Magnitude breakpoints
b = fir2(30,f,m);       % Frequency sampling-based FIR filter design
[h,w] = freqz(b,1,128); % Frequency response of filter
plot(f,m,w/pi,abs(h))
legend('Ideal','fir2 Designed')
title('Comparison of Frequency Response Magnitudes')

What i don't understand is how to find the frequency breakpoints and the magnitude breakpoints?
Also whats the relation of the filter order with the low-high, band filters?

The samples in the 'f' vector in the Matlab Example 1 code are arbitrary, ...just an example of a lowpass filter. The purpose of that code is to merely show you how the "designed" (actual) freq response is similar to, but not equal to, the "ideal" rectangular  freq response.

I don't like that Matlab Example 1 example. It's confusing to some beginners to have two different magnitude ('m') sample values associated with a single 'f' value of 0.6. A better example would define 'f' as f = [0 0.6 0.7 1].

In any case, the Example 1 example's f = [0 0.6 0.6 1] and m = [1 1 0 0], vectors specify that the lowpass filter's "ideal" passband gain is one at $0*f_s/2$ Hz and is one at $0.6*f_s/2$ Hz. (The '*' symbol means multiply.) So the "ideal" lowpass passband extends from $0*f_s/2$ -to- $0.6*f_s/2$ Hz. The code also specifies that the lowpass filter's "ideal" stopband gain is zero at $0.6*f_s/2$ Hz and is zero at $1*f_s/2$ Hz. So the "ideal" stopband extends from $0.6*f_s/2$ -to- $1*f_s/2$ Hz.

Now if you wanted the filter's "ideal" lowpass passband to extend from $0*f_s/2$ -to- $0.5*f_s/2$ Hz (zero Hz –to- $f_s/4$ Hz), you would define vector 'f' as f = [0 0.5 0.5 1].


What windowed Sinc might be closest to a Parks–McClellan low-pass filter? In another question on dsp.stackexchange, a statement was made that the output of MATLAB firpm() or firls() for an anti-alias low-pass filter was "close" to the coefficients of a windowed-Sinc low-pass filter.

If so:

What window (applied to a Sinc "perfect" reconstruction waveform) would provide low-filter kernel coefficients closest (In least squared error or other commonly used metric) to those of a same-length FIR filter designed by the Parks–McClellan or Remez-exchange methods?  Is it one of the well-known named window functions (Nuttall et.al.)?  Or is there a good numerical or tractable equation approximation to this window?

none.

what you're asking for is what window, which multiplies a sinc() function, will get you an optimal impulse response kernel, like you might get from Parks-McClellan.  well, to get that window, you must undo the multiplication.

the problem is that you cannot divide by zero.  any of these "sinc-like" functions that come out of P-McC will pass through 0.  but not necessarily at the same place that the sinc() function goes through zero.


impulse response VS zero-input response I am new in the field of systems and signals, and I have a rather basic for the majority of the group, question:

Can we find the impulse response function of homogeneous ODE, instead of its zero-input response?

for example, we have the following 2nd order homogeneous ODE:

$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = 0$

, where the output is the $x(t)$ given the initial conditions

I understand that if it were:
$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = f(t)$

the $f(t)$ would be its input, $x(t)$ its output given the input, and we could find the impulse response by replacing the $f(t)$ with $\delta (t)$.

Now that the input is zero, how can we find what the output would look like with respect to any input?

Is this:

$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = \delta (t)$

even allowed, for an initially homogeneous equation?

Basically, you can solve this kind of EDO using two different techniques (which also have c different meaning).

Classical Method - you will need the full EDO, containing input and outputs with initial conditions applied over the whole solution y(t), where y(t) = yhomogeneous(t) + yforced(t). The first term, you solve the homogeneous EDO considering NO CONDITIONS. Then, following the pattern of input you exactly find the forced contribution. Then, you add both contributions and apply the initial conditions...
In this technique you totally lose the feeling about the system itself.

zero state response and zero input response. In this case, the whole solution is written as y(t) = y0(t) (zero input response) + x(t)h(t)  (zero state response, as convolution of both signals. The text cutof the asterics).
To find y0(t) you do not consider the input x(t), and apply condition at t= 0^-. It happens because it is completely independent of the input. The system responds just because of the storage energy.
The second term, zero state response, it is the convolution of input and impulse response of the system. So, you need to find the impulse response of EDO. It is done considering the whole EDO , e.g., D4y -D3y +2D2y -2Dy +y = Dx -3x . First you have to solve the homogeneous equation and apply some special conditions generalized... y=Dy=D2y=0 and D3y=1. Lets call it yn(t).... Then you apply the input operator over this yn(t). It means,
D (yn(t)) -3yn(t)... After that, you finally have h(t)...
So, you do the convolution of h(t) and x(t) and find the zero state response
y


Butterworth band pass filter I am trying to build a Butterworth bandpass filter. My intent is to have a pass band between $250\textrm{ Hz}$ and $1000\textrm{ Hz}$. However the plot doesn't look right, as can be seen from the following graph:

Below is my code:

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

#apply bandpass filter
fs = 8000
fso2 = fs/2.
N,Wn = scipy.signal.buttord(wp=[250./fso2,1000./fso2], ws=[200./fso2,1200./fso2],
   gpass=0.1, gstop=30.0)

b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(0, 4, 500))
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Butterworth bandpass filter fit to constraints')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.show()

You need to specify your filter design specifications parameters consistently for either an analog or a digital filter. With your posted code, the butterord computes the required order for a digital filter with cutoff frequencies near 1 (which would make sense as Nyquist-normalized cutoff frequencies, but not so much in Hz), then uses those directly to obtain analog filter coefficients.

Depending on your application, you can either design an analog filter with

N,Wn = scipy.signal.buttord(wp=[250.,1000.], ws=[200.,1200.],
   gpass=0.1, gstop=30.0, analog=True)
b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(0, 4, 500))

which should then produce a graph which looks like:

Alternatively you could design a digital filter with

N,Wn = scipy.signal.buttord(wp=[250./fso2,1000./fso2], ws=[200./fso2,1200./fso2],
   gpass=0.1, gstop=10.0, analog=False)
b, a = signal.butter(N, Wn, 'band', False)
w, h = signal.freqz(b, a)

Note however that the original specifications for a digital filter result in an 18th-order filter which may start to show some artifacts due to numerical accuracy. As such I've reduced slightly the gstop parameter to produce the following graph:


Aliasing at $f_0 + kf_s$ In this question, it has been proved that $x(t) = \sin(2\pi f_0 t)$ and $x_k(t) = \sin(2\pi (f_0 + k f_s) t)$ have the same sample points. So sinusoids with frequencies $f_0$ and $f_0 + kf_s$ will identify as the same signals. In other words, we keep sampling frequency constant and increase frequency by $kf_s$. Apparently, only frequencies $f_0 + kf_s$ lead to aliasing. Why this frequency is special? I tried to see this result in the frequency domain but it wasn't useful. I know by sampling in time domain we get shifted replicas of the original spectrum which is two delta functions at $f_0$ and $-f_0$ but how this corresponds to $f_s$ and aliasing at $f_0 + kf_s$? What I'm trying to get here is a visual explanation of mentioned result. In the case of the sampling theorem it's really easy to see aliasing occurs when there is an overlap.

You've made an error in this statement: "I know by sampling in time domain we get shifted replicas of the original spectrum which is two delta functions at $f_0$ and $−f_0$ but how this corresponds to $f_s$ and aliasing at $f_0+kf_s$?"
The two delta functions at $-f_0$ and $f_0$ are the result of taking the Fourier transform of a sinusoid $x(t)$:
$$ x(t) = \sin(2{\pi}f_0t)$$
$$\mathcal{F}(x(t)) = \frac{1}{2i} \big(\delta(f - f_0) + \delta(f + f_0) \big)$$
We haven't introduced sampling yet, this is simply the result of taking the Fourier transform. When we start sampling, then we introduce replicas of the signal's spectrum.
Let's assume we sample the $x(t)$ above at a frequency of $f_s$. We will have copies of that signal's spectrum at $f_0 + kf_s$ for all integer values of $k$. Assuming we met Nyquist, there's no problem here. The problem now comes when we have the sinusoid
$$ s(t) = \sin(2{\pi}(f_0 + f_s)t)$$
$$\mathcal{F}(s(t)) = \frac{1}{2i}\big( \delta(f - [f_0 + f_s]) + \delta(f + [f_0 + f_s]) \big)$$
When sampled, we have the spectrum copies centered at $(f_0 + f_s) + kf_s$ for all integer values of $k$. If we choose $k = -1$ then we have a copy at
$$ (f_0 + f_s) + kf_s|_{k = -1} = f_0$$
So if we had a combination of $x(t)$ and $s(t)$ and sampled it, a spectrum copy of $s(t)$ would step on $x(t)$'s copy at $f_0$. This is aliasing.


4th order high-pass filter on a DSP: standard or biquads? I have two C implementations of 4th order high-pass filter (fs = 16 kHz, cut-off=100 Hz, designed using Matlab butter() function):

1st. Done as standard difference equation implementation.
2nd. Done as product of two Biquad sections.

The final implementation will be a fixed-point DSP code.

Why is biquad implementation recommended?

The two solutions in a floating point implementation are assumed to be identical, with the two BiQuads being a factored version of the standard difference equation. The BiQuad is the better way to go for fixed point as you isolate two 2nd order systems and in doing so will be easier to keep stable under variations due to the quantization involved.

For more details on this, see the responses to this post:

How does cascading biquad sections for higher order filters work?

and

https://en.wikipedia.org/wiki/Digital_biquad_filter

Note too for those less familiar: cascading two filters is the same mathematically as convolving their coefficients; multiplying two polynomials is done by convolving their coefficients. The transfer function for filters is described by polynomials so it is all the same thing; a 4th order filter can be factored into two 2nd order filters.


Using FFT to demodulated a sine wave envelope with a cosine carrier I am trying to implement AM demodulation in Matlab by taking advantage of the fact that in the frequency domain, the spectrum of the modulated signal is shifted by the frequency of the carrier.

I wrote code in Matlab that shifts the componenets back and does an inverse Fourier Transform to reconstruct the signal. The code seems to work well enough for a cosine with a cosine envelope, but if I try using a cosine to modulate a sine, I do not get the correct signal, and my Fourier specturm seems to have  both real and imaginary components, which has me buffled: If the Fourier transform should be the spectrum of the envelope just shifted, and the envelope is a sine, shouldn't the transform be pure imaginary?

I attach my code along with a graph that shows my resultsAny help is appericiated in getting to the bottom of this.

clear all

f = 50; % Message Frequency
fc = 1000; % Carrier Frequency
Fs = 5000; % Sample Frequency
t = (0:(1/Fs):1);

L = length(t);

s = sin(2*pi*f*t); % Message Signal

mod = s.*cos(2*pi*fc*t); % Modulated Signal

% Fourier Transform
a = fftshift(fft(mod));
w = Fs*(-(L/2):(L/2-1))/Fs; % Frequency axis

%% Perform complex demodulation:

% Zero all components less than the carrier frequency
a_zeroes = a;
a_zeroes(w<fc) = 0;

% Move the spectrum back by the carrier circularly to effectivly shift
% back:
a_shifted = circshift(a_zeroes,[0,-fc]);
% Reflect
a_reflected = fliplr(a_shifted);

a_demoded = a_shifted;
a_demoded(1:floor(length(a_demoded)/2)) = a_reflected(1:floor(length(a_demoded)/2));

% Take inverse fft
demoded = ifft(ifftshift(((a_demoded))));

%% Plot fft and demodulated vs. modulated
figure;
subplot(3,1,1); plot(w,real(a),'r','linewidth',4); hold on; plot(w,imag(a),'g','linewidth',2); title('Original FFT');
subplot(3,1,2); plot(w,real(a_demoded),'r','linewidth',4); hold on; plot(w,imag(a_demoded),'g','linewidth',2); title('Shifted FFT');
subplot(3,1,3); plot(t,s); hold on; plot(t, (sqrt(2*pi)).*(demoded)); title('Original Vs. Demodulated');

Here you can see the original FFT. I reconstructed the envelope's spectrum, then did an inverse FFT, but the result has less amplitude and is out of phase.

Red is the real part, Green is the imaginary part.

As you can see from your "Original FFT" plot, the real part of the FFT has even symmetry and the imaginary part has odd symmetry. When you zero out the lower sidebands, shift, and re-create the negative frequencies, you're forcing the FFT to have even symmetry (look at your "shifted FFT") plot. The signal you re-created corresponds to a cosine; this is what you get in the "Original vs Demodulated" plot, where the recovered signal is 90 degrees out of phase with the original.

The more common approach to demodulating in the frequency domain is as follows:

Zero out the negative frequencies. This is usually accomplished with the Hilbert transform. A signal $s(t)+j\hat s(t)$, where $\hat s(t)$ is the Hilbert transform of $s(t)$, is called an "analytic signal" and its spectrum only exists for positive frequencies.
Shift the spectrum back to baseband by multiplying the analytic signal with a complex exponential of frequency equal to $-f_c$.

The result of this process is the "complex envelope" which, despite its name, is real in the case of DSB-SC modulation of a single sinusoid.


What to do after this last step? I am solving a question from book in which I have to use summation. It is as follows:

$$ \frac{1}{10}\sum_{n=0}^{9}  e^{-jk\omega_0n}          $$
The value of $\omega_0$ is $\frac{2\pi}{10}$.
What I have done so far is:

$$   \frac{1}{10} \left(\frac{1-e^{-9jk\frac{2\pi}{10}}.e^{jk\frac{2\pi}{10}}} {1-e^{-jk\frac{2\pi}{10}}}\right)$$

After this step I am stuck and can not figure out the next step?

This sum appears quite often in DSP.

\begin{align}
\sum_{n=0}^{N-1} \exp(-j\alpha n) &\stackrel{(a)}{=} \frac{1- \exp(-j\alpha N)}{1 - \exp(-j\alpha )}\\
&= \frac{e^{-j\alpha N/2}(e^{+j\alpha N/2}- e^{-j\alpha N/2})}{e^{-j\alpha /2}(e^{+j\alpha /2}- e^{-j\alpha /2})}\\
&\stackrel{(b)}{=} e^{-j\frac{\alpha}{2}(N-1)} \frac{\sin(N\alpha/2)}{\sin(\alpha/2)}
\end{align}
where the last step $(b)$ comes from identity $e^{+j\theta}-e^{-j\theta} = 2j\sin(\theta)$ and be careful at $(a)$ as Matt's comment.

Substitute $\alpha=k\omega_0=k\frac{2\pi}{N}$, $\sin(N\alpha/2)=\sin(k\pi)=0$ if $k\in\mathbb{Z}$ hence so is the sum.

Update for $\alpha=0$

It is trivial that $\exp(-j\alpha n) \stackrel{\alpha=0}{=} 1$ then the sum equals $N$.

It can be thought as the case of the final result after $(b)$  when $\alpha \to 0$ then
$$e^{-j\frac{\alpha}{2}(N-1)} \frac{\sin(N\alpha/2)}{\sin(\alpha/2)} \stackrel{\alpha \to 0}{\to} N$$
where the limit is evaluated by L'Hôpital's rule.

Further interesting discussion can be found in Laurent Duval's answer.


Discrete-time signals increase in frequency and period I am studying Discrete-Time Signal Processing by Alan V. Oppenheim and Ronald W. Schafer. I am confused when it says

periodic discrete signal oscillates more rapidly when its frequency is increased from $0$ to $\pi$

and

frequencies in the vicinity of $\pi$ are referred to as high frequencies

But I don't think it's monotonically increasing, for example, $\omega_0 = 999\pi/1000$ has a way longer period than $\omega_0 = \pi/2$.

You're right that the periods of the two signals are very different. The period of a sinusoidal sequence
$$x[n]=\cos(n\omega_0+\phi)\tag{1}$$
with $\omega_0=\pi/2$ equals $N=4$, whereas for $\omega_0=999\pi/1000$ it equals $N=2000$.
But this is not what "frequency" refers to here. Note that in general the sequence $x[n]$ as defined in $(1)$ isn't even periodic. It is only periodic if $\omega_0$ satisfies
$$\omega_0=\frac{2\pi p}{q},\qquad p,q\in\mathbb{Z}\tag{2}$$
Maximum frequency for a discrete-time signal is achieved if the samples alternate in sign. This is the case for $\omega_0=\pi$, and it is achieved for almost all samples if $\omega_0$ is close to $\pi$. For $\omega_0=\pi/2$ there are always two samples with the same sign before the sign changes, so the frequency is half of the frequency of a sequence for which consecutive samples have alternating signs.


Video processing to detect a logo in frames I am trying to detect a TV channel logo inside a video file, so simply given an input .mp4 video, detect if it has that logo present in all (or most) of frames or not (or maybe what percentage of frames include that logo). My aim is to classify ad and video portions (I noticed ads do not have the logo).

We have that logo in advance and the location is always fixed, so I believe the problem can be trivial to just a matching problem.

What is the best way to achieve this? Is there any tools or libraries that I can use for this?

If you know where the logo is then you can detect its presence, fairly quickly, with cross-correlation...with a little bit more work.

Detecting the presence of "things" using cross-correlation is a way to do template matching, where we are trying to find a specific pattern in a longer "duration" signal.

In this case, the logo is the pattern and the longer "duration" (or surface in the case of an image) is the whole frame.

The question now becomes, how do you recover the pattern? How do you recover the logo, so that you can then search for that particular target?

There are two ways to do this:

It is already available: In which case, you just download a copy of the logo, resize it to the approximate size you expect the logo to be in the frame and then use cross-correlation to find it.
It is not available and it has to be recovered by the frame itself: To recover the logo, use averaging. Take a large number of frames (say 30-50 seconds of video), and average their pixel values. The key idea here is that the logo does not "move" from its position. Therefore, the pixels that correspond to the logo will have a stable value while the pixels that correspond to the background will be taking "random" values from whatever the camera happens to pan to. When you average, the stable values will retain their magnitude and the random values will have considerable lower values. The reason I am telling you to average such a long duration is because the second image is showing a frame from a talk-show. In that one, the camera, the subject or its background are not expected to move too much, therefore, it is likely that if you average for a short duration you might "pick up" a bit of the background too. In that case, you might have to recover the logo from a different part of the video. The bottom line is that this process will give you the form of the logo. You might have to apply thresholding to get rid of those "background" pixels to get a proper "mask" image.

Once you have that small "logo" image, you can then pick up a random frame from your video feed, crop it around the approximate location of the logo and then take that part of the frame and cross-correlate it with the logo. If the logo IS present in the frame, you will get a much stronger cross-correlation peak at its approximate centre than if the logo is not present. You can then apply a threshold to detect this strong peak. If the peak is above that threshold then the logo is present, if not, then the logo is not present. That's the concept and different platforms will have different implementations (e.g. here is the same thing in Python).

Hope this helps.


good free software (or code) to separate multiples sources in a single wav? For a school project I'm working on instrument classification on multi-instrumental tracks. To do so, I'm following Arie Livshin's modus operandi, and I would like to get a good software to perform source separation: given a multi-instrumental short extract, getting n (with n being equal to the number of instruments) tracks with more or less one instrument per track.

The best I've found so far is the FASST toolbox for Matlab, but the result on a saxophone/guitar song was quite disappointing. Have anyone encountered something better I could use?

A great free toolbox is the Flexible Audio Source Separation. Their website contains some demo examples in Python, Matlab, and C++ code.

Another toolbox is the Toolbox for Blind Source Separation. However, this toolbox is not user friendly and you need to do some extra programming in C, Matlab or Python.


If the cosine function is periodic, why does it have a Fourier Transform? As far as I understand Fourier Transforms are for non-periodic signals and Fourier Series for periodic signals.

So why is it we can take the Fourier Transform of a cosine when it is a periodic function, assuming the above paragraph is correct?

Indeed there are two things you have to know.

First, it can be shown that the continuous-time Fourier transform can be obtained from the continuous-time Fourier series by letting the period $T$ go to infinity.

Second, formally speaking the Fourier transform integral for periodic signals do not converge, hence do not exist. The solution is a generalisation of the Fourier transform by the use of Dirac impulse functions.

The result is an interpretation that the Fourier transform of periodic functions is a sum of scaled Dirac impulses at the Fourier harmonic frequencies, and the scale being the corresponding Fourier series coefficients.


Why the error probability is given by this expression? I have seen this equation:

$$
\Pr({\epsilon})  = Q({\sqrt{2E_b/N_0}}) \,.
$$

where these definitions hold:

$$N_0/2  - the \ spectral \ density \ of \ white \ noise$$
$$E_b = PT_b - the \ energy \ per \ transmitted \ bit$$
$$P - Power $$
$$T_b - the \ time\ the \ bit\ transmitted \ in$$

Q - is (this) function, how do I derive this equation?

A whole chapter of a book can be written to answer the "Why" in this question but the
short answer is that if the signals used are an antipodal signal set (meaning
bits $0$ and $1$ are transmitted using signals $s_0(t)$ and $s_1(t)$ respectively
where $s_1(t) = -s_0(t)$ and the signals are of duration $T$ and energy $E$,
then with matched filtering, the receiver makes a decision that $0$ was transmitted
or a $1$ was transmitted according as the sample value from the matched filter
is positive or negative.  The sample value is a random variable $X$ whose
distribution is $\mathcal N\left(\pm \sqrt{E}, N_0/2\right)$ depending on whether
a $0$ or a $1$ was transmitted. If $X$ has mean $-\sqrt{E}$ when a $1$ is transmitted,
then the probability of error is
$$\begin{align*}
P_{e,1} &= P\{X > 0 \mid 1~\text{transmitted}\}\\
&= P\left\{\mathcal N\left(-\sqrt{E}, N_0/2\right) > 0\right\}\\
&= Q\left(\frac{0 - \left(-\sqrt{E}\right)}{\sqrt{N_0/2}}\right)\\
&= Q\left(\sqrt{\frac{2E}{N_0}}\right).
\end{align*}$$
Similarly, $X$ has mean $\sqrt{E}$ when a $0$ is transmitted,
and so the probability of error is
$$\begin{align*}
P_{e,0} &= P\{X < 0 \mid 0~\text{transmitted}\}\\
&= P\left\{\mathcal N\left(\sqrt{E}, N_0/2\right) < 0\right\}\\
&= \Phi\left(\frac{0 -\sqrt{E}}{\sqrt{N_0/2}}\right)\\
&= \Phi\left(-\sqrt{\frac{2E}{N_0}}\right)\\
&= Q\left(\sqrt{\frac{2E}{N_0}}\right).
\end{align*}$$
Finally, using the law of total probability, the probability of error is
$$P_e = P_{e,0}P\{0~\text{transmitted}\}
+ P_{e,1}P\{1~\text{transmitted}\} = Q\left(\sqrt{\frac{2E}{N_0}}\right)$$
regardless of the probabilities of transmitting $0$'s and $1$'s.
For more details such as why the matched filters outputs
are as stated above, see, for example, this lecture
note and
this one
of mine.


Inverse system of sinc? I'm doing some self-study for an important exam I'll have in late March and came across the following question:

So, using the convolution properties, if I want to find an identity system so that the final output is just the input, I would need one of the following, right?

h(t) * h'(t) = $\delta(t)$

H(w) * H'(w) = 1, so H'(z) = $\frac{1}{H(z)}$

My question is, if we go through option 2, then we are going to have a bunch of 1/0 since the Fourier Transform of the sinc will be a rect, so the system would not be stable.
Also, the sinc function itself is not invertible as we have many points that maps to 0 and the even symmetry of the function.
So does that means that for this particular case it is impossible to get such an identity system?
Thanks for the help!

Your $h(t)$ is an ideal low pass filter with cutoff frequency of $75$ Hertz. At the output, everything above $75$ Hertz will be gone and equal to zero. If you're thinking about the inverse system, you need to think about getting back the original input. What can you say about the original input signal at frequencies greater than $75$ Hertz? Isn't it anybody's guess as to what used to be there?


How do you reduce $H\left(e^{j\frac{\pi}{2}}\right)$ further according to a textbook solution I want to know how I could get from the first line to the second. I've been trying to figure it out for a while with no luck. Thank you in advance!
\begin{align}
H\left(e^{j0.5\pi}\right) &= \frac{1 - e^{-j0.5\pi}}{1 - 0.25e^{-j\pi}}\\
&=1.13e^{j\frac{\pi}{4}}
\end{align}

By plotting the unit circle you will find that $1-e^{-j\pi/2} = \sqrt{2}e^{j\pi/4}$.


How to decompose a signal into exponentally decaying sinusoids? I have a system of coupled oscillators. One is a pendulum acting as a tuned mass damper and the other is a wooden structure. When I loosen the screw that allows the pendulum to oscillate I get a signal of the acceleration of the system like the following:

I know for sure that the amplitude of the waves composing this signal decay exponentially, however, I would like to know what their decay rate is. If I perform FFT on the signal I can only get the frequencies of the waves composing the signal. Is there a procedure similar to FFT I can implement in Python/Matlab with which I can find the frequencies, decay constants and initial amplitude of the waves composing this signal?
I found the following library in github (https://github.com/NanoComp/harminv), which would exactly do what I want, but it only runs on Linux systems (I have a Windows PC). Do you know of any programme that I can use to do this?

For a continuous-time signal this decomposition is precisely the Laplace Transform, and for a discrete time signal, the z-transform. Both decompose a time domain waveform into decaying and growing exponentials.
This is clear from observing the formulas and seeing the underlying operation similar to generalized correlation which is a sum of products:
Here is the general form for the Unilateral Laplace Transform:
$$\mathscr{L\{x(t)\}} = X(s) = \int_0^\infty x(t)e^{-st}dt$$
Where $s= \sigma + j\omega$
Note that $e^{-st}$ represents a rotating phasor that is growing or decaying with time!
If it isn't yet clear, $e^{-st} = e^{-(\sigma + j \omega) t}= e^{-\sigma t}e^{-j\omega t}$, and thus we have a rotating phasor ($e^{-j\omega t}$) that is growing or decaying $e^{-\sigma t}$ both as a function of time. If that is confusing, it may further help to know that the general form of the complex exponential $Ke^{j\theta}$ describes a vector (phasor) of magnitude K and angle $\theta$. Rotating phasors are unique frequency tones and fundamental elements (basis functions) of the Fourier Transform. Such a phasor rotating counter-clockwise on a complex plane is a positive frequency  (the phase of the phasor is increasing linearly with time, and frequency is $d\theta/dt$), and one rotating clockwise would be a negative frequency. If we observed the real axis of any such rotating phasor, we would observe a sinusoidal waveform which is how we are typically first introduced to "frequencies". Real sinusoids consist of two such tones, a positive and negative frequency, as given by Euler's equation: $2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$. Thus complex tones of the form $e^{j\omega t}$ are even simpler subcomponents of real tones and once that is well understood, thinking of "frequencies" in terms of complex tones rather than real sinusoids greatly simplifies signal processing and leads to much greater possibilities! Complex need not be "so complex"!!

Let's call $e^{st}$ our hypothesis of what may be a component of $x(t)$ for any given value of $s$. We vary $s$ over all possible values to establish the result $X(s)$. In the Unilateral Laplace Transform, as given by the formula above, we multiply our signal $x(t)$ in time with $e^{-st}$ and integrate for all time, $t\ge 0$. Note the product is negated in the exponent: We spin exactly opposite to our hypothesis, and decay if it is growing and vice versa at the exact same rate as the mechanism to locate all these fundamental components! When a component is found, the rotation is stopped and the growth or decay is leveled out!

For any particular value (point)  $s_1$ of $e^{s_1 t}$ that may be a fundamental component of our signal $x(t)$, meaning $x(t)$ contains as part of it one such fundamental decaying, rotating complex exponential; then for that component, it's product with $e^{-s_1 t}$ will be $e^{s_1 t}e^{-s_1 t} = 1$ for all time (for t>0)!... which then when integrated to infinity as given by the formula for the Unilateral Laplace Transform would go to infinity-- but notice any other value for $s$ where $\sigma > \sigma _1$  (for $s_1 = \sigma _1 + j\omega_1$) would necessarily decay and therefore converge to a finite value. Thus we see how we get distinct points that go to infinity (poles) and a vertical line as the boundary with everything to the right in the "Region of Convergence".  Anything to the left of this vertical line (as given by the left most poles in the solution to the Laplace Transform) will not converge and therefore is an invalid region where no solution exists (and therefore our signal has no components there).
If plotted the Laplace Transform would be a surface of complex results for every combination on the 2 dimensional surface given by $(\sigma, j\omega)$. We often just depict this as poles (points going to infinity) and zeros (points equal to zero) since every other point on that surface is defined and unique by those pole and zero locations (but we may then miss understanding the full picture).
If we restrict $s$ to be only the $j\omega$ axis, then the Laplace Transform reduces to the Fourier Transform, which is a correlation with every possible rotating phasor that doesn't decay or grow (stays constant for all time). Specifically for the Unilateral Laplace Transform where $t\ge 0$, this would also involve multiplying the signal by a unit step ($u(t)$):
$$\mathscr{F\{u(t)x(t)\}} = \int_0^\infty x(t)e^{-j\omega t}dt$$
Notice how similar that equation is to the Unilateral Laplace Transform first given, and observe the surface plot of the magnitude of an example Laplace Transform below (the results of the Laplace Transform are themselves complex values, so each result on the surface has a magnitude and phase) for a low pass filter, and if we slice that 3 dimensional shape along the $j\omega$ axis we end up with the frequency response (Fourier Transform of the filter's impulse response, where $X(s)$ is the Laplace transform of the filter's impulse response), as given in the graphic below in the upper right hand corner (a low pass filter), and we see from the low pass filter that the impulse response consists of a decaying sinusoid (if the poles were on the $j\omega$ axis, this would be a constant level sinusoid, but given they are shifted to the left, that sinusoid is decaying with time.)

For discrete-time systems, the z-transform does the exact same thing, just with discrete phasors rotating in discrete time that are growing and decaying and instead of an s-plane we have a z-plane where the unit circle corresponds to the $j\omega$ axis. With discrete processing we use algorithms that compute the Discrete Fourier Transform (DFT) such as the ubiquitous FFT. The corollary for computing the Z-Transform is the chirp Z-transform (CZT), and for efficient implementations of this see Bluestein's algorithm.


Why does Phase-locking value show spurious effects at the beginning & end of a signal? I am trying to understand the behaviour of the phase-locking-value, as e.g. defined here, by some simple examples.
Basically, I am creating two signals with the same frequency & amplitude, where the phase-offset between the two is random. Now, as expected, at most time points the PLV gives me a low value indicating that there is not a consistent pattern of phase onset. However, at the beginning & end of the signal, there are some strange effects. I suppose that these might be somehow introduced by the Hilbert transform, or is this a problem with my implementation?
Here is how I defined the function:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

def plv(x,y,trial_axis=0):

    hilbert_x = signal.hilbert(x)
    hilbert_y = signal.hilbert(y)

    normed_x = hilbert_x/np.absolute(hilbert_x)
    normed_y = hilbert_y/np.absolute(hilbert_y)

    p = np.absolute(np.mean(normed_x*normed_y.conj(),axis=trial_axis))

    return p

Here we have a plot of the behaviour:
t = np.linspace(0,10,1000)
N = 50
s_1 = np.array([np.sin(10*2*np.pi*t) for i in range(N)])
s_2 = np.array([np.sin(10*(2*np.pi*t)+np.random.uniform(low=0,high=2*np.pi) ) for i in range(N)])
p = plv(s_1,s_2)
plt.plot(p)

What was previously posted as an answer should have been an extended comment -- sorry for that. There simply was not enough space for what I thought I'd write about how the Octave code seemed to behave. Now that the code is cleared, what should have been the answer in the first place (given the rather obvious signs), is now the answer: spectral leakage.
The Hilbert transform is calculated by zeroing half the FFT of the signal and then recomputing the analytical signal through the IFFT. Brute-force zeroing will always lead to some leakage due to the nature of the limited length of the data. This can be easily tested with your code. s1 is a simple sin(), so its Hilbert transform will be a -cos(). Simply plot the difference between -cos()-hx.imag, and you'll see the leakage. The same goes for the second signal.
This is the changed code in Python:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
f = 10
t = np.linspace(0, 100/f, 1000)
N = 50f = 10
nt = 1000
t = np.linspace(0, 100/f, nt)
N = 50
hor = np.ones(nt)
ver = np.ones(N)
r = np.array( [ np.random.uniform(low=0, high=2*np.pi) for i in range(N) ] )
s1 = np.array( [ np.sin( 2*np.pi*f*t ) for i in range(N) ] )
s2 = np.array( np.sin( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
c1 = -np.array( [ np.cos(2*np.pi*f*t) for i in range(N) ] )
c2 = -np.array( np.cos( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
hx = signal.hilbert(s1)
hy = signal.hilbert(s2)
nx = hx/np.absolute(hx)
ny = hy/np.absolute(hy)
p = np.absolute( np.mean( nx*ny.conj(), axis=0 ) )
plt.plot(p)
plt.show()

and this is the test (c2 will show non-overlapping plots due to the randomness):
plt.plot(c1.transpose() - hx.imag.transpose())
plt.show()

The Octave code, for checking:
f = 10
t = linspace(0, 100/f, 1000);
N = 50;
c = ones(N, 1);
ft = f*t.*c;
s1 = sin(2*pi*ft);
c1 = -cos(2*pi*ft);
r = rand(N,1).*c;
s2 = sin(2*pi*(ft + r));
c2 = -cos(2*pi*(ft + r));
hx = hilbert(s1');
hy = hilbert(s2');
nx = hx./abs(hx);
ny = hy./abs(hy);
p = abs(mean((nx.*conj(ny))'));
grid on
subplot(2, 1, 1)
plot(p)
subplot(2, 1, 2)
plot(c1' - imag(hx))


Simple Filter representation in Matlab Suppose I have an IIR filter in the $z$-domain in the following form:
$$
H\left(z\right)=\frac{1}{1-0.2z^{-1}-0.1z^{-2}}
$$
How do I represent this in MATLAB?
I am pretty sure if I just listed the coefficients as H = [1 -0.2 -0.1];
this would be considered wrong. Can somebody help me with this please?

You have to define numerator and denominator polynomials. In your case you have

b = 1;
a = [1,-0.2, -0.1];


Lightweight audio listener for detecting blip The challenge I'm looking at is to design an iPhone component that listens to microphone input and triggers when it hears a particular blip (preferably < 0.1sec duration).

I can design the blip myself, and also false positives are okay.

The key constraint is that the listener is lightweight, so it doesn't drain CPU.

Can anyone suggest a mechanism, or point me in a good direction?

π

EDIT:  this is as far as I have got:

I should be feeding audio into a ~1s buffer

so the first check may throw up lots of false positives,  but this allows a second check to delve into the one second history and determine whether it was indeed a false positive or a valid blip.

so the first check may just be listening for the presence of a particular frequency sustained over 0.1s

to do this I could feed incoming audio into a 100 sample (for example) ring buffer,  and then, 50 times per second multiply the ring-data by a sine wave (giving x) and cosine wave (giving y), effectively manually calculating a single FFT bin x+iy.  then I could flag whenever the phase remains roughly constant over three frames.

or maybe I could use 3 such rings,  length 99, 100, 101 samples.   and monitor the amount of energy getting accumulated in each;  if the 100-sample-ring wins and the neighbours are roughly balanced, that would be an indicator that the 100-ring's resonant frequency has been detected.

but this approach would require feeding all audio data into three rings. it isn't  particularly lightweight...  maybe I could only feed one in 17 (say) samples into each ring, throwing away the other 16.   As 17 is co-prime with 99,100,101  it would fill all cells of the ring before overwriting the first again.

One lighter-weight approach might be to use a narrow-band IIR filter centered at your blip frequency, rectify or square, feed this result to a low-pass filter, and look for a value passing some threshold.  A Goertzel filter is another possibility.

Note that the most current iOS Accelerate framework includes a reportedly power-efficient biquad IIR filtering function.  Also note that on recent iOS devices, even a more heavyweight approach (FFTs, autocorrelation, and such) on a single channel of real-time audio input might use only a tiny percentage of total device power.


Negative frequency of a sinusoidal wave on the magnitude spectrum: Physical explanation? link to the video: Negative frequency

Must a cosine function be expressed as two exp functions with opposite frequencies/directions of rotation and half the magnitude/amplitude(1/2),
therefore resulting in two magnitude spectrum components?

Why can't there be one single component at w?

I know this is very basic but there are only videos about the "how" not the "why".Thank you in advance!

Regarding your first comment, I have a figure that might help, in terms of intuition. It's a bit badly scaled but it might do the work.


what symmetry does system function H(z) have if h[n] is real? assuming h[n] is real...

If frequency response $H(e^{j\omega})$ is Conjugate Symmetric:

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

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

Then, what symmetry does the System Function $H(z)$ have?  Is it the same symmetry as the frequency response  $H(e^{j\omega})$, or is some type of symmetry like conjugate reciprocal symmetry?

$$
H(z) = H^*\left(\frac{1}{z^*}\right)
$$

This is confusing the heck out of me, because they start with something that is conjugate symmetric as a frequency response then convert it to a z transform by substituting $z=e^{j\omega}$, then, suddenly by changing variable to z its "conjugate reciprocal symmetric", instead of just plain "conjugate symmetric"... except it doesn't actually tell you this anywhere... that the system function has this symmetry...

Would this be some weird property that "Conjugate reciprocal symmetry" is the same as "conjugate symmetry" on the unit circle. and the frequency response is really a degenerate case of "conjugate reciprocal symmetry"

it's because $e^{j\omega}$ is such a signal that has the property:

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

then going back from the DTFT to Z-transform is not obvious. So it's better to start from definition of Z-transform to define the symmetry:

$$ X(z) = \sum_{n} x[n] z^{-n} $$

and apply this to a signal which is real; $x[n] = x^*[n]$

$$ x[n] \xrightarrow{ Z-transform} X(z) \implies x^*[n] \xrightarrow{ Z-transform} X^*(z^*) $$

and conclude for real $x[n]$ :
$$ X(z)  = X^*(z^*) $$

Derivation:

$$\begin{align}
\mathcal{Z} \{x^*[n] \} &= \sum_n x^*[n] z^{-n} \\ \\
& = (\sum_n x[n] (z^*)^{-n})^* \\ \\
& = X^*(z^*) \\
\end{align}
$$

As you can see, on the unit circle where $z = e^{j \omega}$, this becomes:

$$ X(e^{j\omega}) = X(e^{-j\omega})^* $$
the conjugate-symmetry property...


What is $F(0)$ is "dc" component in the context of image processing? It has always been said that $F(0)$ is the "DC component" in fourier transform. However, I don't get what it means to say that $F(0)$ is "DC" in the context of image processing.

The zero in this case just meant zero frequency, and hence no changes isn't it? Then in the context of image processing, how should I imagine the picture to be like at $F(0)$? What exactly is that "DC" component?

The zero frequency of an image DCT is the mean gray value of the pixels of the input image (for graylevel images).
$F(0)$ is not an image: it is a single coefficient.


Continuous-time frequency response question Suppose the CTFT of continuous-time input $x_c(t)$ to an LTI system is $X_c(j\Omega)$ and that of its continuous-time output $y_c(t)$ is $Y_c(j\Omega)$. We have,
$$X_c(j\Omega) = 0,\phantom{1}\text{for } |\Omega| \ge \frac{\pi}{T_1}$$.

The LTI system samples the signal (C/D operation at sampling rate $1/T_1$) and passes it through the discrete-time system with response $H(e^{j\omega})$. The result of this DT system is then passed through an ideal D/C system (operating at rate $\frac{1}{T_2} = \frac{1}{2T_1}$) to yield $Y_c(j\Omega)$.

I want to find out the overall frequency response of the system. I can write

$$Y_c(j\Omega) = \begin{cases}\frac{T_2}{T_1} H(e^{j\Omega T_2})X_c\left(j\Omega \frac{T_2}{T_1}\right),\phantom{1}\text{for } |\Omega| < \frac{\pi}{T_2}\\0,\phantom{1}\text{for } |\Omega| \ge \frac{\pi}{T_2}\end{cases}$$

Since I get $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ in the above expression and not $X_c\left(j\Omega \right)$, I am not sure how to find the overall frequency response. Is it possible to express $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ in terms of $X_c\left(j\Omega \right)$ using some continuous-time filter?

Is it possible to express $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ in terms of $X_c\left(j\Omega \right)$ using some continuous-time filter?

an LTI system cannot do that.  continuous-time "filters" need not be LTI but usually we mean LTI when we consider "filters" in analog or digital filters.

you can change the spectrum of $X_c\left(j\Omega \right)$ to $X_c\left(j\Omega \frac{T_2}{T_1}\right)$ by speeding up (or slowing down) the time index of $x_c(t)$ to $x_c\left(t \frac{T_1}{T_2}\right)$.

it appears you want a pitch shifter or similar.  is that it?


Removing the overshoot in step response First of all I'm new to DSP so excuse my simplified words & I'm sorry if it's a duplicated question.

Now, I'm trying to filter a signal but as you all know, the filter step response overshoot the filtered signal at the beginning.

Here it's the signal. It begins with zeros then the signal comes up

Is there any accurate solution to remove this overshoot completely?

** And if there is references or any direct command in Matlab, I would be grateful to know it.**

Thanks in advance..

You said you were new to DSP, but in reality, anyone new to filtering (analog or digital) has the same problem. In other words, we all learn about filtering in the frequency domain but pay little attention to what is taught about the effects of filtering in the time domain.

The problem you are showing, excessive overshoot and ringing, is most easily solved by choosing a diffrent type of filter. You don't say say what type of filter you are using, but it most certainly has a very sharp transition band, and a flat pass band, or a ripple in the pass band, which allows for a sharper transition band.

An extreme example of just the oppisite of your filter is the Gaussian response, which has roll off in the pass band and a wide transition band. Its time domain response will have no overshoot or ringing, but also has "poor" attenuation characteristics.

The solution then, is to design a filter that gives as much filtering as needed, but doesn't cause any more overshoot and ringing than your application can tolerate. You need to learn that with filters, you can have sharp filters, or you can have minimal distortion (overshoot and ringing), but not both.

I am not a MatLab user, so I can't suggest a specific MatLab function to use. Might I suggest some free FIR and IIR software that will clearly and easily demonstrate what I have said here. Both the FIR and IIR programs allow you to adjust the filter's response in the freqency domain, and quickly see the effects in the time domain. In the IIR program, experiment with the Gauss and Elliptic responses as extreme oppisites of one another. The programs are available at  http://www.iowahills.com/

Also, you need to pay attention to the filter's group delay. For 2 filters with essentially the same magnitude response, the filter with the more linear phase (flatter group delay) will cause less distortion.


How does the power spectrum of the amplitude envelope depend on its underlying carrier frequency? when looking at the amplitude envelope of an oscillation (see red trace in Figure below), it seems to me as if the frequency content of the envelope depends on the the frequency of the underlying carrier signal (blue trace). More specifically, it looks as if the fastest fluctuation the envelope can exhibit is related to the frequency of the carrier signal. There should be a straightforward mathematical relation between the carrier frequency and the fastest frequency in the envelope. If I think about the individual peaks of the oscillation as "samples" for the envelope, then the fastest frequency the envelope contains should be carrier_freq/2, but I'm not sure if this is correct. I'd appreciate any help or hints in the right direction.

Your intution is right, the carrier frequency has to be at least twice the highest frequency in your signal. Practically, much higher carrier frequencies are chosen, as you would need ideal or nearly ideal lowpass filters for carrier frequencies this low.


Complex analysis book recommendation for electrical. eng I need recommendation about complex analysis book. As I'm electrical engineering student, it should cover everything one engineer need to know about that mathematical field, but without strict mathematical formalism.

Saff and Snider's Fundamentals of complex analysis is a nice enough book. It's just excessively pricey.

Brown & Churchill's Complex Variables and Applications is a good option as well (you can get an old edition for pretty cheap, and its content-wise the same as the new ones).

Marsden & Hoffman's Basic Complex Analysis is also a good option, but its a tiny bit more advanced than the other options (but still suitable for Eng/Math/Physics students). A nice feature of this is that he has a lot of hints and worked out problems.

The OCW courses linked in Matt L.'s answer are alright -- the second one is based around Ahlfors' book which is more advanced than what most EE's will need (and just has some extra stuff to supplement Ahlfors).

If I were to pick one, I'd go with Brown & Churchill as its the one I used as an undergraduate, with Marsden & Hoffman as a close second.


What's the point of defining the signal over the whole time domain? This question is classic for anyone starting with some signal processing course, suppose $y(t)=x(t/2)$ then the system is noncausal because we have that the output at t=-6 depends on the input at t=-3 which is in the future. Suppose the signal is only periodic for positive t, then the whole signal is aperiodic.
But how does it even make sense to write a negative value for time, shouldn't any "real world" signal start at some instant, which we would call t=0, does negative time do anything other than make us refrain from writing too many piece-wise defined functions?

Bottom line: it makes the math easy.  In addition, there are certainly instances where we care that a signal has existed for some time in the past.

But how does it even make sense to write a negative value for time,

Think of $t=0$ as "now", $t < 0$ as "before now", and $t > 0$ as "later".  Does that make sense?
By extension, $t = 0$ is just any arbitrary moment in time; if you don't like it, shift everything over.

shouldn't any "real world" signal start at some instant

Aside from the big bang, no, why should there be some instant where a signal "starts"?
Note that the concept of signals that exist for all time comes from Fourier analysis.  If you have some reason to start considering signals at an instant, there is Laplace analysis, where time pretty much starts at $t = 0$.

does negative time do anything other than make us refrain from writing too many piece-wise defined functions?

Fourier analysis is a good way to analyze signals that are "just there", like radio signals or speech.  If your analysis revolves around a system that's in steady state, processing a signal that has had pretty much the same statistics that it has now and will have in the future, then treating that signal as existing for all time is both sensible, and makes the math easy.


Positive and negative frequencies: how to prove they are independent? This seems to be quite a basic question but is there some detailed proof why interference power (i.e., power spectral density) coming from the positive and negative frequencies are independent? I am considering two real baseband OFDM signals which are interfering with each other and where $2N-f=-f$ is relation between positive and negative frequencies while $2N$ is total number of frequencies (tones). Assume that two signals are not sync and that is the source of interference leakage coming from the negative spectrum. To simplify, I want to prove why is the received PSD equal to : $PSD_{rec.,k}=A^2_{pos.,k}+A^2_{neg.,k}$ instead of $PSD_{rec.,k}=A^2_{pos.,k}+2A_{pos.,k}A_{neg.,k}+A^2_{neg.,k}$ where $A^2$ is squared amplitude (i.e.,DFT coefficient) of tone $k$. Could anyone recommend me some literature or give me a short answer on this question?

This is a bit waffly, but gives the general idea.

The usual reason for this is because of expectations used when we talk about Power Spectral Densities:
$$
S_{xx}(\omega) = \int_{-\infty}^{+\infty} \gamma(\tau) e^{-i\omega \tau} d\tau
= \int_{-\infty}^{+\infty} E\left[ x(t) x(t+\tau) \right] e^{-i\omega \tau} d\tau
$$

So what you have to look at is how the $E[\cdot]$ works on the cross-products between the positive and negative frequencies. Assuming $x$ is just:
$$
x(t) = \sum_{k=1}^{N_+} a_k e^{j\omega_{+,k} t} +  \sum_{k=1}^{N_-} b_k e^{j\omega_{-,k}t}
$$
then each of the cross terms will look like:
$$
\lim_{T\rightarrow \infty} \frac{1}{T} \int_{-T}^{+T}e^{j\omega_+t}e^{-j\omega_-t}dt = \delta(\omega_+ - \omega_-)
$$
which will be zero if $\omega_+ \not= \omega_-$. That makes all the cross-terms zero, so you're just left with the $x_+(t)x_+(t+\tau)$ and $x_-(t)x_-(t+\tau)$ terms, i.e. your $A^2_{pos,k}$ and $A^2_{neg,k}$ terms.


Time-bandwidth product and uncertainty on harmonics For a signal that is not corrupted with noise, how is the time-bandwidth product affecting the uncertainty or error introduced when extracting the fundamental frequency ? I mean, instead of choosing the fundamental frequency of the signal at $f_0$, you choose a harmonic $f_k$ at $k\cdot f_0$ and divide by $k$ to get $f_0$. How is the uncertainty related to $k$ for a given frequency resolution ?

If you a-priori assume that $k \cdot f$ is a stationary harmonic of periodic waveform of frequency $f$, that implies that you don't have any separable time locality of that harmonic.  But any tighter bandwidth of a frequency estimation might come from the a-priori assumptions regarding the stationarity of the harmonics, not from the data.


Matlab and R functions for fft(): N-point FFT and dimension issue In MATLAB, Y = fft(X,n,dim) returns the Fourier transform along the dimension dim. For example, if X is a matrix, then fft(X,n,2) returns the n-point Fourier transform of each row. However, in R is different. Do you know any function to solve this; so that takes into account n and dim, to get a n length complex vector? See that in the example that n>ncol(X).
MATLAB code
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector

x1 = cos(2*pi*50*t); % First row wave
x2 = cos(2*pi*150*t); % Second row wave
x3 = cos(2*pi*300*t);  % Third row wave

n = 2^nextpow2(L);
dim = 2;
X = [x1; x2; x3];

Y_right = fft(X,n,dim); %Result: 3*1024 complex double

 Y_wrong = fft(X); %Result: 3*1000 complex double (I am looking for the above result)

R code
Fs <- 1000
T <- 1/Fs
L <- 1000
t <- (0:(L-1))*T

x1 <- cos(2*pi*50*t)
x2 <- cos(2*pi*150*t)
x3 <- cos(2*pi*300*t)

n<-2^ceiling(log2(abs(L)))
dim <- 2
X <- t(matrix(c(x1, x2, x3), ncol=3))

Y <- fft(X)
# Result: 3*1000 cplx and it is not the same as `Y_wrong` (operates along each row)
#To get Y_wrong:

Y <- matrix(0, nrow=3, ncol=1000)

for (i in 1:ncol(X)) {
Y[,i] <- fft(X[,i], inverse = T)
}

My aim is to get Y_right, and Y think that by using apply we get the same as Y_wrong. Thank you in advance.

R allows you to loop over rows or columns through the use of the apply family of functions. Note that apply is not necessarily faster than looping but more in line with the functional programming nature of R. There are other apply functions such as lapply, vapply, sapply etc but for your purposes, apply should suffice.  Also, note that you may have to transpose the result using t(result). 1 means go over rows and 2 means go over columns. So, if X is the matrix,

apply(X, 1, fft) does the fft of the rows.

apply(X, 2, fft) does the fft of the columns.

Note that the use of .row or .col is my convention. I could have used x, y , or z because that's just a dummy variable from apply's point of view.

Also, you can type ? apply in R in order to obtain the help for apply. Sometimes they can be kind of obscure but atleast they tell you what arguments are available.


How to set the "duration" parameter for a flanger with LFO oscillator I am reading this book Designing Audio Effect Plugins in C++: For AAX, AU, and VST3 with DSP Theory. Here I implemented the Flanger effect. Quick summary: A flanger effect is basically a echo effect but the delay parameter is changed with an LFO oscillator to give the "flanging" sound.
Here is the graph:

My problem is this version of flanger only has the following properties:

In modern DJ equipments you see flanger with options such as "1beats,2beats,4beats...". What I noticed is The smaller lfoRate is the longer my flanger effect becomes. So If I set lfoRate around 10.26, my flanger becomes 1/8beats flanger. The book does not mention anything related to this. So I wonder what is the relationship between lforate and "Xbeats" of flanger. Lets say that my music is 128 bpm, is there is a special formula which is used to calculate lfoRate so that I have 4 beats flanger?

The LFO rate is the speed of the delay modulation. One full cycle is going from "down" to up" back to "down" again. It's in Hz, i.e. "per second". Let's says the music is 120 bpm and your LFO rate is 0.5Hz. 120 pbm is 2 beats per second or 4 beats per flanger cycle.
Beats Per Cycle = Music_speed/60/flanger_rate


How important is it to use power of 2 when using FFT? Here is the problem. I have a 2D array of data, first column represents the time data and the second column represents the sinusoidal response data, based on the time data. I apply fft and I get my frequency (that I started with) in a specific bin as I expected and I find the amplitude and phase angle from that bin. Now the problem is I have the same set up but with more data points, I apply the fft again and the bin number changes (which is normal and it is where I expect it to be), the amplitude is the same but the phase angle is different) first is this normal? second, what approach should I take? Thank you

PS: neither of the set ups (mentioned above) give data of length of power of 2, say the first one gives 1620 data points and the second one gives 1745 data points, so should be taking the next power of 2 for both from the beginning?

Modern FFT libraries, such as FFTW and Apple's Accelerate framework can do non-power-of-2 FFTs very efficiently, as long as all the prime divisors of the composite length are fairly small (2,3,5,etc.)

A power of 2 makes it simpler (about 1 page of source code) if you have to code your own FFT for some reason, or are otherwise constrained as to max program length (or FPGA gates, etc.)

For phase measurement, it might be easier to do an fftshift (pre-rotate the data by N/2) to reference FFT phase to the center of the data window, where the evenness/oddness ratio, and thus the phase won't change or alternate with bin number (for phase that is the same at the center of that data window) even for signals that are non-periodic in the FFT length, as you vary the length.


HHT : Toning down noisy data I've recently asked a question on stack overflow asking for different topics one can visit to solve the problem of wave re-construction. One of the recommended solutions was to implement the Hilbert-Huang transform, as it can help in understanding non-linear/natural data.

Here are some samples and details surrounding their capture to help me state my question.

I am collecting pressure readings from a pressure sensor at 30Hz. When I put pressure on the sensor, the Y axis above expands. When I release the pressure, it slowly approaches towards the idle state. In this case, the idle state has a base pressure reading of roughly -30.

Samples ranging from 0 to 200, and even those ranging from 400 to 750, have a normal pattern of pressure readings. They are very natural, in that they are not always the same when increasing and decreasing. They are also very natural in that their pressure idles around -30.

The same cannot be said about the samples ranging from roughly 270 to 370. The reason that the samples looks significantly different than the remainder of the sampled data, is because I bend the sensor as I am applying and releasing pressure. Instead of idling around -30, my sample readings drop down to as low as -140. Note: I have highlighted the area of impact in red in the below image.

Ideally, what I'd like to play around and perform is re-construct that wave to be in line alike the rest of the data. The green line in the below graph draws out what i expect to see as a final result of re-constructing the wave. ( It's just a rough sketch to explain the idea )

With that in mind, onto what my question is about.

Below you will find the generated HHT results of my sample data, with MATLAB.

Interestingly enough, as soon as IMF of the 2nd and 3rd round come about, the "abnormal" noise is very evident between samples of 300 to 350ish. Ideally, what I'd like to do is minimize this tone to be in line with the rest of the data.

Questions:

While my noob HHT understanding is poor, would it be logical to set the stop criteria ( in this case ) to IMF of C2 or C3, and use the results to somehow subtract out the IMF from the original wave to re-construct the wave as i desire?

Is there maybe a better way I can expand my utilization of HHT to solve my problem? ( Any topics or brainstorming ideas would help )

UPDATE after Cye's reccoemndation
I managed to filter out CF2 and CF3, however, the results were not as I expected. ( Image below ) The graph only became smoother, but that large dip was still present.

As a pragmatic applied mathematician who has job to do, say, rather than to prove something via theorem or induction, my first approach would be to reconstruct the time series without IMF C2 and C3 and see what happens. I think that's the equivalent of using noise filters in the FFT or deleting outliers in a data set.

As an aside, I'm impressed that you could pick up and utilize the HHT so quickly. Kudos to you.


What is zero-phase filtering and forward-backward filtering? I am pretty familiar with IIR and FIR filtering and I have implemented them in my several projects. However, recently i found something called zero phase filtering. I have tried to understand it but not going well. All i know is that zero phase filtering is filtering that do not produce any phase delay and distortion

Is it right?

I have understood that FIR have a constant delay for a group signal while IIR have a non linear delay based on the frequency, so it will distort the signal phase.so there is a zero phase IIR to overcome the cons

How does the zero phase IIR work? (Needs an overview and maybe some examples)

And i have read about forward backward filtering that also remove the delay and distortion of filtering

How does it work?
Is it one of the zero phase filtering method?
It said that it can be done in offline filtering. What is offline filtering?

Thank you

Zero-phase filtering is a non-causal procedure, so it cannot be done in real time, only offline (or pseudo real-time, i.e., with a sufficient delay). A zero-phase filter needs to have a purely real-valued frequency response, and, consequently, it must have an impulse response that is even with respect to the time index $n=0$, i.e., it is non-causal.

Zero phase filtering with IIR filters is achieved with forward-backward filtering, as implemented in Matlab's filtfilt function. The resulting total frequency response is the squared magnitude of the original IIR filter's frequency response. Since the squared magnitude is real-valued, the resulting filter is a zero-phase filter. Of course there's a delay because you have to feed the time-reversed output of the first filter pass back to the input of the filter, after which you need to time-reverse the output. More details about forward-backward filtering can be found in this answer.


Recover curves from noisy collection of points Background: I'm trying to make a system that tracks a number of bubbles in a video

I'm implementing the bubble detection in the single image case using the Circular Hough Transform. Due to occlusion, blur and other factors, this detection will never be 100% accurate. I am tuning the detection procedure for high recall, possibly at the expense of precision.

Once this is done and applied to a sequence of frames from a video, I will have a number of detections that can be characterised as points in 4D space - x-position, y-position, radius and frame index.

Is there a procedure that can recover curves from this 4D point cloud?

Your problem is very similar to the cell-tracking problem. This has been solved quite well using tracklet based approaches, where each object (bubble, cell etc) is modeled as a Maximum-a-Posteriori (MAP) estimation. The solution to which is through linear programming.
You can get a good idea of the approaches and a web-based implementation here


spectrogram of a leak in a water pipe I'm doing a work on the acoustic emission of a leak. I am using a 1.80 m cylindrical tube with water, and the leak is simulated by means of a tap that opens after I have applied negative pressure inside the tube.
I am taking this spectrogram as a test on a tube filled with water and simulating a leak. I notice that in the spectrogram there are some horizontal lines visible at some specific frequencies throughout the spectrogram, and I have no idea how to interpret those lines at some specific frequencies.
The sharp difference in the middle of the image is when the tap is turned on.

I don't know, but my guess on this is: just like air flowing through a flute, water flowing through a tube causes vibrations, tones. The strong tone with several overtones is pretty consistent with this being similar to an air instrument.


IIR Hilbert Transformer I'm beginning to explore discrete Hilbert transformers - ways to achieve 90°. phase shift across a band of perhaps 6 kHz at a 44.1 kHz sampling rate. I'm trying to stick with IIR filters in order to save computation time - in addition to the fact that I only need to phase shift just a band of frequencies.

I saw this post where the author Ross Wilkinson alludes to creating an IIR based Hilbert transformer by setting the transfer function's zeros to:

\begin{equation}
{(q_n)}^{k}_{n=0},\quad q_n = e^{\frac{\pi}{2^n}}\quad
\end{equation}

The RHS poles inside the unit circle (0 < x < 1) are the inverse of the zeros:
\begin{equation}
{(p_n)}^{k}_{n=0},\quad p_n = \frac{1}{q_n}
\end{equation}

It would be great if somebody could point me towards another resource that would help explain where that comes from in greater detail.

FYI, I'm trying to make a single side band modulator for audio frequencies. It should be able to process in real time on a laptop, with a sampling rate of 44.1 kHz.

This is achievable with two parallel all pass filters.
The two all pass filters synthesize an odd ordered low pass filter whose pass band extends from -90º to +90º in the z-domain. (I will discuss this below).
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
The low pass filter is then rotated by +90º so that its pass band extends from 0º to 180º, which approximates the Hilbert transform. Rotation mathematically is:
$$
H_{Hilbert}(z) =G_{lowpass}(-jz)
$$
As a consequence, one of the all pass filters becomes completely imaginary - this is the Hilbert transformed signal path. By necessity, since the filter is approximating this Hilbert transform, then the imaginary branch must be applying a 90º phase shift relative to the output of $A_0(z)$.
$$
H_{Hillbert}(z)=G_{lowpass}(-jz)=\frac{A_0(z)+jz^{-1}A_1(z)}{2}
$$

Example with Elliptical Filter
The following example was taken from Handbook for Digital Signal Processing pg. 920
You can also use this program that I interpreted to design a new half-band elliptic filter with a higher order.
Starting with a real, half-band elliptical low pass filter $G(z)$ whose frequency response satisfies
$$
1-\delta_1\leq|G(e^{j\omega})|\leq1\>\>\>\>\>for\>\>0\leq\omega\leq\omega_p
$$
$$
\>\>\>|G(e^{j\omega})|\leq\delta_2\>\>\>\>\>for\>\>\omega_s\leq\omega\leq\pi
$$
To obtain half band we choose:
$$
\omega_s+\omega_p=\pi
$$
$$
(1-\delta_1)^{2}+\delta_2^{2}=1
$$
But I just used the examples given seventh order poles or an elliptical low pass half-band filter:
$$
z=0,\>\>\>z=\pm j0.436688,\>\>\>z=\pm j0.743707,\>\>\>z=\pm j0.927758
$$
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
Distribute the poles between the two all pass filters like so:

$$
A_0(z)=\frac{z^{-2}+0.190696}{1+0.190696z^{-2}}\cdot\frac{z^{-2}+0.860735}{1+0.860735z^{-2}}
$$
$$
A_1(z)=\frac{z^{-2}+0.553100}{1+0.553100^{-2}}
$$
Rotate by 90º
$$
H_{Hilbert}(z)=G_{lowpass}(-jz)=\frac{A_0(-jz)+jz^{-1}A_1(-jz)}{2}
$$
$$
A_0(-jz)=\frac{0.190696-z^{-2}}{1-0.190696z^{-2}}\cdot\frac{0.860735-z^{-2}}{1-0.860735z^{-2}}
$$
$$
A_1(z)=\frac{0.553100-z^{-2}}{1-0.553100^{-2}}
$$
Phases of $A_0(z)$ and $A_1(z)$:

Magnitude Response of $H_{Hilbert}(z)$:

Phase Difference Between $A_0(z)$ and $A_1(z)$ Branches:


Best method to extract phase shift between 2 sinosoids, from data provided I have been asking around if the way I was extracting phase shift (lag) was correct, and I ran into some trouble.

So in general, given 2 arrays of data of the same length, representing:

the input sinusoid and
the response, also a sinusoid.

So not knowing the functions that produced the input (but knowing the frequency, starting phase of the input, fs, and so on), what is the best way to find the phase lag?

should I do it "manually" by looking at the graphs and finding the phase lag at each component?
should I use FFT, like I have been trying?

any suggestions, links, books are welcome.

You can use the cross-correlation function to determine the lag between the two signals.


Interpretation of a sampled signal in the frequency domain Suppose we have a signal $y(t)$. For sampling we use a Dirac $\delta$ function.
The sampled function is $\widetilde{y}(t)$.
So $$\widetilde{y}(t) = y(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT).$$ The Fourier transformation is $$\widetilde{Y}(\omega)=\mathcal{F\{\widetilde{y}(t)\}} = \sum_{n=-\infty}^{\infty}y(nT)e^{-i2\pi\omega nT}$$
where $i$ is the imaginary unit and $\omega$ angular frequency. $y(nT)$ are the sampled points but I am having problems in understanding what $e^{-i2\pi\omega nT}$ of $\widetilde{Y}(\omega)$ means? I could write it down with $\cos$ and $\sin$ but this does not make it clearer. I would appreciate it if someone could explain it to me.

First, $i$ is not the imaginary part, but the imaginary unit. This tells you that the signal, under suitable conditions, admits a discrete time Fourier transform.  Its is a linear representation of the discrete signal. It provides you with a spectral interpretation, ie a view of the content of the signal  as a linear combination of waves, pure frequencies  indexed by a continuous index $\omega$. Pure frequencies are cisoids or complex exponentials/sinusoids, by definition

$$e^{-i2\pi\omega nT} = \cos(-2\pi\omega nT)+i\sin(-2\pi\omega nT)\,.$$

The relative importance (weight) and "timing" of these waves is given by $\tilde{Y}(\omega)$, by its magnitude and phase, respectively. Another interpretation of the sum of products is a scalar product between the signal and the different waves. If the scalar product is high, the signal looks like the pure wave  a lot, and the corresponding frequency is very "present" in the signal. If is is low, this component is not really present in the signal. The whole spectrum, made of all the $\tilde{Y}(\omega)$, tells you about the content in "each" frequency in your data.

The two links above drive you to the excellent site of J. O. Smith, and the mathematics of the discrete Fourier transform, an perfect starting point to understand these concepts.


Do filters work on a "per sample basis"? I'm confused because of filter length, whether (such) filters can be used to filter audio on a "per sample basis"?

By per sample basis I mean that I would like to filter audio one sample at a time, but vary the filter parameters even one sample at a time.

What then confuses me is, can the filter have any effect if only one sample is input to it? Or whether the filter needs a longer sample (e.g. that matches the filter length)?

Since the filter type seems crucial, I'm specifically interested in the Parks-Mccellan FIR and the context is dynamic equalization (in which it seems like the parameters such as gain would need to be varied on a close to per sample basis). I.e. I'm interested in whether the PM algorithm can be utilized for dynamic equalization in such way that one recomputes the PM between some samples and then filters those samples, i.e. in a block-like fashion.

Obviously, filters need more than one sample of input – I mean, a single sample is a single number, and how would a single number have something like a frequency? Filtering is something you apply to a digital signal, and signal is defined by being a changing entity – i.e. different samples.

Hence, digital filters always need sequences of samples. The typical (non-rate-changing) filter works in a matter that each time you "push in" a new sample, an output sample "pushes out" of the filter.

You should definitely brush up your knowledge on what a digital filter is, and especially how FIR filters work. The important point here is the understanding that application of a filter to a signal is a discrete convolution.

You're absolutely free to change the system with which you convolve any time – but you'd break the time-invariance that makes LTI (linear, time invariant) systems such as FIRs so "easy" to handle.

EDIT: with your clarification:

Yes, any FIR will need "history" of samples. But also, yes, you can recompute the coefficients anytime, and "switch over" to a different FIR. However, you must make sure your new FIR's "delay"/memory elements were already filled with the previous input samples, otherwise you'll get inconsistencies.