Why can adding delay improve the phase fit in fitting complex transfer functions? When you fit an IIR filter to a complex transfer function you can use a delay to get a stable filter and improve the fitting results in your phase response. Can anyone explain me the reason for this and/or link me to some literature explaining this? I can't find any good literature on this.

Even though an IIR filter of a given order can in principle approximate any given phase response, you will get a better fit if you choose the total delay such that it "matches" the chosen filter order. The problem is that there is no formula to help you choose such an ideal delay; you simply have to find out by trial and error. As you've already mentioned, filter stability is also determined by the choice of the total delay. You have to choose the total delay such that the given desired phase response can be reasonably approximated by a stable IIR filter of the chosen order. Also here there are no formulas for how to choose the total delay. With a bit of trial and error you will gain some intuition for it.

In sum, if $\phi(\omega)$ is your desired phase response, you should specify the phase to be approximated as

$$\tilde{\phi}(\omega)=\phi(\omega)-\tau\omega\tag{1}$$

where the delay parameter $\tau$ (which can be positive or negative) is an additional design variable.

The reason why an additional constant group delay term can make such a difference is that each pole and each zero contribute to the total group delay. The closer to the unit circle they are the bigger is their contribution. So for a given filter order, if the total delay is too small, some poles will move outside the unit circle causing a negative contribution to the group delay. This can be viewed as a non-causal approximation of the desired response. However, since we can only implement causal IIR filters, the corresponding filter is unstable. Increasing the total delay will move all poles inside the unit circle, hence the designed causal filter will be stable. Then the approximation error can be minimized by choosing the ideal total delay such that it best matches the chosen filter order. Note that knowing the filter order is not sufficient for choosing an optimal total delay, because the latter will depend on the exact positions of all poles and zeros.

The following example hopefully illustrates these points. The desired frequency response is low pass with a constant group delay in the pass band. The pass band edge is $\omega_p=0.2\pi$, and the stop band edge is $\omega_s=0.25\pi$. So the desired frequency response is

$$D(\omega)=\begin{cases}e^{-j\omega\tau},&0\le\omega\le\omega_p\\
0,&\omega>\omega_s\end{cases}$$

where $\tau$ is a design parameter. I chose a numerator order $M=16$ (i.e., $16$ zeros), and a denominator order $N=6$ (i.e., $6$ poles away from the origin of the complex plane). It turns out that for $\tau<17$ samples, the designed filter is unstable, i.e., there are poles outside the unit circle. For $\tau=17$ samples, we get a stable filter, but with two poles very close to the unit circle (pole radius $r=0.99576$). This filter is very hard to implement, and its frequency response has a large peak close to the pass band edge (see blue curve in the figure below). Choosing a slightly larger desired pass band group delay of $\tau=19$ samples decreases the pole radius to $r=0.98725$, and hence removes the frequency response peak (see green curve below). Finally, further increasing the desired group delay will further decrease the maximum pole radius, but the approximation error will increase. For $\tau=30$ samples the maximum pole radius is $r=0.97253$, but the approximation error is considerably larger (red curve below).

Also have a look at these answers to related questions: link 1, link 2, link 3, link 4.

You might also be interested in this article, which presents an algorithm that automatically chooses an optimal delay.


"Ensemble averaging ... cannot track dynamic changes"? A book claims this as a motivation for introducing exponential averaging:

  A disadvantage of ensemble averaging is that the resulting estimate cannot track dynamic changes occurring in the observed signal.

  -- L. Sörnmo and P. Laguna, Bioelectrical signal processing in cardiac and neurological applications

The authors seem to use the term ensemble averaging to refer to averaging all observations of a stochastic process. Modifying this method to keep a moving average, i.e. averaging only the past $n$ observations instead of all of them, seems very close at hand.

Such a moving-average smoother would be able to track changes in the stochastic process being observed, just like the exponential-average smoother. Is this the difference in dynamic capabilities that the authors are referring to, or have I missed a deeper difference?

$$ \begin{align}
Y_{\mathrm{ensemble,}N}(z) &= \frac{1}{N}\left(z^{0}+z^{-1}+z^{-2}+\cdots+z^{-(N-1)}\right)X(z) \\
Y_{\mathrm{exponential},\alpha}(z) &= \frac{\alpha}{1-(1-\alpha) z^{-1}}X(z)
\end{align}$$

You're misinterpreting ensemble averages.... this is really in "batch" like after recording a few runs. You average functions on a time interval, e.g. that's where those doctor's charts for kids' height/weight mean & spread over time arise. Time-averaging is totally different...


Does the hold part in sample and hold turns the signal from discrete to continuous? I am trying to link sample and hold to discrete and continuous time domain.
Could someone please confirm:

When sampling but not holding for "an extended" period of time, you have discrete samples.

When sampling AND holding you lose the "discretion" and obtain an approximated continuous signal.

Is this correct?

The sample-and-hold circuit is a physical, analog electrical circuit and, by definition, its output is continuous.

The actual conversion to discrete time happens at the quantizer, not at the sampler. The output of the quantizer is a discrete sequence of numbers associated with time instants, so it can be said to be a discrete signal.


How to Reduce Noise in a Video Stream? I've been using Open Broadcaster to stream the infrared HDMI output from a 4K Sony video camera. It works pretty well, but one thing that I'm finding a little annoying for image analysis is the presence of grayscale salt and pepper noise in the footage.

Does anyone know if there's a way I could reduce this? I've been trying to change the video bitrate (currently at 40,000), but are there potentially other things I could do (such as play with the downscale filter)? Thanks so much.

Usually, Salt and Pepper Noise is done using Median Filter.
There are variations of the Non Local Means to deal with Salt and Pepper but it might be overkill so I'd start with the above.

The tricky thing is to take advantage of having a Video and not just a static single image.
Hence you can also apply the filter in a temporal form.
If you use small kernel you might be able to pull it as it is.
For larger kernels you might need to use Optical Flow to estimate which parts of the frame should be compared.


How to use Kalman filter for altitude prediction based on barometer data? I have barometer noisy data with known variance.
I studied Kalman filter but I did not find an answer to this problem:
My process model is: altitude is changed because of velocity that is changed because of acceleration that is normally distributed.
 s[k+1]=s[k]+v[k]*dt+a[k]dtdt/2
Is my process state just (altitude) or (altitude and velocity) or (altitude, velocity, acceleration)?
When I use velocity and acceleration - how shall I fill the measurement matrix when measuring just altitude? In some examples, I have seen 0 as a measurement of unknown variables, but it makes no sense to me, because velocity and acceleration remain constant, because no correction is applied to them.
My goal is to mainly compute velocity.
Is there some other recommended algorithm to estimate the vertical velocity?

UPDATE:
I will try to ask the more specific question about Kalman filter :-)

I have state X = (s, v, a) - trajectory, velocity, acceleration
I have state transition model F = ((1 dt dt*dt/2)(0 1 dt)(0 0 1)) (where dt is interval between last and current measurement)
I have one-dimensional measurement Z (I measure trajectory only)
So observation model H is (1 0 0)?

I have seen examples where only acceleration was measured and H was (0 0 1), but when acceleration was corrected, velocity and trajectory was updated also in "a priori state estimate" because velocity is dependent on acceleration. But will the velocity and acceleration be updated when they are not dependent on trajectory? See bold zeros in F matrix.

You might get good results when you consider the acceleration as an input $u$, so the model could then be written as

$$
x[k+1] =
\begin{bmatrix}
1 & \Delta t \\
0 & 1
\end{bmatrix} x[k] +
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix} u[k] + w[k]
$$

$$
y[k] =
\begin{bmatrix}
1 & 0
\end{bmatrix} x[k] + v[k]
$$

where $x[k]$ a vector with the position and velocity and $w[k],v[k]$ zero mean Gaussian white noise. Here $v[k]$ has a covariance $R$ equal to the variance of the barometer and $w[k]$ probably has the following covariance $Q$

$$
Q = \sigma_u^2\,
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix}
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix}^\top
= \sigma_u^2\,
\begin{bmatrix}
\frac{\Delta t^4}{4} & \frac{\Delta t^3}{2} \\
\frac{\Delta t^3}{2} & \Delta t^2
\end{bmatrix}
$$

where $\sigma_u^2$ the variance of the measured acceleration (this assumes that there are not disturbance forces acting on the system).


Is phase and amplitude information necessarily lost when undersampling? Is phase and amplitude information necessarily lost when undersampling if you have a constant periodic single frequency sinusoidal?

My second question is: How can one determine the undersampling frequency for such a case? And how to retrieve the amplitude and phase.

I want to know if this is possible for a sampling frequency much lower than the Nyquist rate.

From an RF perspective... since you specified just the one sine wave then yes, Nyquist holds because you are operating at baseband. Any sampling rate less than 2B will result in aliasing, because you are undersampling. In this case your bandwidth is the frequency of the sine wave, since it is the highest frequency you want to recover.

However if you have a low frequency signal modulated onto a high frequency carrier signal, you can get away with undersampling because the carrier has no information. Can use spectral folding to recover the baseband signal and discard the carrier.

Now about phase recovery, there are multiple methods to accomplish it. One is to use complex signals, where you create an I/Q signal by complex mixing. With an I/Q signal you can recover instantaneous phase because you have a reference, the other copy of the signal.


Vinyl Crackling Effect I want to know how can I achieve this effect:

https://www.youtube.com/watch?v=DQg9M_47c3A

Do I need to use a filter? And what can I use for the crackling, poping effect?
Thanks in advance!

Several approaches exist to create a vinyl crackle effect [1]. Usually, this type of disturbance is assumed to be an additive one. Hence, the crackling effect is achieved by somehow generating/obtaining a disturbance signal and adding it to the clean signal, not by filtering the original signal.
If you just want to get hold of such crackle signal, you can find at least one generator online, for example here. If you don't mind to pay money for software that can simulate not only vinyl disturbances but a variety of effects, have a look at iZotope Trash 2. (This company used to offer a vinyl sound generator called "Vinyl" for free, but unfortunately this is discontinued.)
Furthermore, there are a number of free samples online (e.g. at freesound.org) that can be used.
If, on the other hand, you're interested in generating a plausible sounding but artificial vinyl click/crackle signal on your own, I would recommend having a look at an article by Välimäki et al., describing, among other things, how to do this.
[1] S. V. Vaseghi, Advanced Digital Signal Processing and Noise Reduction. Chichester, West Sussex: Wiley, 2008.


What is the difference between natural response and zero input response? I am new to DSP and was going through different responses of a system subjected to an input. My understanding of zero input response is: it is the response/output of the system when the input signal is set to zero. In other words if a system is described by a linear constant coefficient difference equation the zero input response is the homogeneous solution.

However if the $\mathcal Z$-transform of the input is a rational function $X(z)=N(z)/Q(z)$ and that of the LTI system function is $H(z)=B(z)/A(z)$ and the system is initially relaxed, then $Y(z)= H(z)X(z) = N(z)B(z)/A(z)Q(z)$. Assuming distinct zeros(real  only) and poles(real only) of $X(z)$ and $H(z)$ then

$$Y(z) = \sum_{k=1}^N \frac{A_k}{1-p_kz^{-1}} + \sum_{k=1}^L \frac{Q_k}{1-q_kz^{-1}}$$

which gives

$$y(n) = \sum_{k=1}^N A_k(p_k)^{n}u(n) + \sum_{k=1}^L Q_k(q_k)^{n}u(n)$$

where $p_k$ and $q_k$ are the poles of system $H(z)$ and input signal $X(z)$ respectively and $u(n)$ is the unit step function. Now the first term is referred to as the natural response of the system $H(z)$. It's very confusing to grasp the difference between zero input and natural response.

Edit: The reference of the question is to book DSP : Principles , Algorithms and Applications by John Proakis and D Manolakis
pdf of the book is here
Page no 203 and 204. The two paragraphs after formula 3.6.4 explains the difference between zero input response and natural response

Thank you Peter and Matt for your answers and comments.

First it's important to realize that many authors use the terms zero-input response and natural response as synonyms. This convention is used in the corresponding wikipedia article, and for instance also in this book. Even Proakis and Manolakis are not entirely clear about it. In the book you quoted you can find the following sentence on page 97:

[...] the output of the system with zero input is called the zero-input response or natural response.

This suggests that the two terms can be used interchangeably. Further down the page, we find the following sentence:

Thus the zero-input response is a characteristic of
the system itself, and it is also known as the natural or free response of the svstem.

Again, this strongly suggests that the authors believe that both terms are equivalent.
However, on the pages you mentioned they appear to make a difference between the two. And the difference is as follows. The zero-input response is the response which is caused by non-zero initial conditions. It only depends on the system properties and on the values of the initial conditions. The zero-input response is zero if the initial conditions are zero.
The natural response is the part of the total response the shape of which is only determined by the poles of the system, and which doesn't depend on the poles of the (transform of the) input signal. The natural response does depend on the input signal in terms of constants but its form is entirely determined by the system's poles. Unlike the zero-input response, the natural response does not vanish for zero initial conditions.
The total response of the system can be written as the following two sums:

zero-input response + zero-state response
natural response + forced response

The zero-state response is the response for zero initial conditions, and the forced response is the part of the response the form of which is determined by the form of the input signal.
I hope this becomes clear in the following example. Let's investigate the following system:
$$y[n]+ay[n-1]=b^nu[n],\qquad y[-1]=c\tag{1}$$
where $u[n]$ is the unit step sequence. The total response can be computed using $\mathcal{Z}$-transform techniques:
$$y[n]=\left[\frac{1}{a+b}b^{n+1}+\left(c-\frac{1}{a+b}\right)(-a)^{n+1}\right]u[n]\tag{2}$$
The zero-input response is the part of the total response that is determined by the initial condition and that does not depend on $b$:
$$y_{ZI}[n]=c(-a)^{n+1}u[n]\tag{3}$$
Obviously, $y_{ZI}[n]=0$ for $c=y[-1]=0$, i.e., for zero initial condition.
The natural response is the part of the total response the shape of which is determined by the system's pole:
$$y_N[n]=\left(c-\frac{1}{a+b}\right)(-a)^{n+1}u[n]\tag{4}$$
Note that it depends on the initial conditions as well as on the input signal (via the constant $b$).
Also note that it is the shape of the zero-state response that depends on the poles of the system as well as on the poles of the input signal transform. All other responses mentioned here only depend on one of the two sets of poles. The shapes of the zero-input response and of the natural response depend only on the system's poles, whereas the shape of the forced response is determined by the poles of the input signal. The expression for $y[n]$ quoted in your question from Proakis and Manolakis is the zero-state response (because the system is initially at rest), and the first sum is the forced response, and the second sum is the natural response. Since the zero-input response is zero in this case, the sum of natural response and forced response (i.e., the total response) equals the zero-state response
In mathematical terms, the natural response is the homogeneous solution of the difference equation, where the constants are determined such that the sum of the particular solution (the forced response) and the homogeneous solution satisfy the given initial condition. Clearly, the zero-input response is also a solution to the homogeneous equation, but the difference with the natural response is that the zero-input response alone satisfies the initial conditions, because it is combined with the zero-state response, which assumes zero initial conditions. On the other hand, the natural response alone does not satisfy the initial conditions. The initial conditions are satisfied only by combining the natural response with the particular solution of the difference equation (the latter being the forced response).
As mentioned above, we can write the total solution as
$$y[n]=y_{ZI}[n]+y_{ZS}[n]$$
(zero-input response plus zero-state response)
and as
$$y[n]=y_N[n]+y_F[n]$$
(natural response plus forced response). For the given example, we have
$$y_{ZI}[-1]=y[-1]$$
i.e., it is $y_{ZI}[n]$ that takes care of the initial condition. That's also why $y_{ZI}[n]=0$ if the initial condition is zero. $y_{ZI}[n]$ must satisfy the homogeneous equation
$$y_{ZI}[n]+ay_{ZI}[n-1]=0,\qquad y_{ZI}[-1]=y[-1]$$
So if $y[-1]=0$, $y_{ZI}[n]=0$ for all $n$. The natural response also satisfies the homogeneous equation, but not with the initial condition $y_N[-1]=y[-1]$. What is satisfied is $y_N[-1]+y_F[-1]=y[-1]$. This is why the natural response is generally non-zero, even for zero initial conditions. And the natural response is the homogeneous solution which we need to combine with the particular solution (forced response) found in the standard way. We usually have no direct means to find the specific particular solution which, when combined with the special homogeneous solution represented by the zero-input response, will give the complete solution of the difference equation. For this we need another homogeneous solution, and this is the natural response.
Again using the above example will hopefully clarify this. For an exponential forcing signal, the standard (and most straight-forward) way to obtain a particular solution is choosing a scaled version of the forcing function:
$$y_p[n]=Ab^n\tag{A1}$$
(for the sake of simplicity I leave out the unit step $u[n]$, assuming that we consider $n\ge 0$, unless we talk about the initial condition). The constant $A$ is determined by plugging $(A1)$ into the difference equation:
$$Ab^n+aAb^{n-1}=b^n$$
giving $A=\frac{b}{a+b}$. The general form of the homogeneous solution is
$$y_h[n]=B(-a)^n\tag{A2}$$
Of course $y_h[n]=0$ (i.e., $B=0$) is one specific solution, but that's not the one we're looking for. We need to determine the constant $B$ in such a way that the sum of the particular and the homogeneous solution satisfies the initial condition:
$$y[-1]=y_p[-1]+y_h[-1]=\frac{A}{b}-\frac{B}{a}$$
From this equation we get
$$B=\frac{a}{a+b}-ay[-1]$$
which shows that the homogeneous solution we need is non-zero if $y[-1]=0$. $y_p[n]$ and $y_h[n]$ found in this way are identical to the forced response and the natural response, respectively, as shown in $(4)$ and - implicitly - in $(2)$.


realtime software 90 degree phase shift a sinusoidal signal I would post this to stack overflow but it seems to me more like a signals and systems problem, which EE's thrive on.

I need a realtime system implemented in software that can generate a sinusoid that matches the input signal's frequency but with a phase lag of 90 degrees. The input signal will be a pure tone but can drift a bit in frequency (0.4 to 0.8hz) over time. Eventually I'll want output amplitude to scale with input amplitude as well but that's for another day. I'm aiming to be able to respond to input signal changes such that I'm under 5% phase error within about 10-15 seconds (a few input signal cycles).

My first attempt at this has used a hilbert transformer FIR filter to give me the "instantanous" relative phase between my input signal (100hz sample rate) and output, then a PI controller with set point =-90 to drive frequency f_out of output_signal = math.sin(2*math.pi*f_out*t). My hilbert transofmer FIR has a latency of 2.5 seconds and the dead time seems to prevent me from having any kind of a P constant without instability.

I'm currently doing a deep dive into phase locked loops (I'm new to them) but it's not clear to me yet if this is the right tree to bark up.

I'd appreciate any manner of direction/advice.

If you have it in a buffer and have enough computing power, you can track it to quite significant accuracy using two DFT bins and the formulas shown here:

A Two Bin Solution

For minimum latency I would recommend a frame length near 2 and 1/2 cycles, though for a clean signal 1 and 1/2 cycles should also give good results.

As you do a series of sliding windows, this will generate a sequence of parameters.  A moving exponential average of the parameters should give you a good current value (when they are slowly varying) to use.  More sophisticated techniques are available.

Once you have the current parameters, you can reconstruct a signal with any phase lag and/or amplitude adjustment you want.

Or if you are tracking the phase value frame to frame, the most accurate fit is at the center of the frame ($t_f$ is at $n=N/2$).

$$ x(t) = A( t ) \cos( \Phi( t ) ) $$

Where $ \Phi( t ) $ is a (predictive) near linear function and $A( t_f )$ comes directly from the two bins solution.


What factors should I consider when choosing an IIR filter realization? IIR filters may be implemented as direct form 1 or 2. But how do you decide which form to use? What are the pros and cons of each?

First of all, a bit from wikipedia on Direct Form I and II implementation.

Direct form I requires more memory, but is a somewhat simpler strategy, and is less likely to have round-off and resonance problems.

Direct form II requires less memory, but it has the potential for unusual interactions, larger numbers, and more round off error. Much of this can be reduced by cascading smaller filters, especially 2nd order filters.


Sampling and reconstruction of signal in Matlab I'm trying to write a program in Matlab that samples (using Nyquist theorem) and recovers signal. I got stucked on recovery part...recovery signal doesn't match with the original one (see photo).

And I have to make graph that shows every sinc separately (before the sum) like on photo. Also I have to use formula from photo.

This is my code:

clear; clc;
subplot(3,1,1);
f = 30e3;
fs = 2 * f;
Ts = 1/fs;
t1 = 0:1e-7:5/f;
x1 = cos(2 * pi * f * t1);
plot(t1,x1);
hold on;
t2 = 0:Ts:5/f;
x2 = cos(2 * pi * f * t2);
stem(t2,x2);
xlabel('Time');
ylabel('Amplitude');
xr = zeros(length(t1));
for i=1:length(t1)
    for j=1:length(x2)
        xr(i)= xr(i) + x2(j)*sinc(2*fs*t1(i)-j);
    end
end
plot(t1,xr);
xlabel('Time');
ylabel('Amplitude');
legend('x(t)','x(nT)','x_r(t)');

What I need to change and add?
Thanks in advance!

You are pretty close. Here is a hint: you need to make sure that your sinc pulses are lining up with your samples. You can check this by breaking it down and plotting individually the sinc pulse train that you are getting. Make sure that these line up as shown in your picture and you should be good to go.

%% Sampling and reconstruction demo
clear,clc,close all;

%% Parameters
F = 30;     % frequency of signal [Hz]
Fs = 2*F;   % sampling rate [Hz]
Ts = 1/Fs;  % sampling period [sec]

%% Generate "continuous time" signal and discrete time signal
tc = 0:1e-4:5/F;        % CT axis
xc = cos(2*pi*F*tc);    % CT signal
td = 0:Ts:5/F;          % DT axis
xd = cos(2*pi*F*td);    % DT signal
N = length(td);         % number of samples

%% Reconstruction by using the formula:
% xr(t) = sum over n=0,...,N-1: x(nT)*sin(pi*(t-nT)/T)/(pi*(t-nT)/T)
% Note that sin(pi*(t-nT)/T)/(pi*(t-nT)/T) = sinc((t-nT)/T)
% sinc(x) = sin(pi*x)/(pi*x) according to MATLAB
xr = zeros(size(tc));
sinc_train = zeros(N,length(tc));
for t = 1:length(tc)
    for n = 0:N-1
        sinc_train(n+1,:) = sin(pi*(tc-n*Ts)/Ts)./(pi*(tc-n*Ts)/Ts);
        xr(t) = xr(t) + xd(n+1)*sin(pi*(tc(t)-n*Ts)/Ts)/(pi*(tc(t)-n*Ts)/Ts);
    end
end

%% Plot the results
figure
hold on
grid on
plot(tc,xc)
stem(td,xd)
plot(tc,xr)
xlabel('Time [sec]')
ylabel('Amplitude')

%% Sinc train visualization
figure
hold on
grid on
plot(tc,xd.'.*sinc_train)
stem(td,xd)
xlabel('Time [sec]')
ylabel('Amplitude')


Calculating the extent of spectral leakage As a follow up to this question -

Is there a way that I can calculate the amount of spectral leakage in a given bin, i.e., given an input amplitude, what percentage of the amplitude is present in the adjacent bins?

Or given an input frequency and the frequency range of each bin, can I calculate the amplitude present in each bin?

Windowing in the time domain results in a convolution with the transform of that window in the frequency domain.  So you need to compute the result of this convolution to see the effect (which some call "leakage").

You can sample the transform of the window function used (rectangular, Von Hann, etc.) to determine the FFT frequency response of a windowed pure single frequency sinusoid.  For a rectangular window, the magnitude of the apparent "leakage" in an FFT result will be samples of a normalized Sinc function (more precisely a pair of normalized periodic Sincs or Dirichlet kernels), sampled at integer offsets from the frequency bin center offset, multiplied by the sinusoid's magnitude.

An attempt at the equation for the full Dirichlet pair frequency response result is on my dsp web page here: http://www.nicholson.com/rhn/dsp.html#4

But using just a normalized Sinc as the window transform equation is usually close enough, except for the bins right next to DC or 0 Hz and to bin N/2.

For windows other that rectangular, use their transform instead of a Sinc function, up-sampled and/or interpolated if needed.


How to calculate a significant value of a definite curve/shaded region in a graph? The attached graph below consists of $4.5$ seconds of pattern which I have recorded via a hardware device.

The $x$-axis is the timestamp value in seconds and the $y$-axis is the hardware value. The values are taken on $0.01$ seconds of interval.

As seen in the graph, there is a definite pattern/curve from $0.5$ to $3$ seconds. I tried calculating a definite value using the below steps:

Taking 2.5 seconds block starting from $0$ Blocks: $0 - 2.5$, $0.01 - 2.51$, $0.02 - 2.52$, etc.
Now divide this $2.5$ second block into $5$ parts of $0.5$ seconds interval blocks. $(2.5/5 = 0.5)$
So, first block of $0 - 2.5$ will give: $0, 0.5, 1.0, 1.5, 2.0, 2.5$
Now I have calculated the slope of these values $(x_2-x_1/y_2-y_1)$. This will give me $5$ slope values.
Lastly I have added these values, $(\textrm{Slope}_1 + \textrm{Slope}_2+\ldots + \textrm{Slope}_5)/5$ to get a single average value.

The thing is this value is definite for the time block of $0.5 - 3$ seconds (the curve which I want to identify), but sometimes it comes for other time blocks also.

Also, second issue is, In the current graph, most of the values start from near $0$ or at $0$ on the $y$-axis. In real time, values can start from $-1$ or $1$ also, but the curve will still be the same, I have tested it.

I think the only way is to calculate the value of the curve, as the $y$ values are not definite, but the curve shape is. In real time I would have around $30$ seconds of data and I need to check if any $2.5$ seconds from those $30$ seconds contains a curve like this.

You can use cross-correlation to detect a known signal in a received/measured signal. Here is an example using python and numpy/scipy:

pattern = lambda t: 0.8*(abs(t-1) < 1).astype(int) * np.sin(-2*np.pi*t/2)
Fs = 100.0
t = np.arange(-1, 15, 1/Fs)
t_samp = np.arange(-0.5, 2, 1/Fs)

sigma = 0.02  # add some noise to the measurement

rx = -1+0.1*(t-3)+pattern(t) + pattern(t-4)+pattern(t-9) + np.sqrt(sigma) * np.random.randn(len(t))
patsamples = pattern(t_samp)

plt.subplot(221)
plt.plot(t, rx)
plt.ylim((-4,4))
plt.title('The received signal')

plt.subplot(222)
plt.plot(t_samp, patsamples)
plt.title('The pattern you look for')

plt.subplot(212)
metric = np.correlate(rx, patsamples, 'same')
plt.plot(t, metric)
plt.title('The correlation metric')
plt.grid(True)

peaks = scipy.signal.find_peaks_cwt(-metric, np.array([2.5])/Fs, min_snr=0.5)
peaks = np.array(peaks)
peaks = peaks[metric[peaks]>40]
plt.stem(t[peaks], metric[peaks])
plt.tight_layout()

As you see, this works also when there is an offset or a constant slope in the received signal. In the correlation metric you see peaks, which you need to detect. There is the scipy.signal function find_peaks_cwt, however it is not very intuitively usable. For better/easier implementation of find_peaks you can have a look at this page.


How does salt and pepper noise occurs in an image How does salt & pepper noise occur? And, what are some techniques to remove it?

Salt-and-pepper noise is a form of noise sometimes seen on images. It presents itself as sparsely occurring white and black pixels.

In another words ( in the sense of pixels), salt and pepper noise means that are high frequencies, so for salt noise the values of this noise type is high (255 ... 200), and for the pepper noise the values of this noise type is low (5 ... 0).

To salt noise there are good filters like harmonic mean filters
to salt and pepper noise.

$$
   \hat{f}(x,y) = \frac{mn}{\sum_{(s,t) \in S_{xy}}\frac{1}{g(s,t)} }
   $$

There are the contra harmonic filter that depends of a $Q$ argument, e.g., for salt noise $Q$ should be negative and for pepper Q should be positive.

$$
   \hat{f}(x,y) = \frac{\sum_{(s,t) \in S_{xy}}g(s,t)^{Q+1}}{{\tiny\sum_{(s,t) \in S_{xy}}}g(s,t)^{Q} }
   $$

Another aproach is the order statistics filters, e.g., $\operatorname{max}$
 filter are good for pepper noise, and $\operatorname{min}$ filter are good for salt noise.

$$ \hat{f}(x,y) = \operatorname{max}_{(s,t)\in S_{xy}}{g(s,t)}$$
$$ \hat{f}(x,y) = \operatorname{min}_{(s,t)\in S_{xy}}{g(s,t)} $$

And another one is alpha trimmed filter that also are a order statistics filter type.

$$ \hat{f}(x,y) = \frac{1}{mn-d}\sum_{(s,t) \in S_{x,y}} g_r(s,t)$$

Beside this there are Adaptive Mean Filter.


Lorentzian signal fit A plot showing two S21 peaks through a resonator is given below:

As you can see, the first resonance peak occurs at 4.43GHz. The fit applied to this peak is given below.

$$y = A - 10 \log_{10} {\left[ 1.0 + \left( 4Q^2 \left(1-\frac{\omega_d}{\omega_0}\right)^2\right) \right]}$$

This fit gives a nice curve and physically sensible values of $A$, $Q$ and $\omega_0$.

The second peak occurs at a frequency close to 4.99GHz. This peak is different to the first, it is a peak sharply followed by a trough.

In the case of the second resonance peak at 4.99GHz, are there any tricks for fitting such a peak properly? My main concern is determining the true $A$, $Q$ and $\omega_0$ values for this peak.

Cheers for any ideas.

Also, if anybody has any ideas why such a resonance peak takes the form of that at 4.99GHz, please share.

What you have for your 2nd 'resonance' is actually a vibrational mode with a transfer function of the form ('s' being the Laplace variable)

H(s) = (s^2 + 2*zeta1*omega1*s + omega1^2) / (s^2 + 2*zeta2*omega2*s + omega2^2)

The denominator 'omega2' term corresponds to the peak at 4.99GHz, and the numerator 'omega1' term corresponds to the deep trough that follows the peak at 4.99 GHz.

One can demonstrate this behavior in Matlab very simply using the code below

zeta1 = .01; omega1 = 1.02; zeta2 = .011; omega2 = 1.0;

num = [1 2*zeta1*omega1 omega1^2];

den = [1 2*zeta2*omega2 omega2^2];

freqs(num,den)

The frequency response is shown in the figure below which mimics your 2nd 'resonance'.

What you need to do for your 2nd 'resonance' is to estimate 4 (not 2) parameters, viz., omega2 and zeta2 (for the peak) and omega1 and zeta1 (for the trough). Depending on your data a nonlinear least squares optimization scheme should do. As a start, you could assume zeta1 and zeta2 to be equal (say zeta) and get an initial estimate of the terms 'omega1' and 'omega2' from your graph (or data). This leaves you with solving for just 'zeta'. Then you could iterate the model to get zeta1 and zeta2, etc., etc.

Note that the 'Q' term that you used is directly related to the quantity 'zeta' in my exposition here.


Solving a Convolution Problem I was unable to solve a convolution question.Question is attached herewith.

I don't know how to initiate for solving the problem and what is the final expression that prove the convolution of  a top-hat function with itself is the triangle function.

When you slide a rectangular signal over another rectangular signal and take the area under the curve, i.e., perform convolution, at first there is no overlap. Then, the overlap gradually starts increasing in a linear fashion. When the two rectangular signals are exactly on top of each other, that is the point of maximum overlap. After than, there is a gradual decrease until there is no overlap. The resulting convolution signal is then a triangular signal.

Update: A rectangular signal is given by
\begin{equation*}
\prod (t) =
            \begin{cases}
1 &  0 \le t \le 1 \\
0 & otherwise
            \end{cases}
    \end{equation*}
Its convolution is given by
$$z(t) = \int _{\tau=-\infty}^{\infty} \prod(\tau) * \prod(t-\tau) d\tau$$
which turns out to be
\begin{equation}
z(t) = \begin{cases}
0 & -\infty \le t \le 0 \\
t & 0 \le t \le 1 \\
2 - t & 1 \le t \le 2 \\
0 & 2 \le t \le \infty
\end{cases}
= \Lambda(t)
\end{equation}


How to approximate gaussian kernel for image blur From wiki, a $3 \times3$ gaussian kernel is approximated as:
$$\frac{1}{16}\begin{bmatrix}1&2&1\\2&4&2\\1&2&1 \end{bmatrix}.$$
Applying this kernel in an image equals to applying an one-dimensional kernel in x-direction then again in y-direction, so one-dimensional kernel is
$$\frac{1}{4}\begin{bmatrix}1&2&1\end{bmatrix}.\qquad \text{for 3*3 kernel}$$
$$\frac{1}{16}\begin{bmatrix}1&4&6&4&1\end{bmatrix}.\qquad \text{for 5*5 kernel}$$
My question is how to derive this approximation?

The continuous Gaussian, whatever its dimension (1D, 2D), is a very important function in signal and image processing. As most data is discrete, and filtering can be costly, it has been and still is, subject of quantities of optimization and quantification/quantization schemes. In one 1D, the three most direct for a finite-length filter are illustrated below:

In solid blue, the continuous Gaussian. For an $L=7$-length, I choose a Gaussian with parameter  $s=\sqrt{L/2}$, in consistence with the first approximation $(P)$ given by the binomial formula or the Pascal triangle, as detailed by @Olli Niemitalo.
Each filter coefficient (red crosses) is of the form $$h_P[k] =  \frac{\binom{L}{k}}{2^L}\,.$$
Their sum is naturally equal to unity, and the coefficients are dyadic-rationals, in the form of $a/2^b$. They are practical for integer operations, on integer-valued pixels. As you can see, they are slightly above the Gaussian around the center, and below on the tails. When $L\to \infty, the shape tends to a certain bell function shape, see From Pascal's Triangle to the Bell-shaped Curve.
A second approximation (green circles), often called the truncated Gaussian, consists in sampling ($S$) $x\mapsto G(x)=\frac{1}{s\sqrt{\pi}}e^{-\frac{x^2}{s^2}}$ at integer values of $x$: $[-3,-2,\ldots,0,\ldots,3]$. They need to be normalized to 1-unit sum:
$$h_S[k] =  \frac{G(k-(L-1)/2)}{\sum_{l=1}^L G(l-(L-1)/2)}\,.$$
The third approximation ($A$), blue circles and bars) is area dependent: once again after normalization, the coefficients are proportional to the area under the Gaussian, on the interval $[k-1/2,k+1/2]$:
$$h_A[k] \propto \int_{k-1/2}^{k+1/2} G(x)dx\,.$$
The exact value depends on the quadrature formulae.
In our settings (with the chosen $s$), coefficients are, respectively from left to right (pascal, Sampling, Area):
$$
 \begin{Bmatrix}
0.0156 & 0.0232 & 0.0255\\
    0.0938 & 0.0968 & 0.0998\\
  0.2344 & 0.2282 & 0.2262\\
 0.3125 & 0.3036 & 0.2970\\
 0.2344 & 0.2282 & 0.2262\\
 0.0938 & 0.0968 & 0.0998\\
 0.0156 & 0.0232 & 0.0255\\
 \end{Bmatrix}
$$
Only the first column admits dyadic/integer operations. To work with integer arithmetic, one can round to the nearest integer. As $2^7=128$, this quantization choice (another approximation)  leads to:
$$\frac{1}{128}
 \begin{Bmatrix}
2 &  3 &  3\\
12 & 12 & 13\\
30 & 29 & 29\\
40 & 39 & 38\\
30 & 29 & 29\\
12 & 12 & 13\\
2 &  3 &  3\\
 \end{Bmatrix}
$$
With  $2^6=64$, you get:
$$\frac{1}{64}
 \begin{Bmatrix}
1 &  1 &  2\\
6 &  6 &  6\\
15 & 15 & 14\\
20 & 19 & 19\\
15 & 15 & 14\\
6 &  6 &  6\\
1 &  1 &  2\\
 \end{Bmatrix}
$$
In the case you mention in comments, the formula for $\sigma$ is a bit different: $\sigma=(3L+1)/20$, equal to $1.4$ for $L=7$. I don't know its origin yet. As said above, those are the main  principles. In the Gaussian Kernel Calculator demo, with those values, you get:
$$\begin{Bmatrix}
0.031251\\  0.106235    \\0.221252\\    0.282524    \\0.221252\\    0.106235\\  0.031251
\end{Bmatrix}
$$
which is quite close (up to the integral quadrature precision) to your:
$$\begin{Bmatrix}
0.03125\\0.109375\\0.21875\\0.28125\\0.21875\\0.109375\\0.03125
\end{Bmatrix}
$$
as the former, multiplied by 64, gives:
$$\begin{Bmatrix}
2.0001  \\   6.7990 \\  14.1601 \\  18.0815 \\  14.1601 \\   6.7990 \\   2.0001
\end{Bmatrix}
$$
The ($A$) approximation is possibly more accurate with respect to the  integral version of the convolution, with a
continuous Gaussian and a discrete or piecewise constant signal.
In 2D, one can tensorize two 1D filters, in which case the 2D filter has rank one. Or one
extends the above method to 2D.  The 2D coefficients can be obtained by two-dimensional integral estimations:
$$h_A[m,n] \propto \iint_{[m-1/2,m+1/2]\times  [n-1/2,n+1/2]} G(x,y)dxdy\,.$$
See for instance Gaussian smoothing, esp. Figure 3, Discrete approximation to Gaussian function with $\sigma=1.0$:


How do DAWs align the phase when working with plugins? I'm no signal processing guru, so I'm looking for a pretty basic explanation that I'll understand. I've seen this question about Phase Locked Loops, but I have a hard time understanding that, especially in relation to digital audio workstations (DAWs).

What I want to know is how DAWs align the phase of signals when sending them through AU or VST plugins. How do they know how much the plugins will delay a signal they send through it?

I may be mistaken in my assumption that any correction is actually necessary, but then how do DAWs correct for delay of each channel in parallel processing?

As the author of a few AUs/VSTs an AU hosting app, I'll chime in here:

There is usually no phase alignment of signals passed through effects processing (AU/VST)
The "latency" of the effect is normally reported by the effect to the "hosting" environment on initialization. This way the hosting app can keep all channels in sync (i.e "wait" until all effects have returned the same time-range of audio before mixing, etc). Core Audio SDK states that AUs should implement this property kAudioUnitProperty_Latency if their latency is > 0 seconds.
Sometimes effects do not report their latency to the host app, so there are tools like Voxengo to work around that problem.

If the hosting environment doesn't have an automatic way to handshake & compensate for latency, then its compensation is a manual affair. For example Ableton Live has a feature (17.5 Device Delay Compensation) to help work around these timing issues.


Why convolution is required, or what is the philosophy behind convolution? I am working in digital image restoration field. I have read all things about convolution, that, for an LTI system, if we know its impulse response, then we can find its output by just using convolution between input and impulse response.

Can anyone tell me that what is the main mathematical philosophy behind it? Your experience with will tell more to me than just internet surfing about it.

The Idea of Convolution

My favorite exposition of the topic is in one of Brad Osgood's lectures on the Fourier Transform. The discussion of convolution begins around 36:00, but the whole lecture has additional context that's worth watching.

The basic idea is that, when you define something like the Fourier Transform, rather than working directly with the definition all the time, it's useful to derive higher-level properties that simplify calculations. For example, one such property is that the transform of the sum of two functions is equal to the sum of the transforms, i.e.

$$
F\{f + g\} = F\{f\} + F\{g\}.
$$

That means if you have a function with an unknown transform, and it can be decomposed as a sum of functions with known transforms, you basically get the answer for free.

Now, since we have an identity for the sum of two transforms, it's a natural question to ask what the identity is for the product of two transforms, i.e.

$$
F\{f\}F\{g\} = \space ?.
$$

It turns out that when you calculate the answer, convolution is what appears. The entire derivation is given in the video, and since your question is mostly conceptual, I won't recapitulate it here.

The implication of approaching convolution in this way is that it's an intrinsic part of the way the Laplace Transform (of which the Fourier Tranform is a special case) turns linear constant-coefficient ordinary differential equations (LCCODE's) into algebraic equations. That fact that such a transform is available to make LCCODE's analytically tractable is a large part of the reason why they're studied in signal processing. For example, to quote Oppenheim and Schafer:

  Because they are relatively easy to characterize mathematically and
  because they can be designed to perform useful signal processing
  functions, the class of linear shift-invariant systems will be studied
  extensively.

So one answer to the question, is that if you're using transform methods to analyze and/or synthesize LTI systems, sooner or later, convolution will arise (either implicitly or explicitly). Note that this approach to introducing convolution is very standard in the context of differential equations. For example, see this MIT lecture by Arthur Mattuck. Most presentations either present the convolution integral without comment, then derive its properties (i.e. pull it out of a hat), or hem and haw about the strange form of the integral, talk about flipping and dragging, time-reversal, etc., etc.

The reason I like Prof. Osgood's approach is that it avoids all that tsouris, as well as providing, in my opinion, deep insight into how mathematicians probably arrived at the idea in the first place. And I quote:

  I said, "Is there a way of combining F and G in the time domain, so
  that in the frequency domain the spectrums multiply, the Fourier
  transforms multiply?" And the answer is, yes, there is, by this
  complicated integral. It's not so obvious. You wouldn't get out of bed
  in the morning and write this down, and expect that this was going to
  solve that problem. How do we get it? You said, suppose the problem is
  solved, see what has to happen, and then we'd have to recognize when
  it's time to declare victory. And it's time to declare victory.

  Now, being an obnoxious mathematician, you cover your tracks and say,
  "Well, I'm simply going to define the convolution of two functions by
  this formula."

LTI Systems

In most DSP texts, convolution is usually introduced in a different way (that avoids any reference to transform methods). By expressing an arbitrary input signal $x(n)$ as a sum of scaled and shifted unit impulses,

$$
x(n) = \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}, \tag{1}
$$

where

$$
\delta(n) = \begin{cases}
0, & n \ne 0 \\
1, & n = 0,
\end{cases} \tag{2}
$$

the defining properties of linear time-invariant systems lead directly a convolution sum involving the impulse response $h(n) = L[\space \delta(n) \space]$. If the system defined by an LTI operator $L$ is expressed as $y(n) = L[\space x(n)\space]$, then by applying the repective properties, namely linearity

$$
\overbrace{L[\space ax_1(n) + bx_2(n)\space ]}^\text{Transform of the sum of scaled inputs} = \overbrace{aL[\space x_1(n)\space ] + bL[\space x_2(n)\space ]}^\text{Sum of scaled transforms} , \tag{3}
$$

and time/shift invariance

$$
L[\space x(n) \space] = y(n) \space \xrightarrow{\text{implies}} L[\space x(n-k) \space] = y(n-k), \tag{4}
$$

the system can be rewritten as

$$
y(n) = \overbrace{L\left[ \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}\right]}^\text{Tranform of the sum of scaled inputs} = \overbrace{\sum_{k=-\infty}^\infty{x(k)L[\delta(n - k)]}}^\text{Sum of scaled transforms} = \overbrace{\sum_{k=-\infty}^\infty{x(k)h(n-k) }. }^\text{Convolution with the impulse response}
$$

That's a very standard way to present convolution, and it's a perfectly elegant and useful way to go about it. Similar derivations can be found in Oppenheim and Schafer, Proakis and Manolakis, Rabiner and Gold, and I'm sure many others. Some deeper insight [that goes further than the standard introductions] is given by Dilip in his excellent answer here.

Note, however, that this derivation is somewhat of a magic trick. Taking another look at how the signal is decomposed in $\left(1\right)$, we can see that it's already in the form of a convolution. If

$$
\overbrace{(f * g)(n)}^\text{f convolved with g} = \sum_{k=-\infty}^\infty{f(k)g(n-k)},
$$

then $\left(1\right)$ is just $x * \delta$. Because the delta function is the identity element for convolution, saying any signal can be expressed in that form is a lot like saying any number $n$ can be expressed as $n + 0$ or $n \times 1$. Now, choosing to describe signals that way is brilliant because it leads directly to the idea of an impulse response--it's just that the idea of convolution is already "baked in" to the decomposition of the signal.

From this perspective, convolution is intrinsically related to the idea of a delta function (i.e. it's a binary operation that has the delta function as its identity element). Even without considering its relation to convolution, the description of the signal depends crucially on the idea of the delta function. So the question then becomes, where did we get the idea for the delta function in the first place? As far as I can tell, it goes at least as far back as Fourier's paper on the Analytical Theory of Heat, where it appears implicitly. One source for further information is this paper on Origin and History of Convolution by Alejandro Domínguez.

Now, those are the two of the main approaches to the idea in the context of linear systems theory. One favors analytical insight, and the other favors numerical solution. I think both are useful for a full picture of the importance of convolution. However, in the discrete case, neglecting linear systems entirely, there's a sense in which convolution is a much older idea.

Polynomial Multiplication

One good presentation of the idea that discrete convolution is just polynomial multiplication is given by Gilbert Strang in this lecture starting around 5:46. From that perspective, the idea goes all the way back to the introduction of positional number systems (which represent numbers implicitly as polynomials). Because the Z-transform represents signals as polynomials in z, convolution will arise in that context as well--even if the Z-transform is defined formally as a delay operator without recourse to complex analysis and/or as a special case of the Laplace Transform.


FSK analysis: only ever find one frequency I've trying to write this software as part of my amateur radio and SDR hobby for analyzing various two-level FSK signals.  I'm very new to DSP but have been having a lot of fun so far learning more about this subject.  My long-term idea is to try and automatically detect the mark and space frequencies, as well as the baud rate, in arbitrary signals (maybe this is incredibly naive -- I don't know).

Right now I'm trying to get this software working using a known signal, known mark and space frequencies, and known baud rate.  I've tested it with samples from sources like the Emergency Alert System, Bell 103, and even a pure $2600\textrm{ Hz}$ sine wave just to make sure that my frequency detection algorithm (very basic zero-crossing detection) works correctly.

I get back a $2600\textrm{ Hz}$ frequency when I test against the pure $2600\textrm{ Hz}$ sine wave sample so I'm at least somewhat confident about my ability to detect one particular frequency.  However, I can never extract the mark and space frequencies when I use two-level FSK samples.

What I'm using right now is a sliding window that starts at the beginning of a sample and moves along the sample, changing both the starting position of the window and the overall length of the window, in an effort to try and find a good fit.  But no matter what, I never see the mark or space frequency (actually, sometimes I will get a somewhat near read on one of the two frequencies but never both).

The project is pretty small if you want to have a look for yourself and give me some pointers (no pun intended).  Apologies for the amount of code in the Cli project right now.  I've put the code that I'm pretty confident about into their own classes and behind interfaces but my "work-in-progress" code is currently just living in the Main() method of the command-line interface.

Here is my current "window algorithm" in a nutshell:

Start with window positioned at $0\textrm{ ms}$ into the sample
Start with window length of a fairly minimal value
Count zero crossings across the current window to determine the frequency of the signal in that window
Output various information about the window and the frequency it contains
Increase the window length and repeat until the window length hits a maximum (maybe $1-10\textrm{ ms}$ in $0.25-1.0\textrm{ ms}$ increments -- all easily configured)
Once the maximum window length is reached, reset to the minimum window length, increment the window position a little, and repeat the varying window length frequency samples again
Repeat this until the maximum window position is reached

And here is a screenshot of the software running, where it was able to accurately detect one of the two FSK frequencies but not the other:

EDIT 1: Turns out I had a bad sample.  When I looked at a frequency spectrum plot I only saw one frequency peak, whereas a 2-level FSK spectrum should have two peaks. I found another sample and now my window algorithm is working much better.

EDIT 2: Still using the same algorithm but now with the new sample and a few very minor tweaks to the software am successfully decoding at least the first non-trivial amount of this NWR/NOAA SAME signal.  You can very clearly see the 1010 1011 preamble repeating at the beginning.

EDIT 3: Refactored a lot of the code and cleaned-up the Cli project so that it is actually just a Command-Line Interface.

You can make a very simple frequency discriminator by multiplying your received signal with a delayed version of itself, followed by a low pass filter. Adjust the delay to set the frequency range of your discriminator. This works because a multiplier followed by a low pass filter is a phase detector (multiply two sinusoids with same frequency and different phase and you will get the sum and the difference of the frequencies: the sum is the doubled frequency to filter out with the low pass filter and the difference is the phase). A change in phase versus a change in time is frequency, therefore the low pass filtered multiplication of the signal with a delayed version of itself will be proportional to the frequency over that particular interval. Choose your low pass filter and delay based on your symbol rate.

Optionally you can make a tuned discriminator with two pass band filters with your signal fed into the input of each one. The BW of the filter depends on your symbol rate. Follow each filter with a magnitude detector, and subtract one detector output from the other. The output of the subtractor will be your demodulated FSK.


Problem designing a specific FIR filter Consider an LTI system whose impulse response is $$h[n]=\frac{1}{2^n}u[n]+\frac{1}{3^n}u[n]$$
The input signal to this system is $x[n]$ and is null for $n<0$ but may or may not be null for $n=0$.

I found the difference equation, which turned out to be $$y[n]-\frac{5}{6}y[n-1]+\frac{1}{6}y[n-2]=2x[n]-\frac{5}{6}x[n-1]$$

Then the instruction goes as follows:

  It is desired to compute the output $y[n]$ for $n=0,1,...,1000$ with the input signal being the same $x[n]$. Determine the impulse response $\hat{h}[n]$ of the minimum-sized FIR system whose output $\hat{y}[n]$ equals $y[n]$ for $n=0,1,...,1000$ and find the difference equation of that FIR filter.

I don't know for certain how to deal with this. My first attempt consisted of establishing that

$$\hat{h}[n]=\left\{
 \begin{array}{ll}
  h[n]  & \mbox{if } 0\leq n\leq 1000 \\
  0 & \mbox{otherwise }
 \end{array}
\right.$$

I don't know if that is correct. I just thought that, due to the causality of both $h[n]$ and $x[n]$, the convolution between them would be exactly the same for the first 1001 values if $h[n]$ were cut off at $n=1000$. I named $\hat{h}[n]$ that cut-off version of the first impulse response. Nevertheless, I'm really doubtful about this reasoning. On top of it, I don't even know if that response is the shortest possible as I am asked to find in the exercise.

Assuming that the expression above for $\hat{h}[n]$ is okay, it can be expressed as $$\hat{h}[n]=\frac{1}{2^n}(u[n]-u[n-1001])+\frac{1}{3^n}(u[n]-u[n-1001])$$

After some algebra, that equation becomes impossible to deal with. The $\hat{H}(z)$ has a 1001st grade polynomial in its numerator and a 2nd grade one in its denominator. Because the filter is a FIR one, those poles should be cancelled by two of the zeros of the numerator, but it's just impossible to find the zeros of a 1001st grade polynomial, so I really don't think that this is the way to solve this problem.

Any ideas?

Thanks for your time!

Think about how the output signal $y[n]$ is computed:

$$\begin{align}y[0]&=x[0]h[0]\\
y[1]&=x[1]h[0]+x[0]h[1]\\
y[2]&=x[2]h[0]+x[1]h[1]+x[0]h[2]\\\vdots\end{align}$$

From the first equation above you know that the first coefficients of the FIR filter must equal $h[0]$. The second equation gives you $h[1]$ as the second coefficient, etc. So your intuition was correct: the FIR filter's impulse response is just a truncated version of the original impulse response.

For finding the difference equation, you don't need to compute any zeros or poles, you just need the coefficients:

$$y[n]=x[n]h[0]+x[n-1]h[1]+\ldots +x[n-1000]h[1000]=\sum_{k=0}^{1000}h[k]x[n-k]$$


MOD-N Circular Convolution How to find MOD-2 Circular convolution for the two sequences $h =[-1,3,-2,1]$ and $x = [1,-1,-2,1,3,2,1,2]$.

I know the answer is $7$ $0$ from matlab but I don't know how to find it graphicly or mathematically

Write $h(z) = -1 + 3z -2z^2 + z^3$ and compute $h(z) \bmod (z^2 - 1)$, that is,
divide $h(z)$ by $z^2 - 1$ and take just the remainder.  While this looks
very complicated, if you think about it a bit you will see that all you
are doing is dividing $[-1\quad3\ -2\quad 1]$ into $[-1\quad 3]$ and
$[-2\quad 1]$ and
adding the shorter vectors to get $[-3\quad 4]$.  Repeat for
$x = [1 \ -1\ -2\quad 1 \quad 3 \quad 2\quad 1 \quad 2]$ to add
four vectors of length $2$ to get
$[3\quad 4]$.  This is presumably not too hard to do in MATLAB though I am
not familiar enough with the syntax to suggest
specific commands.  Next, compute the cyclic convolution
of $[-3\quad 4]$ and $[3\quad 4]$ preferably without invoking MATLAB functions.
The result is
$$[\{(-3)\times 3 + 4\times 4\}\quad \{(-3)\times 4 + 4 \times 3\}]
= [7\quad 0]$$

Mathematically, what you are doing is computing $h(z)x(z) \bmod (z^2-1)$
which can be done fancily by first finding $h(z)x(z)$ using FFTs and what
have you followed by the $\bmod (z^2-1)$ computation (this effectively
chops up the long vector into short pieces and adds them), or more simply
by first computing $\hat{h}(z) = h(z) \bmod (z^2-1)$ and
$\hat{x}(z) = x(z) \bmod (z^2-1)$ (the chopping up into shorter vectors
and adding them) and then computing the cyclic convolution
$\hat{z}(z)\hat{x}(z) \bmod (z^2 - 1)$ which is easy to do.

  Chop-add-convolve is easier than convolve-chop-add


Is it possible to spread signal using sum of more than one orthogonal code I was wondering if we can spread the signal using more than one spreading code. Assume we have the Walsh code used to spread the transmitted signal, the Walsh code is, for example, gotten based on Hadamard matrix:

V = [1  1  1  1
     1 -1  1 -1
     1  1 -1 -1
     1 -1 -1  1]

where every column of matrix $V$ represents one spreading code we can choose and they are all orthogonal. where, $V_1$ = [1 1 1 1], $V_2$ = [1 -1 1 -1], $V_3$ = [1 1 -1 -1] and $V_4$ = [1 -1 -1 1].  So, if I have $V_{21} = V_1 + V_2$, what's the advantages and disadvantages of using $V_{21}$ as spreading code since it's orthogonal also on all codes except $V_1$ and $V_2$

thank you

The disadvantage is that the signal V21 = V1 + V2 = [2, 0, 2, 0] disappears every other chip interval which makes maintaining carrier and phase synchronization more difficult. Also, while the total energy in V21 is equals to the sum of the energies of V1 and V2, V21 uses 4 times the instantaneous power needed by either V1 or V2 and so the transmitter needs to have greater power-handling capacity which increases the cost. Also, one problem that received a great deal of attention forty years ago was the "near-far" problem which had to do with the fact that different spread-spectrum users received with different power levels caused problems when synchronization was not perfect and the strict orthogonality used in the OP's question was not quite valid; practical stuff that modern engineers don't need to deal with in their MATLAB programs.

Advantages? I will let someone else take on the task of enumerating the advantages of the OP's scheme.


frequency response of a time domain subroutine in matlab I implemented a sigal processing structure using very basic operation like multiply and add, I wonder is it possible to find a built-in function in matlab to treat this as a black box and obtain its frequency response. I imagine the process being:

Feeding in windowed sinusoid with magnitude 1 at different frequency.
Get output sinusoid fft at the corresponding frequency location and obtain magnitude and phase.

The easiest would be to do an FFT of the impulse response. Something like this
%% calculating the transfer function of a time domain LTI model

% create an FFT grid
fftLength = 1024;
% make a unit impulse
diracDelta = zeros(fftLength,1); diracDelta (1) = 1;
% run your "black box" function, to create the impulse response
h = myFunction (diracDelta); % plus whatever arguments you need
% FFT
transferFunction = fft(h);

You can certainly encapsulate this in a function but you need to pass in your black box function as a function handle and need to figure out how to manage additional arguments and state keeping (for filters with state). Given that the code above is so simple, it doesn't seem worth the bother.
Keep in mind that this ONLY works for LTI system. For anything else, the transfer function is undefined.
The choice of the FFT length is important: to avoid time domain aliasing it needs to be long enough to contain the entire impulse response or at least where it's decayed to "low enough".


How were windows originally conceived? I am aware of the common types of windows, (Hamming, Hanning, Kaiser, Tukey, etc etc). However while many books describe them - almost none tell me just how exactly they were derived.

What is so holy about the hamming window? What about the hanning? I understand that they all play on the ratio of mainlobe width VS sidelobe attenuation, but how exactly were they derived?

The motivation for my question, is because I am trying to figure out if one can design their own windows, that also play off main lobe width and sidelobe energy.

This is only a partial answer, but there's a lecture online where Hamming talks about how he came up with his eponymous window. Starting at roughly 15:15 gives the full context.

With a reasonably entertaining story, he credits John Tukey with inventing the theory of windows (for spectrum analysis). However, he introduces the whole subject in the context of using Lanczos sigma factors to reduce Gibbs phenomenon. In addition, in The Art of Doing Science and Engineering (based on the same lectures), he describes how his window is a variation on the Hann window, which he claims was used by von Hann in economics (long before its application in signal processing). That suggests the history goes much further back, depending on how you want to define it.

The book where Tukey first named the Hamming window is The Measurement of Power Spectra from the Point of View of Communications Engineering. Given Hamming's assertion that Tukey invented the theory of windows, it would probably be a good place to start for a deeper understanding of how to design new ones. I think the book is just a reprint of Part I and Part II of his Bell System Technical Journal article, so it's available online.


How to apply watershed to segment images using matlab? How to segment this image using watershed to retrieve only the people in the image ?

I have done the following so far :

Calculated a gradient

Calculated the watershed transform

My code:

  clear;
I=imread('inpaint.jpg');
I=rgb2gray(I);

hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);

figure, imshow(gradmag,[]), title('Gradient magnitude (gradmag)')
L = watershed(gradmag);
% Lrgb = label2rgb(L);
 figure, imshow(L), title('Watershed transform of gradient magnitude (Lrgb)')

I have been successful in apply the watershed

Wanted to know can i outline my objects in the original image so that it appears segmented ?

Recall that the Watershed transform treats its input as a topographic map, and simulates flooding that topography with water. The "catchment basins" or "Watershed regions" are then the parts of the map which "hold water" without spilling into other regions.

The gradient magnitude is a poor segmentation function as-is; the noise and open contours lead to an extreme oversegmentation of the image. We can try a series of morphological operations with the intent of creating approximate foreground and background markers, and use these to remove the spurious parts of the gradient.

%# Normalize.
g = gradmag - min(gradmag(:));
g = g / max(g(:));

th = graythresh(g); %# Otsu's method.
a = imhmax(g,th/2); %# Conservatively remove local maxima.
th = graythresh(a);
b = a > th/4; %# Conservative global threshold.
c = imclose(b,ones(6)); %# Try to close contours.
d = imfill(c,'holes'); %# Not a bad segmentation by itself.
%# Use the rough segmentation to define markers.
g2 = imimposemin(g, ~ imdilate( bwperim(a), ones(3) );
L = watershed(g2);

This works OK. You get both groups of people and their shadows as regions, with a little bit of noise.

Can you elaborate on your goals? That is, will you be segmenting many different images or just images highly similar to this example? Do you need to ignore the shadows and separate the two overlapping people?

I will try to update the answer if you respond to these questions.

Segmentation Overlay

You asked how to overlay a segmentation. One way is to use the watershed lines to specify pixels in the original and set them to a bright color.

boundaries = L == 0;
I(boundaries) = 255;


Is it possible to approximate given skewness and kurtosis values using an IFFT? I have time domain data for a signal that behaves randomly but has non-random (non-Gaussian?) skewness and kurtosis values of ~0.9 and ~7.4, respectively.

The FFT of the signal shows that it's not broadband random -- it's generally flat below 100 Hz and tapers off at about -3 dB/oct above 100 Hz up to a Nyquist of 400 Hz.

I would like to create a synthesized signal that is statistically equivalent (or approximately so) to my reference signal, with nearly identical

Frequency content
RMS
Skewness, and
Kurtosis.

Additionally, I have extrapolated PSD amplitudes up to 10 kHz that I would like to add into the mix and invert to the time domain.

I have done this in the past satisfying only (1) and (2) by using, for example, MATLAB's randn() function to create random white noise (namely to create random phase angles), applying an FFT to the white noise, scaling its frequency components to match a desired FFT spectrum, and performing an IFFT to convert to time domain.  The rand() function approaches a skewness and kurtosis of 0.0 and 3.0, respectively, so my initial thought was to create white noise with a specified skewness and kurtosis.

I have tried this same procedure using MATLAB's pearsnd() function, which allows one to specify skewness and kurtosis values, but the above procedure doesn't seem to preserve the skewness and kurtosis of the scaled time domain signal.  I believe the IFFT assumes a Gaussian distribution.  Is there a work-around for this?

Thanks!

When you know that the process that generated the data is non-linear (and you are lucky to be in full control of the acquisition), you can try Attractor Reconstruction.

This technique attempts to reconstruct the trajectory a system may be taking through phase space, that results in the complex signal that is being recorded. The signal itself, may be showing chaotic behaviour with seemingly unpredictable spikes or a trace that looks like broadband noise. However, through Taken's theorem it is possible to find a "delay embedding" within which this behaviour is "explainable".

In practice, this looks like two (for example) pointers over the signal, separated by some distance (or delay) $k$. Given a one dimensional signal $x[n]$, this results in a two dimensional signal composed of $x[n], x[n+k]$ with $n,k \in \mathbb{N}, k>0$. If $k$ (and the number of dimensions, here 2) is chosen correctly, the subsets of points (here $x[n], x[n+k]$) as plotted in a Cartesian plain, would reveal an attractor.

So, given this analyical framework, the task of "synthesizing a signal" becomes that of recovering the attractor, capturing it in a model and then playing it back at will, to generate the data.

As you may suspect, the choice of $k$ and number of dimensions is the main headache here. Theoretically, this should be equal to the fractal dimension of the attractor (which is not always known) but there are many different metrics and methods for estimating that.

An example of the whole process is available in this paper or, much more extensively in Kantz's book.

Finally, in more recent years, that "modeling" part of the attractor from above could not have escaped the Machine Learning treatment, so for the more up to date approach, you could refer to something like this one or even a GAN approach. These are everywhere these days.

You can see the...attraction here, especially when you are in control of the experimental conditions. I think that something like this would capture the dynamics you are after in a better way and manage to reproduce the conditions you are (were (?)) facing.

Hope this helps.


Controlling audio high/mid/low in frequency domain Is it a realistic mapping to divide the frequency bins in 3 parts, high/mid/low frequencies?

How can these bins then be manipulated to control the level of these bands in the time domain?

Is it sufficient to adjust these bin segments and do an iFFT to time domain?

while you could do this in the frequency domain, there are many problems with this approach, some of which are outlined here.

A better approach is to use a time domain filter. There is a tutorial here.


How to choose good features for my machine vision application? I am trying to create an application (with OpenCV3) that can distinguish sidewalks from streets by a real-time camera that is embedded in a glass. I have 3 question about this project:

Which features are better to use? (For example cars and road lines for street, people for sidewalk).
Which OpenCV functions may be usable for this project (I know this is a very bad question, but I am newbie and seeking some idea and guide).
How many hours of camera capture do I need for this? Can one do a rough estimation?

Thank you in advanced!

The choice of the algorithm is dependent on the constraints of the application. If you are able to collect (or find) thousands of labeled images covering many possible scenarios, then the popular approach would be to train a deep neural network for semantic segmentation (sidewalk, road, car and etc.). The features are learnt and are not hand-crafted. There are many packages such as Tensorflow or Torch for doing this, as well as trained models. This is one of the technologies powering autonomous driving these days.

If you would need a faster solution and do not have the resources to invest in such data collection, then I would say it is still doable because roads and sidewalks can look very different (depends which sidewalks of course). As the texture of the two classes will be quite dissimilar, I would say that any classical algorithm, capable of capturing texture features could do a decent job. If we are thinking of using OpenCV, then I would immediately try with a HoG (histogram of oriented gradients) + SVM (support vector machines) pair to extract features and classify, respectively. There is already an OpenCV demo (link 1, link 2) for performing object detection in a similar way. Similarly, one could do a sliding HoG extraction + binary classification in order to do a segmentation.

Of course one needs to gather some amount of data for both approaches, but I would say much less for the second one. Rough estimation is difficult, as I haven't seen your sidewalks or roads and cannot predict the variation in the data.


Binary classifier for a small range of the audio spectrum I am trying to detect a pitch in a narrow range of the audio spectrum with minimal samples. This corresponds to the rattling of certain mechanical systems. For example, the hum of an engine.

To do this, I am hoping to identify if there is a frequency at at a certain range such as 550hz-555hz.

I spoke with a learned fellow about this problem and he mentioned that there is a variation of the FFT that only targets the desired frequency range. What is this called?
I was hoping to avoid an amplitude (volume) training period by distinguishing peaks from white-noise in the frequency domain. Can anybody point me to a tunable parameter that will enable me to identify if a certain range has a peak? I was thinking of comparing the max height to the average.
I was wondering if anybody knew of a good method to do this.

It could have been Goertzel's algorithm, though it looks at a single frequency rather than a specific band.

Another approach is to apply modulation techniques to shift the central frequency of your range of interest into the baseband, aka "zoom FFT".

Your intuition about the max to average ratio is good. Another "peakedness" metric is the ratio of geometric mean to arithmetic mean.


How to find the phase spectrum of a rectangular pulse? (Fourier Transform) How to find the phase spectrum of a rectangular pulse?

The Fourier transform of a rectangular pulse

$$ x(t) =
\begin{cases}
1,  & \text{for $|t| \le \tau /2$ } \\
0, & \text{otherwise}
\end{cases}
$$

is given by:

$$F[x(t)]=\tau [\frac{\sin \omega(\tau /2)}{\omega (\tau /2)}]$$

In general, the Fourier Transform, $X(\omega)$ is id a complex valued function of $\omega$. Therefore, $X(\omega)$ can be written as:

$$X(\omega)=X_R(\omega)+jX_I(\omega)$$

The magnitude of $X(\omega)$ is given by

$$\vert X(\omega) \vert= \sqrt {(X_R(\omega))^2+(X_I(\omega))^2}$$

The phase of $X(\omega)$ is given by

$$\angle{X(\omega)}=\tan^{-1}\frac{X_I(\omega)}{X_R(\omega}$$

  Question:

  1)How can we find the phase spectrum of a rectangular pulse as there doesn't appear to be any imaginary part in $X(\omega)$?

  2)Why is the phase spectrum changing the way it is when we are time
  shifting the rectangular pulse from the origin?

  3) Any other insight regarding the magnitude and phase spectrum is
  welcome.

For the rectangular pulse without shifting, the phase spectrum its quite simple to obtain, since its Fourier Transform $F(\omega)$ is a real function. Real functions may only lay on the real axis, thus having $0$ or $\pm \pi$ phases, depending on the positive or negative sign of the function, respectively. Try visualizing the complex plane to realize this.

Now, if you're new to DSP, for $g(t)$ it might be tricky. First, take into account the scale change when dividing by 2 in $g(t)$. This expands the spectrum by a proportion of 2, i.e. $f(t/2)$  transforms as: $2~F(2\omega)$. That's why the $sinc$ function in $g(t)$ is twice as wide as in $f(t)$. Then the amplitude is multiplied by 3, that's why $G(\omega) = 6$ for $\omega = 0$.

Due to time shifting, you'd need to multiply $F(\omega)$ by a complex exponential. Recall the time-shifting Fourier transform theorem: $f(t - t_0)$ transforms as $F(\omega)~e^{-j\omega t_0}$.

The amplitude spectrum isn't modified by a time shift (since $|e^{-j\omega t_0}| = 1$), but the phase spectrum is added to $-\omega t_0$, which is the phase of the complex exponential (i.e., $\angle e^{-j\omega t_0} = -\omega t_0$).

Finally, you'd only need to add to the phase spectrum of $f(t)$ the linear function of $\omega$: $\angle H(\omega) = -4\omega$ (since $t_0 = 4$ in this case), which has a negative slope.

That's why you have in Fig. (d) the sum of a negative-slope linear function and a rectangular pulse train of $\pm \pi$ values.

This is basic DSP, I'd recommend to take a look to these theorems at Wikipedia for a very first glance: https://en.wikipedia.org/wiki/Fourier_transform#Properties_of_the_Fourier_transform


Illumination Normalization Paper: Help understand equation I hope this is the right place to post this question.
Im reading a paper titled "Illumination Normalization Based on Weber’s Law With Application to Face Recognition"
Im only interested in Illumination normalization part of paper, the algorithm for the illumination normalization part is shown below;

WLD is Weber local descriptor.Alpha is "tune-able" parameter as described in paper. Algorithm is pretty self explainatory but here is the pdf link to paper.
Q) In the WLD part, I dont understand (highlighted part) i(del-of-x) and similarly, i(del-of-y). Isn't the del gradient ? Why would one take gradient of a pixel position ?
It is my first time reading a paper like this, it would be highly appreciated if someone would take their time out and help me understand this. I already have coded this but im not sure if I implemented it correctly,I can post code here, if needed.

It is indeed a gradient but in particular a relative gradient or normalized gradient as the paper authors put it, they also give its definition:
$$
\nabla_rF(x,y) = \dfrac{\nabla F(x,y)}{|F(x,y)|+c}
$$, where $c$ is a constant to avoid division by zero (since $F(x,y)$ represents pixel intensities it is reasonable for a zero intensity value to occur, encoding a black pixel).
Because an image is represented in a discretized 2D spatial domain the gradient can be computed through finite differences around a neighborhood for each pixel's coordinates. A pixel in a rectilinear grid has at most 8 and at least 3 (for corner pixels) immediate neighbors, thus computing its relative differential involves taking 8 discrete differences, not including the difference with the center pixel itself.
Now the authors claim that this relative gradient is illumination insensitive because the image, in the context of lambertian illumination model, can be analyzed into two components: a reflectance and an illuminance component, with the latter varying slowly in a small neighborhood around a pixel and the first depending heavily on the characteristics of the face's surfaces.
$$
F(x,y) = R(x,y) \times I(x,y)
$$, where $R$ is the reflectance component and $I$ the illuminance component.
Combining the above, this leaves us with:
$$
WLD(F'(x,y)) = \arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{F'(x,y)-F'(x-i\Delta x, y-j\Delta y)}{F'(x,y)+c}\right) \\
= \arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x-i\Delta x, y-j\Delta y)}{R(x,y)\times I(x,y)+c}\right)
$$
Slow illuminance variations mean that we can locally approximate $I(x-i\Delta x, y-j\Delta y) \approx I(x,y) $.
Thus, by factoring out $I(x,y)$, ignoring $c$ and then adding a new constant $c'$:
$$
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x-i\Delta x, y-j\Delta y)}{R(x,y)\times I(x,y)}\right)\\
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x, y)}{R(x,y)\times I(x,y)}\right)\\
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)-R(x-i\Delta x, y-j\Delta y)}{R(x,y)+c'}\right)
$$
which approximates the facial pixel intensity variations while being relatively invariant to intensity variations regarding the illuminance component of the image.
The final step of taking the arctangent of the product of the double sum and the scaling parameter $\alpha$, is used as a saturation step to avoid large values of the descriptors and to measure the direction of the variations.


Downsampling an Oversampled packet I have series of amplitudes of a signal recieved at center frequency of 433.9MHz, sampling at 250000 Hz. I can see the information received is the packet I am after, just sampled too high. I believe the transmitter transmits at 3918bps. What is the best method to downsample from 250000 Hz to 3918 bps?

I have tried using a butterworth filter, which doesn't seem to work and I don't understand why. I've tried picking every 63 samples (because 250000/3918 = ~63), but didn't seem to work either. Please bare in mind I am a bit of a newbie regarding DSP.

Based on your information a packet length is about 77.2 millisecond. With 3918 bps, each packet contains about 302 bits. Frequency is shifting about every 0.255 millisecond. you get about 256 samples per bit, so 128 point FFT seems to be correct. I assume you have a real signal, so you need to evaluate only 64 bins. If you are synchronized properly, each two consecutive FFTs will provide you with one frequency, switching at the other. If you are not synchronized with the start of the first bit, only every other bin will provide you with a nice single frequency (assuming that you do not have noise and spurious signals)


Convolution in frequency domain, what is the output bandwidth If two signals do convolution in frequency domain, would the output width/BW be related to the width/BW input signals?
Would it be addition or any other?

Following a time-domain convolution analysis which relates the output signal's domain of support to the input signals' domain of supports, we could reach a similar conclusion:

Given two frequency domain signals with domain of supports (for positive frequencies) such that $X(\omega) \neq 0$ for $\omega_1 < \omega < \omega_2$ and $H(\omega) \neq 0$ for $\omega_3 < \omega < \omega_4$, then their convolution $Y(\omega) = X(\omega) \star H(\omega)$ will have the domain of support as $Y(\omega) \neq 0$ for $$ \omega_1 + \omega_3 < \omega < \omega_2 + \omega_4 $$

Therefore the bandwidth of $Y(\omega)$ is related to the bandwidths of $X(\omega)$ and $H(\omega)$ as:

$$ BW_Y =  (\omega_2 + \omega_4) - (\omega_1 + \omega_3) $$
$$ BW_Y =  (\omega_4 - \omega_3) + (\omega_2 - \omega_1) $$
$$ BW_Y =  BW_H + BW_X $$

This assumes that the domains were compact and finite. If they are either noncompact or non-finite the result should be modified accordingly.


Understanding the over sampling factor effects before a Matched Filter? I am trying to design a matched filter at the receiver side to a half sine pulse shaping function sent by the transmitter of an OQPSK transceiver. I am not sure I understand the effect of varying the over sampling factor of the ADC right before the matched filter (the minimum is twice the symbol rate). I think it may have something to do with the quality of the timing synchronization. I'm using an early-late gate synchronization algorithm, which needs at least 3 samples per symbol. I've also read in literature that the over sampling is typically a factor of 8 for QPSK systems. Does it have other effects, maybe the BER of the system?

Perhaps the answer to a similar question may be helpful.

I've spent most of my career in SONAR, so I want to see a lot of oversampling After a matched filter.

A matched filter should be sampled at the exact time corresponding to it's delay, which is unknown but can be within a range of time that is often called a range gate. Any filter output that corresponds to a maximum deflection inside of that range gate is a possible return.  The nearer those samples are in time the more likely a match in delay will occur, and the nearer to the actual time, the greater the expected deflection. So, in this case matched filter miss match can be mitigated with a higher sample rate.

The hard part is calculating the performance because the range gate is moving in time so there are several chances to pick a maximum that is a signal and several chances to produce a false alarm and the tests are on overlapping data, so one's tests are not independent.  The detector also has memory because the probability of a threshold exceeding is different if the previous test exceeded the threshold or did't exceed a threshold.  It does not produce simple closed form expressions.

So oversampling should increase performance but showing by how much, is typically done by simulation, and results will depend on the details of the filter.  Doppler introduces more complication.

Hopefully, the answer to a question you didn't ask may help.


doing sine regression to recover the paramters from a given signal I have a discrete signal. I want to do sinusoidal regression to estimate the parameters such as amplitude, phase, frequency, etc.

It is important to note I only have the signal values and have no time when these values were recorded.

How could I possibly do the sine regression?

Thanks.

If your samples are equidistant you can use a model function like
$$\hat{x}(n)=a\sin(bn+c)+d$$ where $n$ is your sample index. You can compute the parameters $a,b,c,d$ by setting up an overdetermined nonlinear system of equations.

Let $x(n)$ be your signal. Then you'll have

$$a\sin(b\cdot 0 + c) + d = x(0)\\
a\sin(b\cdot 1 + c) + d = x(1)\\
a\sin(b\cdot 2 + c) + d = x(2)\\\vdots
$$

(Of course, the equalities are only approximate equalities.) You can use a nonlinear method such as the Newton method to solve this nonlinear overdetermined system. The only problem is that you need a sufficiently good initial solution. Try to find a simple initial guess by inspection (if your data looks at all like a sinusoid ...). For $a$ and $d$ you could simply use

$$a_0 = \max(|x(n)|),\quad d_0 = \frac{1}{N}\sum_{n=0}^{N-1}x(n)$$

as initial guess ($N$ is the number of data points).


imshow command for Gaussian function I'm trying to plot 2D gaussian function using "imshow" in MATLAB. However I could not see any white points in image. I tried to change sampling values and limits of x and y values. In addition to this, I tried z.^2 and higher orders but nothing changed.Thank you for any help.What can be the problem? My code is below:

sigma = 10;
[X,Y] = meshgrid(-40:0.5:40,-40:0.5:40);
z     = (1/(2*pi*(sigma^2)))*exp(-((X.^2)+(Y.^2))/(2*(sigma^2)));
figure
surf(X,Y,z);%shading interp
figure
imshow(z);

Why not to simply use either:

imagesc(z);

Or:

h=pcolor(z)
set(h, 'LineStyle','None')

Although if you really want to use imshow then provide a set of extra parameters to scale the plotting range. It's more computational and memory efficient.

imshow(z, [min(min(z)) max(max(z))])


Pitch class profiling Currently, I am developing a piece of software for chord recognition. There are FFT and pitch class profiling on it. I tried using the pitch class profile proposed by Takuya Fujishima, but the result is not quite accurate. This may be because I misinterpreted the pitch class profile algorithm.

I don't understand the folowing formula formula:

$p(k)= [12 \;log_2( \frac {l}{N} \cdot \frac {f_{sr}} {f_{ref}})]$  $mod$  $12$

What are $f_{ref}$,and $N$ here?
Does $N$ mean the sample size of FFT (4096,1024,etc.) or the total bins of the spectrum?

I know this is an old thread, but for reference for future enthusiasts in this field:

The equation for the full PCP algorithm as described by Fujishima is as follows:

$$ PCP(p)=\sum_{ls.t.M(l)=p}^{l}||X(l)||^2$$

$$M(l)=\mathrm{round}\left[12\log_2 \left(\dfrac{f_\text{s} \ l}{N f_\text{ref}} \right) \right]\operatorname{mod} 12 $$

The confusion that OP has is in regards to the 2nd equation:

$N$ is how many samples are used to calculate the FFT, or typically known as the FFT size.  Note that a larger $N$ (more samples, longer time) would mean higher resolution of the FFT.

$f_\text{s}$ is the sampling frequency of the audio file you are processing.  Typically this is  44 100Hz.

$f_\text{ref}$ is the reference frequency for the pitch class you are trying to match.  This is typically a set of 12 frequencies, representing the fundamental frequencies of the 12 semitones of a classical keyboard.  For your reference, they are
[16.35, 17.32, 18.35, 19.45, 20.60, 21.83, 23.12, 24.50, 25.96, 27.50, 29.14, 30.87], starting from C0 to B0.

To understand this equation better, the collective term of $f_\text{s} l/N$ is equivalent to converting the FFT bins into the actual frequency that bin represents.  $l$ here is just the bin counts for the FFT output.  So for every bin of the FFT, the associate frequency can be calculated with this term.

Because every octave is effectively a double in frequency to the previous octave (if C0 is 16.35Hz, then C1 = 2*C0 = 32.70  C2 = 2*C1).  So to reverse this, we take the $\log_2$ of the calculated frequency from the previous step.

Now to the first equation, note the summation if for all $l$, such that $M(l) = p$.  This means for all of the different pitch classes p (0-11), we must only use bins in the FFT such that $M(l) = \textrm{the pitch classes}$.  To find such conditions, we sweep the fref for all frequencies in the FFT so we find the frequencies that contributes the most to the individual pitch classes.

Note that the $M(l)$ equation should not give you a negative number.  The reason for this is the FFT creates both positive and negative frequencies.  But for all real signals (such as audio signals), the negative frequencies is just a reflection of the positive frequencies.  You only really need to deal with the positive frequencies, which resides in the first half of the output from the FFT.  Or you can use a rFFT (such as the one in python - numpy) algorithm that specifically deals with real signals, and will only return to you the positive frequencies.

If you were to implement this algorithm, you would do the following:

Pick a pitch class profile (lets say you want to look at $PCP(0)$ to find all notes that falls under the musical C)
Now, according to the 1st equation, you need to find all $l$ such that $M(l) = 0$
You go to the 2nd equation, going through all the bins of your FFT (which is l).  You also set your $f_\text{ref}$ to be the reference frequency for the PCP bin you're after (at this case, since we are looking for C, we can set the $f_\text{ref}$ to 16.35Hz).
Find all values of $l$, such that $M(l) = 0$.  Go back to the 1st equation, now with the subset list of $l$, take the magnitude of the FFT at those specific bins ($l$), to find the final sum for pitch class 0.


Fourier transform of $ne^{-an}u[n]$ I need to find the Fourier transform of the following signal:

$$ne^{-an}u[n]$$

The answers start by using the rule of the basic signal:

$$a^nu[n] \rightarrow \frac{1}{1-ae^{-j\omega}} $$

and then uses the property of differentiation in the frequency domain. However I think there's something wrong with the way the rule is used. The answers assume that:

$$e^{-an}u[n] = a^nu[n] \rightarrow \frac{1}{1-ae^{-j\omega}} $$

Is there something wrong or am I missing something?

Here's the answer for reference:

You're right, the first Fourier transform correspondence in your reference is wrong. It should be

$$\mathcal{F}\{e^{-\alpha n}u[n]\}=\frac{1}{1-e^{-\alpha}e^{-j\omega}}\tag{1}$$

You just need to substitute $a=e^{-\alpha}$ and use the formulas you know.


Reducing noise from the same frequency band as signal Sorry for the vague question (as I'm not even quite sure what I want to do is possible), just asking for some general direction to take my research.

For a brief description, my signal resembles exponentially distributed noise, bandlimited, with a triangle shaped spectrum around the baseband. I'm looking to remove additive white Gaussian noise within the same frequency band as the signal, ideally without affecting the signal bandwidth.

Below are some simulated examples of a section of my signal and noise:

and their bandwidth

So the question is whether I can, somehow, attempt to improve the SNR, given that they share the same frequency band, and my signal is not known beforehand (though it differs from the noise in distribution and spectral shape).

EDIT Averaging successive acquisitions is not an option, as what we want to measure is a local time delay of the blue trace. As such, just averaging would interfere with the measured delay.

EDIT 2 I see people commenting that the SNR is fairly high in this case, just want to clarify: This is the current use case and we are currently happy with the performance.
In the case of our system, without going too much into irrelevant details, if we increase the spatial resolution the SNR will deteriorate (up to 10 or 20dB), so maybe the representation isn't accurate of what I was trying to convey.
I will post new figures as soon as I am back in the office, or if I get an opportunity tomorrow of a noisier case.

Once you resample to the common Nyquist rate (or nearly), you can consider it a mixture model. Iirc, the optimal detector for a Laplacian signal in Gaussian noise is the soft threshold function.


Calculating the perimeter in binary image After pre-processing, I can obtain the image like this:

Then, I want to calculate the perimeter of the bottom region (annotated with grey line).

I am making an ImageJ plguin. So I have to implement on my own. (firstly, I have to identify front region (in this case, bottom) and then calculate the perimeter)

Is there any idea or suggestion?

Provided that these are representative images, you can apply a region growing operation, somewhere at the bottom of the image (or with the seed point being user defined) which will help in extracting just the contour you want to measure. It might look something like this:

Then, apply the high pass convolution matrix of:

$$\begin{bmatrix} 1 & 1 & 1 \\ 1& -8 & 1 \\ 1 & 1 & 1 \end{bmatrix}$$

Which will get rid of areas of constant "amplitude" and leave only the sharp transition of the perimeter. The image now might look like:

Now you can either just count the white pixels (there are 2658 in this one) or trace the white curve (again, a form of region growing) by starting from one of the points for which the above convolution kernel provided a high output (that is, one of the points that is exactly on the perimeter) recursively by looking for 8-connected components.

Obviously, tracing the curve will provide a much more accurate estimate because there, you can take into account the $\sqrt{2}$ distances of pixels resting on a diagonal. If you then know your pixel size, you can also convert the estimate from pixel space to actual units of length.

Hope this helps.

(All processing done in GIMP, some white pixels outside of the perimeter where detected as well but can be easily taken out manually (even in ImageJ) or, if tracing the curve would not affect the measurement.)


MATLAB: Issue with passing data through AWGN Channel I am trying to send QAM data through AWGN channel and I have 256 samples/symbol, I have the following code:

data = randi([0 3], 1, 2^12); % Generate data
dataMod = qammod(data, 4);  % Modulate data
dataMod_u = upsample(dataMod,256); % Upsample
b = rcosdesign(0.25,8,256,'sqrt');
dataMod_f = filter(b, 1, dataMod_u); % Tx Filter

channel2 = comm.AWGNChannel('EbNo',15,'BitsPerSymbol',2,'SamplesPerSymbol',256);
ch_out2 = step(channel2, dataMod_f); % Pass through channel
ch_out_f2 = filter(b,1,ch_out2);    % RX Filter
ch_out_f2 = downsample(ch_out_f2,256); % Downsample
scatterplot(ch_out_f2)

I get the following scatter plot after passing through channel:

Why I don't see the modulated output, what am I doing wrong?

UPDATE: Scaling was not adjust for upsample/downsample.

data = randi([0 3], 1, 2^12); % Generate data
dataMod = qammod(data, 4);  % Modulate data
dataMod_u = 16*upsample(dataMod,256); % Upsample
b = rcosdesign(0.25,8,256,'sqrt');
dataMod_f = filter(b, 1, dataMod_u); % Tx Filter

channel2 = comm.AWGNChannel('EbNo',15,'BitsPerSymbol',2,'SamplesPerSymbol',256);
ch_out2 = step(channel2, dataMod_f); % Pass through channel
ch_out_f2 = filter(b,1,ch_out2);    % RX Filter
ch_out_f2 = (1/16)*downsample(ch_out_f2,256); % Downsample
scatterplot(ch_out_f2)


Changing the cutoff of an IIR lowpass filter with minimal transients? I have a 1st order lowpass IIR filter of the form

$ H(z) = \frac{b_{0} +b_{1} z^{-1}}{1+a_{1} z^{-1}} $

implemented using the direct method my code looks something like this

state2=state1;
state1=signal x-state2*a1;
output=b0*state1+b1*state2;

that has a cutoff frequency X that has been calculated from the fundamental frequency of signal x. Now given that the fundamental frequency of signal x changes and it takes me time "t" to detect this new frequency, is there any way I can reduce the transient response that I see from resetting the states to zero and changing the a and b coefficiets (by logging the states, old filter output or other data over time "t" and using this to calculate the new states rather than resetting them to 0)   ? Or would someone be able to provide me with a signpost as to what topics i should be googling?

How about running two filters parallel with each other. Both filters get the same input signal but only one filter will output the filtered signal because there is a crossfader behind them which will pass only one of the filters output.
You can change the parameters of the currently offline filter, and fade the crossfader to it. Parameter changing without any transient.


Optimum Filter Signal Detection for Non AWGN Channels I have been reading this question and it confirms that the matched filter is the maximum-likelihood receiver in the presence of additive white Gaussian noise. So  in the AWGN channel it maximizes the signal-to-noise ratio at the decision instant at the end of each symbol.

However,according to Wikipedia

  AWGN channel is not a good model for most terrestrial links because of
  multipath, terrain blocking, interference, etc.

I would like to know:

If the transmission channel could not be considered as affected by AWGN. Which filter would yield the optimum SNR? I'm particularly  interested in the noise present in sound channel (sound channel through air not underwater).
Is it possible to "adapt" matched filter to non AWGN channels?

AWGN is a specific noise model that is to varying physical circumstances, of varying utility.

Not AWGN is not a specific noise model, but you can look at it as detection not dependent on Gaussian assumptions.  There is no one kind of optimal detector for every kind of non Gaussian noise.

A good book is

https://books.google.com/books/about/Signal_Detection_in_Non_Gaussian_Noise.html?id=tvpSAAAAMAAJ&hl=en

Multiplicative noise is another class as well.