Where does the "Desired Response" come from in a Wiener filter? I'm confused right off the bat in the introduction to Wiener filters, because how can you know what the "Desired response" is? Where does it come from? And if you know it already, then why build a filter to process an input signal in the first place? Can you provide an example of a typical input $u(n)$ and desired response $d(n)$ (such that that the error $e(n) = d(n) - y(n))$?

I'm reading Haykin's Adaptive Filter Theory, start of Chapter 2.

$u(n)$ = filter input
$y(n)$ = linear filter output
$d(n)$ = desired signal
$e(n)$ = error signal: $d(n) - y(n)$

thanks,
js.

In many applications a Wiener filter is used to decorrelate the desired signal from the input to the filter. The desired system output is then the error signal, which is decorrelated from the filter's input signal.

As an example, in adaptive noise cancellation we assume that we have some noise reference signal (e.g., noise picked up by a microphone). And we have a noisy signal which we want to clean up. The noise in the signal is assumed to be correlated with the noise reference. In terms of Wiener filtering, the noisy signal is the desired signal, the noise reference is the input to the filter, the filter output is hopefully a good estimate of the noise in the desired signal, and the error signal is the actual output of the noise canceler, which is a cleaner version of the originally noisy signal. The figure below shows the set-up (from S.J.Orfanidis: Optimum Signal Processing):

There are many other applications of (adaptive) Wiener filters, such as adaptive equalizers in communication systems (see this answer), and echo cancelation (see this answer).

In sum, the "desired signal" is often a signal contaminated with noise or interference, and the output of the Wiener filter is an estimate of that noise or interference, which is subtracted from the "desired signal". This is possible if at the input of the filter there is some reference of the contaminating signal. The "error signal" is then the actually desired output, which is the noisy "desired signal" minus the noise-estimate at the output of the filter.


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

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

Designing the filter

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

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

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

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

The final taps I got can be found there.

The problem

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

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

Some additional information:

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

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

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

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

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

    # Phase response
    ax2 = fig.add_subplot(212)

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

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

    plt.show()

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

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

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

Explicitly making the imaginary part of your FFT coefficients 0 will change your magnitude response.

It is very clear if you consider this:
$$
X_k = a_k + jb_k\\
|X_k| = \sqrt{a_k^2 + b_k^2}
$$
Now, you're explicitly setting new coefficients to be $X'_k = a_k$, and hence $|X'_k| = |a_k|$.
Taking the IFFT of your new frequency response will obviously not be the same as the impulse response you expect.

If instead you made $X'_k = |X_k|$ and then performed an IFFT, you would get your expected filter taps. This is because your magnitude response is enforced to be the same and your phase is enforced to be zero.


How to decide whether to use AR or MA for smoothing data? Imagine I've got some offline data that I want to smooth. I could use an auto-regressive or moving-average filter of some appropriate order for conducting the smoothing. On which criteria should I base my filter decision?

MA (Moving Average) filter is a FIR filter.
In general it is given in this form -
$$ Y(n) = \sum { w_i * X(n-i) } $$

In inherently it is always stable.

AR (Auto regressive) filter is an IIR filter.
In general, it is given that AF filter is in the form
$$ Y(n) = X(n) * w_0 + \sum { w_i * Y(n-i) } $$

If you want to design a filter who's frequency response window of a very accurate shape it is more manageable or easier to achieve through FIR filters but might be difficult with IIR filters. For example if you need to do equripple filters or Linear phase filters - FIR (MA) filters could be straight in design. However, if the exact shape is not so crucial, using AR one could always produce similarly effective filter (in terms of how steeply gain drops after required cut off), with relatively lesser number of taps.

Here is an excerpt from "Digital Signal Processing" Proakis & Manolakis; this is about the comparison of FIR vs. IIR filter design. Page: 724.

  As a general rule, FIR filters are used in applications where there is
  a need for a linear-phase filter. This requirement occurs in many
  applications, especially in telecommunications, where there is a
  requirement to separate (demultiplex) signals such as data that have
  been frequency-division multiplexed, without distorting these signals
  in the process of demultiplexing. of the several methods described for
  designing FIR filters, the frequency sampling design method and the
  optimum Chebyshev approximation method yeild the best results.

  IIR filters are generally used in applications where some phase
  distortion is tolerable. Of the class of IIR filter, elliptic filters
  are the most efficient to implement in the sense that for a given sent
  of specifications, an elliptic filter has a lower order or fewer
  co-efficients than any other IIR filter type. When compared with FIR
  filters, elliptic filters are also considerably more efficient.


Oversampling a signal I’ve built a circuit with an analogue sensor (70 Hz bandwidth) connected to a 24 bit ADC, which is oversampling the signal at 10 kHz. Nyquist sampling theorem implies I should sample by at least 140 Hz to reconstruct the signal.

Currently (if I’ve understood oversampling correctly) I’m oversampling by a factor of $N = \frac{10000-140}{140} = 70$. As the minimum sampling frequency is the Nyquist rate, and oversampling is a factor $N$ times this minimum rate. My questions are:

Does this oversampling only reduce the noise introduced from the output of the analogue sensor (given the low bandwidth of the sensor) to the end of the ADC process?
How much can I reduce the oversampling (which introduces a processor overhead) by for a negligible increase in the amount of noise to the signal?
Could I do anything else to reduce the amount of oversampling, e.g increase the ADC resolution?

From that same page you link to:

  A signal is said to be oversampled by a factor of N if it is sampled at N times the Nyquist rate.

So, you're sampling at $1\cdot10^4\,\text{Hz}$, that's $N=\frac{10^4}{140}\approx 71$; no need to subtract anything.

  Does this oversampling only reduce the noise introduced from the output of the analogue sensor (given the low bandwidth of the sensor) to the end of the ADC process?

Oversampling alone doesn't reduce noise. It's what you do afterwards with the signal. Oversampling can be necessary to e.g. to be able to filter out the noise that can be shaped to your advantage.

  How much can I reduce the oversampling (which introduces a processor overhead) by for a negligible increase in the amount of noise to the signal?

You need to define "negligible"; also, you need to think about the noise is distributed in frequency domain during your sampling process; there's no general answer.

  Could I do anything else to reduce the amount of oversampling, e.g increase the ADC resolution?

I'm confused: The only way to reduce oversampling is to reduce the sampling rate.

Maybe you mean

  to reduce the necessary amount of oversampling

Well, that fully depends on your signal, and what you want to do with it, and you forgot to tell us about either!


Why Are There Two Different Common $ 3 \times 3 $ Kernels for the Laplacian? I find both of these 3x3 Laplacian kernels to be commonly used:

 0 -1  0
-1  4 -1
 0 -1  0

and:

-1 -1 -1
-1  8 -1
-1 -1 -1

Is there a reason to choose one or the other? What's the difference between them?

(My math is not good enough to understand the Wikipedia article, so I can't figure it out myself from reading that. Googling is failing me for this one.)

The Laplacian kernel is used for edge detection, that is a high pass filter.
The 2D convolution of an image with it will "detect" or highlight big differences between pixel intensities (that is why it has negative values around a bigger positive one).

The basic difference between the kernels you wrote is that the first one will only highlight horizontal and vertical intensity differences in the image, while the second one will also detect diagonal differences.

Another important characteristic is that all its values should add to zero, so that portions of an image which had uniform values should stay uniform after the filtering.


Can you increase frequency resolution of FFT without increasing window size? I'd like to use STFT for multipitch analysis.
I realise detecting the partials existing in the signal is just the beginning.
Still I have problem with it.

Let's say I have signal sampled with 'CD' frequency 44100Hz.
With window of 1024 samples I get frequency bin resolution of 22500Hz/512=43Hz.
This is enough only to discern high piano notes like:
C5 = 523.251Hz and C#5 = 554.365.

I used to think 1024 is quite a large window.
But maybe it isn't and normally larger windows are used for detecting partials?

Can frequency resolution be increased with some other method than increasing window size, which worsens time resolution?
I thought of two methods:

Method1:

Divide signal into frequency bands with bandpassfilters (for example 0-11.25Hz and 11.25-22.5Hz).
Downsample higher bands so that original high frequencies will be now low frequencies
(so do it for second band 11.25-22.5Hz -> 0Hz-22.5Hz) - not sure this is possible.
Concat resulting bins sets with adjusted labels.

Method2:

Use series of lowpass filters with increasing limit.
Perform FFT on increasing frequency ranges.
For each frequency use best possible resolution (bins from the first FFT in which this frequency was included).
This will cause low frequencies to have better resolution but I think this is ok because for higher notes
the frequency difference is grater.

I will be grateful any remarks on this issues.

I also read here:
How do window size, sample rate influence FFT pitch estimation?
about method of improving peak picking results. I think will try to use it.

If you really insist on using FFT (rather than parametric methods, which wouldn't suffer from time/frequency trade-offs), you can fake a much better resolution by using the phase information to recover the instantaneous frequency for each FFT bin. Partials can then be detected by looking for plateaus in the function giving instantaneous frequency as a function of FFT bin index. The common implementation of this technique as described in this paper will "cost" you one extra STFT (instantaneous frequency is recovered by operations on the STFT of the signal, and STFT of the derivative of the signal).

See for example the ifgram function in this Matlab implementation of sinusoidal modeling of audio signals.

Note that this won't help resolving two partials falling into adjacent FFT bins. It will just provide a much more accurate frequency estimate than just converting into a frequency the FFT bin index of a spectral peak.


Why is fourier analysis used in doppler effect? I am a high school student that carried out the experiment on the Doppler effect to therefore analyze the graph but after doing some research, I saw one of the articles saying “ Fourier analysis is a mathematical tool that allows us to extract information about the frequency spectrum of any signal.” I don’t understand, could this be used in my analysis and what does it actually mean?
This experiment was carried out under normal conditions in other words at home and the procedure of this experiment is quite simple where I would be rotating a sound source producing constant frequency in a circular and a smartphone would be recording the data using the app. The results I got shows sharp peaks and frequency ranging from 425 when away from sound source and 441 when in-front of sound source with of course ups and downs due to errors of measurements. The idea is that I am creating a report on how math is used in a real life scenario and so I will be analyzing the graph I got by using different mathematical analysis techniques and therefore this is where I am stuck and looking for what math can be applied to my modeling change in frequency report.
(For my experiment I got sharp peaks of minima and maxima incase you wanted to know what my results were like because unfortunately I have no clue=e how to post my results here)
I was also told by one of the professors from math stack exchange that it would be a good idea to ask this question here
Reference: http://www2.york.psu.edu/~kxt7/Labs/Doncheski/doppler/lab01.pdf

The fourier transform of a signal (like sound) will show you the sound in the frequency domain. In other words, the independent variable, or x axis, will be in terms of frequency (such as hertz or radians per second), and the y axis will show the magnitude of those frequencies.
One reason that the fourier transform is useful is that it will show what frequencies are present in a signal and in what quantities. When used in the experiment you described (if I understand it correctly), you could help demonstrate the doppler effect by showing that the sound source is outputting a tone of a constant pitch, and in the recording of the sound source moving in the circular pattern, the fourier transform of the recording will show that the frequency is not constant, and appears in significant quantities in the sound.
If you record the sound as the source moves in a circle, then the fourier transform of the sound should appear to be high in a range close to the pitch that the source is outputting.
A tool I frequently use for computing the fourier transform of an audio file is matlab. I am sure there are many tools online to help you compute this as well if you do not have access to matlab.


why butterworth bandpass filter changes signal in the passband? I want to apply a butterworth bandpass filter to EEG signals and it's not working so I tried with a simple signal:

t=[0:0.01:10];
s3= 9*sin(2*pi*4*t);

Still there is a problem(probably obvious) that i cannot detect. :/
Here is how I design and apply the filter with passband [3 5]Hz:

%Filtering
Fs=100;
Wp1=3;
Wp2=5;
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 4, 40);
[z, p, k] = butter(n, Wn);
[sos,g] = zp2sos(z,p,k);
Xfilt_s3= filtfilt(sos, g, s3);

This is the initial signal and the signal after butterworth application in time and frequency domain.

I expected that the signal wouldn't change after the filtering, because its frequency is in the passband.
Also, why is this amplitude change occurring?
Lastly, since I want to apply it at EEG signals, I would like to know if filtfilt function is appropriate for EEG filtering.
Thank you! :)

You can solve your problem in a number of ways:

First of all you shall better use the a,b coefficients form instead of the z,p,k pole-zero form for designing and applying the filter. Then it can be seen that problem is about the unstability of the designed filter due to very high order with the order returned from the order estimation function.

In your program, your order returns as $n=14$ from the buttord function, with the supplied pass-stop band specifications. However when this order $n$ is used in the butter filter design function with a bandpass type (implicitly specified in the Wn variable which is also returned from buttord function) then the designed filter will have an order of $2n=28$ which is too high and as you can see from its pole-zero plot it's unstable as well:

% S0 - Filter specs
Fs = 100;   % Sampling rate in Hz
Wp1 = 3;    % Pass Band freq1
Wp2 = 5;    % Pass Band freq2

% S1 - Generate the signal (4 Hz sinusoidal)
t=[0:0.01:10];
x= 9*sin(2*pi*4*t);

% S2 - Estimate the order using buttord
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 3, 40);

% S3 - Design the filter by the "estimated" order
[b,a] = butter(n,Wn);       % using Wn returned from order estimate

figure,zplane(b,a);         % display the pole-zero plot of the filter
                            % order = 28 , it's untable.
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below:

The solution is to reduce this order as much as possible. At this point it may not meet the specifications anymore but it will be stable nevertheless. I just divided $n$ by $2$ before calling the butter design function.

[b,a] = butter(n/2,Wn);
figure,zplane(b,a)          % Now it's stable...(14th order)
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below: Note that poles are dangerously close to unit circle. Yet still we call it stable.

Lets apply the stable filter to see if it works as expected:

y1 = filtfilt(b,a,x);

The result is:

Note that you don't have to use the filt filt function unless you want zero phase filtering property.


For image quantization, is it necessary to sample image? When we want to quantize an image, is it necessary to do sampling? When we know that image is a discrete signal, is it necessary to find the best samples with Nyquist rate?

Quantization of a continuous signal is the process of representing an infinite set of input values, $x(t)$, that fall into an interval $I_k=[a_k,b_k]$, with a single level $y_k(t)$, the output, that represents all the inputs to that interval with some optimality criterion.

You can quantize an analog signal without sampling it. (however don't mix this with binarization of its amplitude). To do so, you can pass the analog signal from a set of reference thresholders, whose output will be constant as long as the input signal remains in a given interval. Binarization is not necessary if you do not want to convert this into a digital signal.


FFT showing evenly spaced signals when tuned to ~56 MHz down to 25.680 I'm running GQRX on Mac OS X 10.11.6 hooked up to an RTL-SDR dongle.  I noticed when I'm tuned anywhere between ~56 MHz on down to 25.680 MHz I see evenly spaced signals on the FFT.  But when I tune to 25.679, they disappear completely.

On the higher end, the signals are faint (here I am tuned to 56 MHz):

They become more pronounced as I work my way down (here I am tuned to 43 MHz):

They are most pronounced when I'm tuned to 25.680 MHz.  But when I tune to 25.679 MHz they disappear entirely.  Here I am switching back and forth between 25.679 and 25.680 MHz:

What's causing this?  Since I see signals below 25.3 MHz when I'm tuned to 25.680 but they disappear when 25.679 I tend to think this is a bug.  If so, is it hardware or software related?

EDIT:

The signals are more pronounced tonight.  Each one has a valley surrounded by two peaks:

The one I'm tuned to in the image above has peaks at 27.03 and 27.05 MHz.  The first peak is at about -60 decibels and the second one is at -50.

These periodic spectral lines look like they might be caused by something like a sparking power supply, DC-DC converter (for LED lighting or USB power, etc.), or high powered digital drivers running at approx. 150 kHz. You might want to walk around the house or neighborhood with a pocket AM radio to see if you can DF this RF noise (determine the direction its coming from).

RTL-SDRs USB dongles based on R820T2 + RTL2832U also have a lower frequency limit somewhere around 52 MHz (where they go deaf or alias or won't lock).  26.68 appears to be near a folding point in the sample rate aliasing, e.g. right around Fs/2, where any spectrum plus its own complex conjugate can partially cancel.


What is the "waterbed effect" in control system design? I recently stumbled across some notes on the "Waterbed effect" in some notes by A. Megretski for an MIT course on "multivariate control systems". Here's an excerpt:

  A common effect, usually associated with unstable zeroes and poles of the open loop plant, makes it theoretically impossible to make certain closed loop transfer functions “small” simultaneously at all frequencies: if amplitude of the frequency response is reduced in one part of the spectrum, it may have to get larger in the other part. This effect, sometimes called the waterbed effect, can be explained mathematically in terms of integral inequalities imposed on the closed loop transfer functions. In the basis of such results is the affine characterization of all possible closed loop responses, as well as the Cauchy integral relation for analytical functions.

I don't think I've ever heard of this before. Could someone explain the effect in more practical terms?  When am I likely to encounter this effect in practice?

If I am understanding this paper please correct me if I am wrong:

A common effect, usually associated with unstable zeroes and poles of the open
loop plant, makes it theoretically impossible to make certain closed loop transfer
functions “small” simultaneously at all frequencies:

This is talking about Pole Zero Cancellation in realizable control systems. Essentially:

$$\frac{1}{s-\alpha}$$

is unstable for a step response however:

$$\frac{s-\alpha_1}{s-\alpha_2} = 1$$ where $$\alpha_1 = \alpha_2$$

which is stable; however, due to paramter variation (resistor/capacitor tolerance)  it is impossible to cancle out an unstable pole. alpha_1 and alpha_2 may never perfectly align to cancel each other out. (maybe through digital controls)

if amplitude of the frequency
response is reduced in one part of the spectrum, it may have to get larger in the other
part. This effect, sometimes called the waterbed effect, can be explained mathematically
 in terms of integral inequalities imposed on the closed loop transfer functions.

Basically, if alpha_1 increses then this "water bed effect" is caused by alpha_2 draging down the frequency response longer becofre alpha_1 zero kicks in.

essentially the Frequency repsonse would look like this if they werent matched:

--------\
         \
          \-------------

instead of this when they are matched exactly which looks like this:

----------------------------------

(That is, a flat response)

If the opposit happens (alpha_2 is made larger, you should see the opposit effect of this response)

             -----------------
             /
            /
      -----/

.

In the basis of such results is the affine characterization of all possible
closed loop responses, as well as the Cauchy integral relation for analytical
functions.

Is answered by this paper:


Windowed Fourier Transform I need some help with applying the properties of Fourier Transform.

We define the windowed Fourier Transform of $f \in L^2(R)$ as
$$Sf(\mu,\xi)=\int_\mathbb{R}f(t)g(t-\mu)e^{-i\xi t}dt$$

Prove that for $f=e^{i\eta_0t}$, we have $Sf(\mu,\eta)=e^{-i\mu(\eta-\eta_0)}\hat{g}(\eta-\eta_0)$

All I can obtain till now is this.
$$Sf(\mu,\eta)=\int_\mathbb{R}e^{i\eta_0t}g(t-\mu)e^{-i\eta t}dt=\int_\mathbb{R}e^{-it(\eta-\eta_0)}g(t-\mu)dt$$

I know that I need to apply some properties of fourier transform but I am not sure how. Thanks in advance.

HINT:

Just substitute $\tau=t-\mu$ in the integral and see what you get. It's just very basic math, no complicated properties of the Fourier transform. I assume that you know the definition of $\hat{g}(\eta)$.


Why does the unit step function not have infinite poles? Let's say I have a transfer function which is a unit step function.

$H(s) = \int_{0}^{∞}e^{-st}dt$

But when we write,

$H(s) = 1/s$ it is only true when $Re(s) > 0$

So after we derive the Laplace Transform for the unit step, we forget all about how we derived it and how $H(s) = \infty$ for all $Re(s) \leq 0$

What makes the pole $s=0$ so special. It does not lie in the ROC like every other point in the left half plane, so why is it the only pole for our system?

Shouldn't the system have infinite poles then on the left half of the s-plane? I guess the more generalized question would be isn't every point outside the ROC of a transfer function its pole?

EDIT:

If you can accept that the above proposition is TRUE then consider the following contradiction.

Consider a simple RC filter whose transfer function is

$H(s) = \frac{1}{s+1}$

This system is said to have 1 pole at $s=-1$

Applying the inverse Laplace transform, we now the corresponding time domain function

$h(t) =  \exp(-t)$

But we know that Laplace Transforms for exponential such as $\exp(at)$ is only defined when $s>a$

In our case the Laplace transform could only be defined for $s>-1$. Therefore, to say that there is only one pole at $S=-1$ would be incorrect as everything to the left of $s=-1$, including the line at $Re(s)=1$ would also be undefined and out of the ROC.

So the question again becomes why do we, when solving such systems only consider them to have ONLY one pole at $s=-a$?

Is that not an incomplete answer?

No. Every point that's not in ROC doesn't have to be a pole. ROC is a region where the transform attains a finite value, but this ROC also has to "produce back" the time domain signal via the inverse transform, or else it does not converge to that sequence in time domain. So you have to take both the finiteness of the transform as well as the reproducibility of the time domain signal via the inverse transform into account to define a "region that converges and is consistence for both aspects" while defining an ROC for a signal. I am sure you know that the same transfer function with different ROCs produces different time domain signals. So specifying a transfer function does not uniquely define a time domain signal. Only when an ROC is defined along with the transfer function can we associate to a unique time domain sequence. So the left over region that's not part of the ROC doesn't have a pole, since clearly if we shift the ROC to be that side now, we can define a different sequence.


Poster detection in OpenCV? I'm really new in image processing, so please sorry I'm a newbie. I tried to use the squares.cpp for detecting posters (since they usually are rectangles) without using expensive feature detectors (like SIFT). Unfortunately, the results are pretty much disappointing (as it was pretty predictable, results below).

Actually I don't care that only posters are detected, since statistically the posters are the biggest (or second biggest) rectangle in the image (decent heuristic).

The last image is the result of this code using this Hough Transofrm code (which seems working even worse!).

Any idea how to improve this code?

[

It seemed to work relatively well where there was a strong lightness difference at the edge of the poster, compared to background.

This matches what squares.cpp seems to do: scale the image, do a canny edge detection, simplify the resulting contours and look for squares.

There's two things that, looking at your example images, go wrong here:

The contour detection / simplification doesn't seem to deal well with reflections crossing the edge (see Cpt America poster in last photo, and the batman vs superman poster in the 4th image)
The canny edge detector will fail (and so does my brain) on posters when the frame/background difference isn't clear. I had to look twice myself at the second picture to make sure where the chrome frame ends and where the poster begins

So, this points to it being helpful to pre-process your images in different way, try to detect the rectangle, change your preprocessing, repeat, and so on.

A few observations:

maybe you can just reduce the aggressiveness of the contour simplification
Try different edge detectors! Canny is not always the best choice.
the weaker reflections just seem to be relatively white light – maybe convert your image to a brightness/hue/saturation colorspace, and reduce the variance of the brightness channel

As a side note: this approach (pyramid scaling, edge detection, contour extraction, contour simplification, rectangle detection) sounds more complex than just a normal Feature extraction by means of a transform of the image. For example, there are several transforms in OpenCV, some of which are very usable to detect dominant straight lines in the picture. If you can find parallelograms in those, you should be one step ahead, especially if you can really assume the poster is the biggest rectangle.

I think the Hough transform is quite exactly what you need.


resilient codes to phase shift I need to define an encoding to beacon IDs. The codes are going to be periodic but in the receiver side I can have phase shifts so I' wondering the best way to decode:
- 8 bits code groups that represent the same ID:

Group1: 00000001, 00000010, 00000100 .....  10000000
Group2: 00000011, 00000110, ......11000000, 10000001
Group3: 00000101, 00001010, .............., 01000001
...

How can I verify that a code belongs to a group? maybe the Real Part of a 8 bit DFT?
Is there a group/codes generator?

How can I verify that a code belongs to a group?

Assuming these are really just sub-words that are common to all elements in a group, simply trying out all 8 possible (cyclic) bitshifts of the received sequence to see whether you can match any group representative (e.g. the first column in your GroupN table) would very likely be the most efficient approach – because 8 bit is really not that much, and bit shifts are very efficient on CPUs or on FPGAs (whatever you're using), and even more importantly (since "beacons" surely doesn't sound like high data throughput), the code is easy to understand, and premature optimization is the root of all evil.


What is the purpose of wrapping the negative times of a response function in discrete convolution? I am trying to rationalize a figure given in the Numerical Recipes in C in the section of Fourier based convolution and deconvolution.
The authors show the example of a continuous convolution with a zero centered response function r(t). Convoluting the continuous signal blurs out the edges and there is no shift.
What is somewhat confusing is that they suggestion for the discrete convolution analogue of the continuous convolution, we should wrap the negative time section of the response function to the extreme right end of the array. What is the mathematical justification of doing N/2 circular shift, where N is even number of samples. The authors do not mention as N/2 circular shift, but it seems this is what they are doing.
These are their actual figures.

It seems the authors are suggesting circular shifting the zero centered response by N/2 (N is even) to get an equivalent convolution to continuous convolution. Why can't we have discrete version of the response function as shown by the squared markers which I drew in MATLAB? I do not get equivalent convolution results if we convolute a square wave with the sampled response function on the right and the left as suggested in Numerical Recipes.

It's about zero-phase filtering.
First, the continuous-time response function $r(t)$ is an even-symmetric, noncausal filter so that its Fourier transform will be a zero-phase function; i.e.,
$R(\omega)$ is a real-valued function of frequency $\omega$.
Note that there's almost a mistake in the figure 13.1.1; the output should begin before $t=0$, even for a zero-phase filter. Because this is a non-causal filter, it will anticipate future inputs before they happen, and put out the response before $t=0$.
I think, however, they wanted to justify this condition, by shifting the input $s(t)$ sufficiently to the right (or assuming that $s(t)=0$ until $t=t_1$) so that the output response before $t=0$ will be zero anyway, and thus eventually we are interested at the output of the filter beginning from $t=0$.
Now, sample the signals $s(t)$ and $r(t)$ (assuming they are bandlimited) and obtain a sampled version of the continuous-time output by performing a discrete-time linear convolution, or equivalently use a DFT based efficient implementation of linear convolution by circular convolution of proper length.
If you choose the first approach, and perform direct discrete-time linear convolution using the sequences $s[n]$ and $r[n]$, then you won't need to shift anything. But remember that a practical implementation will still need some shifts to represent the noncausal impulse response $r[n]$ in the proper index frame.
But if you choose the second approach; using DFT to implement a linear filtering, you have two choices for the filter response $r[n]$ to be registered; (because in the DFT domain, there's no negative time-index $n$, and you are restricted to the range $0 \leq n < N$) first you can linearly shift $r[n]$ sufficiently right to make it a causal filter $r_c[n]$, and follow the procedure as usual; with this choice, the effective discrete-time filter loses its zero-phase (and zero group delay) property and your effective output begins at the group delay position index.
If, instead, you wrap the negative-index portion of $r[n]$ to its far right edge, in the DFT frame of $N$ samples, where $N$ is chosen properly to aviod aliasing in the time-domain to ensure that computed circular convolution will match the desired linear one, then you will maintain the zero-phase property of $r[n]$ while performing the DFT based convolution; your first output at $n=0$ will be what you would get at $t=0$ in the continuous-time case.


Is it correct to subtract a low-pass filtered signal from the original signal and use the result as a "high-pass"? I am having a hard time trying to find documentation to implement band-pass or high-pass filter with python/scipy/numpy.

I can easily create and apply a low-pass filter, though, so I ask:

Would it be conceptually correct to low-pass-filter a signal, then subtract the result from the original signal, in order to get just the high-frequencies?

Also, if anyone has a simple example of a naive bandpass filter in Python (preferrably using the numpy and scipy libraries), I'd be very thankful.

What I look for is something like:

filtered_signal = band_pass(original_signal, rate, low=20, high=500)

Thanks for any help!

EDIT: with scipy, I'm using this as low-pass, with good results:

import numpy, scipy.signal

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    firstpass = scipy.signal.lfilter(b, a, interval)
    ## second pass to compensate phase delay
    secondpass = scipy.signal.lfilter(b, a, firstpass[::-1])[::-1]
    return secondpass

In theory you can do this, but in practice it is difficult to do because the time and phase alignment must be pretty good for it to work.  If the alignment is good you will get the destructive interference that you are seeking.  If they aren't, you will get some constructive interference.  Even worse, whether they are destructively or constructively interfering will depend on the frequency- i.e. you can get both constructive and destructive interference at the same time.  It can work, though, if you are only filtering out fairly low frequencies, since their timing requirements are the loosest because they change so slowly.

Short story- it's possible to do but is difficult enough that it generally makes sense to just do a high-pass filter.

A relatively straightforward way to create a bandpass filter is to create a low-pass filter, and then modulate it to the center frequency you want by multiplying it with a sinusoid of that frequency.


Reasoning Behind Additional Phase Shift in DPSK In MATLAB's Comm System Toolbox there is an option for 'Phase Rotation', defined as the phase difference between previous and current modulated symbols when the input is zero. This is the first time I've heard of this concept.

Let's take DQPSK as an example. Say my previous symbol had an absolute phase of $0$, (1+j0) and I am trying to encode '01'. Under the DQPSK scheme I'm familiar with, '01' corresponds to a phase shift of $\pi/2$. Therefore the absolute phase of this symbol would be $\pi/2$, (0+j1). MATLAB would describe this scheme as having a phase rotation of 0.

If the same example were to use a phase rotation as $\pi/2$ then the absolute phase of the symbol would be:

$0$ (phase of last symbol) + $\pi/2$ (shift due to modulation) + $\pi/2$ (additional phase rotation) = $\pi$

This brings me to my question. Performing the additional phase shift requires increased complexity at the receiver and transmitter. So what is to be gained from the process of adding an additional phase shift? Why do it in the first place?

A QPSK constellation could in theory have any rotation. For example, $$\mathcal{S}=\lbrace 1, j, -1, -j \rbrace $$ is a QPSK constellation, and $e^{j\phi}\cdot\mathcal{S}$ is also a QPSK constellation. Matlab's phaserot parameter in the dpskmod command is $\phi$ in the equation above.

What is to be gained depends on the particular application I guess. In particular, note that $\phi=\pi/4$ is nice because then the ML decoding rule is equivalent to finding what quadrant the received point is in.


how does the Frequency spectrum of a signal change in the logarithmic space? I am trying to implement basic homomorphic filtering but I can't seem to understand what happens in the frequency domain when you take the natural logarithm of the function.

In Matlab I'm getting not a number or infinity everywhere.

Plus is the simplest operation : $z = x+y$. Fourier is inherently linear, and good at addressing it. However, most processes and data combination are nonlinear, and they should be dealt with. The second simplest operation is multiplication. Homomorphic filtering deals with $z = x\times y$. It is more complicated, especially because of zeroes.

Logarithms were invented to linearize products. Because $\log z =  \log x + \log y$, and then we get back to linear. The problem is that logarithms are not defined everywhere in a simple way, especially with negative numbers. So, it is necessary to deal with that, and I do not know of a naturally sound method. But classically, one offsets and scales data, such as with a  modified logarithm:

$$ \mathrm{l}_m (x) = \mathrm{sign}(x) \frac{\log(1+a|x|)}{\log(1+a)} $$

related to companding. When $a$ and $|x|$ are small, this is close to $$\mathrm{sign}(x) \frac{\log(1+a|x|)}{\log(1+a)}\approx \mathrm{sign}(x) \frac{a|x|}{a}\approx x$$ thus back to the linear case.

You can also check: Why do we substract a background image and not divide it?, talking about LIP (Logarithm Image Processing).


Down sampling an EEG signal I have a set of 10 minute EEG signals that were sampled at 400 Hz and have 16 channels which corresponds to a 16x240000 matrix. These EEG signals belong to two different classes. I am trying to classify these 10 minute segments using a Neural Network, in particular a LSTM.

Since the size of the matrix is very large, even after I split them into time segments, I would like to first reduce the number of samples. My current approach for pre-processing is to use wavelet transform for denosing  and down sample to at least 200 Hz. This results in a 16x1200 matrix.

I was also looking at ICA as a feature extractor and reduce the size.

My question is, how do I tell which method is better suited for my task? Will downsampling the signal lead to a significant loss in information?

I will appreciate any suggestions. Thank you.

Edit:
 - The data corresponds to two stages of seizure:

Preictal: Right before the seizure
Interictal: in between seizure

-The data was captured from 3 different patients and using 16 electrodes, sampled at 400Hz.

Actually I am not quite sure if you should/need to downsample it with conventional methods. Did you consider using CNN (with stride larger than 1) first and then feed segment of the output to LSTM? The former step gives you some downsampling effect but it is probably better than any conventional preprocessing methods since those layers will be trainable and the whole model can be trained end-to-end (assuming that you have sufficient data).


Regarding MFCC feature of a speech signal I am trying to implement speech recognition but I'm struck with a few questions:

There are 12 Coefficients of MFCC. What are their names?
What is the range of values of MFCC coefficients?
If I want to relate MFCC coefficients graphically, what will be unit of weights on their edges?

You can have as many MFCC coefficients as you want, 12 is just only a widely used number. As you probably know (or if not then please refer to the old answer), coefficients are being obtained via fitting of cosinusoids to log energies in your filter banks. Usually 12 is enough, more is not improving recognition rate too much. It is better to calculate their derivatives ($\Delta $'s and $\Delta\Delta$'s across time). Regarding names - they are simply MFCC's (first, second, and so on). Although coefficient $c_0$ is usually called the energy.
Question about range is not straightforward to answer. Usually in my applications they are in $[-12, \ 12] $. Although if you perform the Cepstral Mean Normalisation on yours, then you can expect them to have values of Gaussian Distribution with $95\%$ of them being in range of $[-2, \ 2] $.
I do not really understand that question - you must be more specific about what you are asking. Edges of what?


Support Vector Machine: A non-probabilistic binary linear classifier I read that SVM is a supervised learning method, it is also a non-probabilistic binary linear classifier. I understand why it is binary because it classifies our training pattern to two classes $w_1$ and $w_2$ which are labeled as $y_i=+ 1$ and $-1$. And I understand why it is linear because it classifies when the classes are linearly separable. But I don't understand the non-probabilistic part, could someone clarify?

It can be said that, it is not based on a generative model which models the generation process of the data and itself a probabilistic model. SVM is based on a minimization problem, which seeks for a hyperplane that optimally separates data points of two clusters given the labels. There is nothing probabilistic in this formulation. You just formulate an optimization problem and solve it via Lagrange multipliers.

I heard (but not studied) that, there are probabilistic SVMs, maybe if you study them, you can see the difference.

Note that, at the fundamentals, authors talk about probabilistic notions, such as risk minimization, however at the macro scale, there is no probabilistic notion.


Continuity and its relationship with asymptotic spectral decay The asymptotic decay of the magnitude of the Fourier transform of a function appears always to be determined by its continuity properties as follows, with examples given in Fig. 1:

Continuous integral, -20 dB / decade;
Continuous function, -40 dB / decade;
Continuous 1st derivative; -60 dB / decade;
...

How can this be expressed mathematically in a rigorous way, and proved?

Figure 1. dB magnitude of the Fourier transform of a function as function of the base-10 logarithm of frequency for (from slowest to fastest decay) a rectangular pulse (turquoise), a single lobe of a cosine (purple), and Hann function (blue).

I suppose it would be sufficient to show that a discontinuous function in the time domain decays as $1/\omega$ in the frequency domain (e.g., $\textrm{rect}\Leftrightarrow\textrm{sinc}$), and then use the fact that each integration in the time domain adds a multiplicative factor of $1/j\omega$ in the frequency domain.


Harmonics to Noise Ratio Estimation I'm willing to estimate the Harmonics to Noise Ratio (HNR) of a speech signal x[k] and using autocorrelation method.

Theoretically, HNR is given as,

$$ \ HNR =  \frac{R_{xx}[T_0] }{R_{xx}[0]-R_{xx}[T_0]} $$

where $\ R_{xx}$ is the autocorrelation function (ACF) and estimated as,

$$ R_{xx,biased}[l] = \frac{1}{N} \sum_{k=l}^{N-1}{x[k]x[k-l]} $$

and,

$$ R_{xx,unbiased}[l] = \frac{1}{N-l} \sum_{k=l}^{N-1}{x[k]x[k-l]} $$

My question is how can I choose ACF function? Should it be biased or unbiased?

Constraints are signal x should be stationary (so I use the small part of the signal to be stationary) and noise is white and should be uncorreleated (N should be large enough). But these constraints has nothing to do the with type of the ACF estimation.

For large number of samples both will be indistinguishable.

The Biased Version is the Maximum Likelihood Estmator (MLE) of the problem.
It means it has many nice properties for $ N \to \infty $.
The Unbiased version is Minimum Variance Unbiased Estimator (MVUE) of the problem.
Yet in practice when $ N $ is large, as in your case, it won't matter much and there won't be any difference between the two.

Yet, you should remember those are for the case the mean isn't known.
In case it is known, one should always use the Unbiased Form.


Detect quasi-parallel lines with spurs in a binary image I'm looking for a way to remove quasi-parallel curves from the net, as depicted in the figure. Not only the original nets (upper row in figure) contain many spurs, but also the curves to be detected are not necessary straight lines. Ideally, I'd like to keep only the border of each such a region (lower row in figure).

I've came across solutions like hit-or-miss with line structuring element, Hough transform and Radon transform - but they were aimed at either detecting long straight lines or arbitrary shapes (like ellipses).

My idea was to represent a net as a graph, linearly approximate each edge (and split it into smaller segments if it deviates from the approximated line too much), for each pixel store the slope of an approximated edge in a matrix, and finally for each pixel I'd search such a matrix in a given neighborhood for other pixels (belonging to a different edge) with a similar slope.

Could you recommend me any method already used to accomplish such a task or tell me whether my approach seems reasonable?

I would start with a pruning algorithm (see skeleton pruning) in order to erase small segments, then use:

the Hough transform, but segments are not really straights.
if the pruning was a success, then you separate each segment and use the first and last point to compute the orientation, which will define when line are parallel. You will have to introduce a tolerance because you want quasi-parallel lines.


Enlarging a picture by adding pixels in matlab I was wondering is there a way(some code or functions) to convert an image with specific number of pixels into a picture with more pixels in matlab(enlarging the size or the number of the picture)?

Yes, you can use imread to get the image values, interp2 for the interpolation, and imwrite to make a new file.  Pay attention to data format and scaling.  Read the docs for these 3 functions and you should be ok


Reduce "flicker" in spectrum analysis / spectrogram I'm performing a series of FFTs on windowed extracts of an input signal, to get a frequency spectrum over time.  This is for a visualisation - I'd like to keep the window small, so the visualisation is fast-moving and responsive.

I'm only interested in the amplitude.  Occasionally, two or more frequencies are close enough that the amplitude flickers wildly with the interference pattern between the two.

What ways are there to reduce this flickering, that will disturb my time-domain responsiveness less than increasing the window length?

(For reference, this is the visualisation - I'm not concerned with the "flashes" when something happens in the audio, but the juddering is not good.)

My suggestion is to low pass filter the FFT slices.

Suppose your current FFTs are:
$$
X(n,k) = {\tt FFT}\{x[n] \cdot w_n[n]\}
$$
where your signal is $x[n]$ and the window $w_n[n]$ enforces how much of the signal you want to include in the FFT at time $n$.

Then form:
$$
Y(n,k) = \alpha \cdot Y(n-1,k)  + (1-\alpha) \cdot X(n,k)
$$
where $\alpha$ is a fraction between 0 and 1 (probably around 0.7 to 0.99).

This will effectively apply a first order lowpass filter to your FFTs.

Another possibility:

$$
Y(n,k) = {\tt median}\{X(n,k)\}
$$

where the median is taken over the last $M$ values of $X(n,k)$ i.e. $\{ X(n,k), X(n-1,k), \ldots, X(N-M+1, k)\}$.


Fourier transform of a periodic/aperiodic signal Generally speaking, I know that periodic signals (continuous-time domain signals) with period 2pi/wo have a spectrum with equidistance Delta-impulses of distance w0.

My question is that, if we have a spectrum with equidistance Delta-impulses of distance w0, does it necessarily means that the time domain signal is periodic? Or for example, if our time domain signal itself is not periodic, like a sin(wt) when w is not a rational ratio, will the spectrum be periodic?!

Also, does a one-sided spectrum with equidistance Delta-impulses of distance w0 correspond to a periodic signal in time domain?

Equidistant Dirac impulses in the spectrum imply a periodic time domain signal.

As pointed out in a comment, in continuous time, the signal $x(t)=\sin(\omega_0t)$ is always periodic, regardless of the value of $\omega_0$.

Your question about the spectrum being periodic is unclear to me.

A one-sided spectrum with equidistant Dirac impulses also implies a periodic time-domain signal, but that signal is necessarily complex-valued because its spectrum has no Hermitian symmetry, as required for real-valued signals.


Linear correlation of two vectors? We had a lecture on digital watermark detection using linear correlation. Here it was explained that the linear correlation of two vectors is equal to their scalar product divided by the dimension. So, the larger the scalar product, the larger the linear correlation.
This is not true, is it?
Simple example:
u = [1,2,3,4]
v1 = [10,10,10,10]
v2 = [0, 1,2 3]

The scalar product u*v1 would be way bigger than u*v2, but obviously, there is a perfect correlation between u and v2 ?
Can anyone guess how the statement was to be understood?
Thanks!

You can use corrcoef to compute the Correlation coefficients as a measure of linear dependence of the random vectors (i.e you transpose your row vectors, the variables are column vectors, and the rows are observations). With the vectors
u = [1, 2, 3, 4];
v1 = [10, 10, 10, 10];
v2 = [0, 1, 2, 3];

You have
>> R = corrcoef([u', v1', v2'])

R =

     1   NaN     1
   NaN   NaN   NaN
     1   NaN     1

There is a perfect match between u and v2 as v2' = u' - 1 showing that the two variables u and v2 are directly correlated. In the results above R is
$$
\mathbf R = \begin{pmatrix}
\frac{\operatorname{cov}(\mathbf u, \mathbf u)}{\sigma_{\mathbf u}^2}& \frac{\operatorname{cov}(\mathbf u, \mathbf v_1)}{\sigma_{\mathbf u}\sigma_{\mathbf v_1}} & \frac{\operatorname{cov}(\mathbf u, \mathbf v_2)}{\sigma_{\mathbf u}\sigma_{\mathbf v_2}}\\
\frac{\operatorname{cov}(\mathbf v_1, \mathbf u)}{\sigma_{\mathbf v_1}\sigma_{\mathbf u}}& \frac{\operatorname{cov}(\mathbf v_1, \mathbf v_1)}{\sigma_{\mathbf v_1}^2}  & \frac{\operatorname{cov}(\mathbf v_1, \mathbf v_2)}{\sigma_{\mathbf v_1}\sigma_{\mathbf v_2}}\\
\frac{\operatorname{cov}(\mathbf v_2, \mathbf u)}{\sigma_{\mathbf v_2}\sigma_{\mathbf u}}& \frac{\operatorname{cov}(\mathbf v_2, \mathbf v_1)}{\sigma_{\mathbf v_2}\sigma_{\mathbf v_1}} & \frac{\operatorname{cov}(\mathbf v_2, \mathbf v_2)}{\sigma_{\mathbf v_2}^2}
\end{pmatrix}
$$
The NaN results between (v1', u1'), (u', v1'), (v1', v1'), (v1', v2') and (v2', v1') is because v1' is constant, see MATLAB documentation on corrcoef, more specifically because
$$\begin{cases}
\operatorname{cov}(\mathbf v_1, \mathbf u) = \operatorname{cov}(\mathbf u, \mathbf v_1) = 0\\
\operatorname{cov}(\mathbf v_1, \mathbf v_2) = \operatorname{cov}(\mathbf v_2, \mathbf v_1) = 0\\
\operatorname{cov}(\mathbf v_1, \mathbf v_1) = 0
\end{cases}\quad\text{and}\quad\sigma_{\mathbf v_1} = 0
$$
In the end, each element of the matrix $\mathbf R$ is the covariance of two vectors divided by the product of their standard deviations.
EDIT:
Now let's say v1 isn't a constant vector and you have the three vectors as
u = [1, 2, 3, 4];
v1 = [11, 12, 11, 12];
v2 = [0, 1, 2, 3];

You then get
R =

    1.0000    0.4472    1.0000
    0.4472    1.0000    0.4472
    1.0000    0.4472    1.0000

Note that in the case here $\sigma_{\mathbf v_1} \neq 0$,  $\operatorname{cov}(\mathbf v_1, \mathbf u) \neq 0$, and $\operatorname{cov}(\mathbf v_1, \mathbf v_2) \neq 0$


Filter design with zero - pole placement method I have some question for you about Filter design.

Α. Calculate the transfer function in order to stop the frequency of $300\textrm{ Hz}$ for sampling frequency at $12\textrm{ kHz}$. Use the
zero – pole placement method
Β. Calculate the absolute value for the designed filter at $f_1=300\textrm{ Hz}$ and at $f_2=200\textrm{ Hz}$.

Can someone show me how can i solved this problem.

Here let me show you a simple procedure very similar to pole zero placement which will be helpful for your notch filter design.

First, lets analyse the frequency response of a single zero and let
$$ H(z) = 1 - b z^{-1} $$
be a first order system with a single zero at $z = b$ where $b$ is a complex constant with a radius $r$ and phase angle $\phi$ radians; i.e.,
$$ b = r e^{j\phi} $$

Lets see this zero on the z-plane and the corresponding frequency response magnitude for the values of $r = 0.9$ and $\phi = \pi/4$ radians:

Note that since there is a single zero, this is a non-symmetric frequency response (due to a complex impulse response $h[n]$). To get rid of this non-symmetric frequency response, we shall force the impulse response to be real and the way to achieve this is to add a second zero at the complex-conjugate location resulting in the filter:
$$ H_{az}(z) = (1 - bz^{-1})(1 - b^{*} z^{-1}) $$
as the second figure shows:

Now we have a pair of complex-conjugate zeros, whose frequency response is conjugate symmetric (and therefore the magnitude response is even symmetric as displayed)

next, we display the same figures, for a pair of complex-conjugate poles with a radius of $r=0.85$ and the same angle $\phi = \pi/4$ radians:
$$ H_{ap}(z) = \frac{1}{(1 - az^{-1})(1 - a^{*} z^{-1})} $$

So far we have displayed the pole-zero locations and corresponding frequency responses for individual pairs of poles and zeros at the same frequency. Lets combine them into a single filter and display the result:
$$ H(z) = \frac{(1 - bz^{-1})(1 - b^{*} z^{-1})}{(1 - az^{-1})(1 - a^{*} z^{-1})} $$

Nice! We have reached a system whose frequency response now resembles that of a notch, but a weak one? The solution comes by recognizing that we should better put the zero on the unit circle (i.e., set its radius $r=1$) for an infinetely deep nulling at that frequency;

Now this system is what can be called as a 2nd order notch filter with a pole radius of $r=0.85$ and radian frequency of $\omega = \pi/4$ radians per sample. This is a quite satisfactory notch filter.

Finally lets display a much sharper notch by moving the pole closer to the zero (closer to the unit circle) by setting its radius to $r = 0.99$. But note that it's very dangerous for a pole to wander about unit circle, as it can easily fall over it, thereby making the system unstable...

For your convenience, choose the pole radius between 0.9 and 0.99 depending on your numerical accuracy and how sharp a notch you need. And adjust the frequency of the pole-zero pair according to which frequency you want to null.


Trouble with inverse Z-transform and calculating of samples I have a little problem. I have to solve this task but I can't.

  Z-transform of sequence $\{x(k)\}$ describe by the formula:

  $$X(z) = \frac{2.5 -3.15z^{-1} + 1.2 z^{-2}}{1-2.3z^{-1} + 1.2z^{-2}}$$

  Calculate of samples $x(-2), x(-1), x(0), x(1)$ and $x(2)$
  double-sides sequence, which correspond this Z-transform.

Thank you very much.

The poles in the equation are Z=1.5,0.8;hence to be two sided the roc is $0.8<=z<=1.5$
so the 0.8 term is right sided and 1.5 term is left sided let us do partial fractions now $ X(Z)=1+\frac{1.5-0.85z^{-1}} {(1-0.8z^{-1})(1-1.5z^{-1})} $
after doing some maths i got $x(n)=\delta(n)-0.5(0.8)^nu(n)-2(1.5)^nu(-n-1)$
so x(-2)=-0.88888,x(-1)=-1.3333,x(0)=0.5,x(1)=-0.4,x(2)=-0.32.


Finding phase from group delay I know that group delay is defined as :

$$\begin{align} \tau(\omega) = grd [H(e^{j\omega})] &=-\frac{d[\arg[H(e^{j\omega})]]}{d\omega}\end{align}$$

But if what is available is group delay itself, is it possible to recover the phase?

According to Discrete-Time Signal Processing 3rd edition by Alan Oppenheim equation 5.86 the expression would look like this:

$$\begin{align} \arg[H(e^{j\omega})] = -\int\limits_{0}^{\omega} \end{align} grd [ H(e^{j\phi})]\, d \phi + \arg[H(e^{j0})] $$

and my doubt is:

Where did the second term $$\arg[H(e^{j0}) ]$$ came from?

Note that if $F(x)$ is an antiderivative of $f(x)$, i.e., if $F'(x)=f(x)$, then from the fundamental theorem of calculus we have

$$\int_a^bf(x)dx=F(b)-F(a)\tag{1}$$

Since the phase $\phi(\omega)$ is an antiderivative of the negative group delay $-\tau(\omega)$ we get from $(1)$

$$-\int_a^b\tau(\Omega)d\Omega=\phi(b)-\phi(a)\tag{2}$$

With $b=\omega$ and $a=0$ we get from $(2)$

$$-\int_0^\omega\tau(\Omega)d\Omega=\phi(\omega)-\phi(0)\tag{3}$$

By adding $\phi(0)$ on both sides of $(3)$ we obtain the given formula:

$$\phi(\omega)=-\int_0^\omega\tau(\Omega)d\Omega+\phi(0)\tag{4}$$


Sensor rise time compared to experimental step input rise time I have a dynamic calibration system. With it, I am applying a experimental step input to a pressure sensor. The goal is to measure the step response of the sensor and use it to calculate the frequency response.

The problem I have is that the sensor has a fast rise time. Due to the nature of what I doing, I cannot make the rise time of the input signal much faster than it already is. So I am looking for some reference on how the non-ideal input rise time of a signal affects the measured step or frequency response. Would it be reasonable to quantify error based on the difference in the frequency domain between ideal and non-ideal step inputs?

I have looked my school's library a number of times for papers on this, but I could just be using the wrong words.

So, when you remember how the step response, impulse response and frequency response are related, you'll notice that if you, instead of an actual impulse integral (i.e. a step) use something that is wider, then the frequency response simply gets windowed to lower frequencies.

In other words, and to put it as an information-gathering problem:

if the bandwidth of your excitation signal is limited (e.g. by a slow rise time of your step), then the information you can get about a system reacting to that is limited to the frequencies in the excitation.

It's often physically hard to make a large jump excitation. But maybe smaller jumps work, and you'd be better of

Generating a pseudorandom (band-limited) white signal (and save it!)
exciting your system with that
correlating the output of your system with the known white signal

The idea is this: you get the frequency response as the transform of the impulse response.

A white signal $s$ has the property ($*$ representing convolution here, ${}^*$ meaning complex conjugation, which is identity for real-valued signals)

$$s(t) * s^*(-t) = \delta\text,$$

i.e. you can "generate" an impulse $\delta$ by correlating a white signal with itself.

Knowing that your system has a step response and impulse response $h$, and knowing that convolution is commutative, we can do the following to find that:

$$\delta * h = s(t)*s^*(-t) *h = s(t)*h*s^*(-t)$$

In other words: instead of exciting your system with an impulse, which can be hard to do, to get the impulse response, we excite it with a white signal, which we can (withing the bounds of our time and bandwidth restrictions) produce digitally.

Then, we just correlate the output of the system with the same white signal, and get a result that is the same as had we put on an impulse of the same bandwidth as the white signal – but the white signal can be made arbitrarily long, so that very small amplitudes of excitation suffice to deliver the same energy as a really high, short pulse.

This is something that is commonly done to identify systems (especially: communication channels) in digital communications, but it applies to any LTI system.


Extended Kalman filter gain I'm using an EKF and some gains are bigger than 1, is that possible? In the KF I understand that the typical values are [0-1], but in EKF??

I really appreciate your response. I was studying the equations for a simple system of position of a particle at constant velocity that I show in the figure, discretized with euler. In the Figure I show the equations to make the estimation for the three possible cases for this system, where we only measure x1, only x2 or both.

In the equations of the KF I see that for example I have expressions like xk1 + = (1-W11) * (xk1 -) + W11 * (z1), ie if W11 changes from zero (the filter uses only the prediction in this case to perform the estimation) to one (the system uses only the measurement in this case to make the estimation) and in the middle we have a mixture. That's why I thought that concept was unique. I do not understand when W should be greater than one. ("Mido" means meassure).

Ok, so the question is not "Can the Kalman gain be bigger than 1?" but "Can the Kalman gain for this system be bigger than 1?"

In that case, I suggest you look at the evolution of the Kalman gain given your system and see if there are any conditions where it does. Perhaps for this particular system it can't be. You'll need to specify all the parameters of your signal model to do that.

Of course it's possible!

Why do you think the KF gain is less than 1? What the Kalman gain is depends on the system model and the data being processed.

There are plenty of examples showing gains greater than 1. See, for example, the image below taken from here.


How to decode Morse code (CW) I have two related questions on this topic.  First, about how wide should the filter be around the carrier frequency?  What is used to determine this?

Second, in Python, I must use 'numtaps' when implementing a filter.  What exactly are numtaps?  How can I determine what to set this value to?  Setting it to about the carrier frequency seems to make decoding more successful.

From another stackexchange question (https://ham.stackexchange.com/a/2018) you should use a filter approximately 4 times as wide in Hz as the expected words per minute of the morse code.  That answer also includes some links and more info.

The value "numtaps" says how long to make the filter.  You are using an FIR filter.  The sharpness and the stop band attenuation (how much the stuff you don't want to hear gets reduced) are both dependent (amongst other things) on the value of numtaps.  Higher values gets you better filter performance, but it is more work for the computer.
If you find that you need high values, then you are not using an optimal filter.  There is another factor that influences the quality of the filter, and that is the "window" that is used to shape the filter.  If you don't use a window, then the stop band will not be reduced far enough.

The length of the filter (numtaps) is also relevant to the lowest frequency you can work with.  If you try to make a band pass filter, for example, then numtaps needs to be at least as long as the lowest frequency (in samples.)  Say you are sampling at 1000Hz and want to make a band pass from 1Hz to 3Hz.  You will need for numtaps to be a least 1000.  If you used instead only 333 for numtaps then you would end up with a lowpass that lets through everything below 3 Hz.

You seem to be using scipy - here is the description of the FIR-filter designer "firwin." http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.firwin.html

Also refer to the windows page: http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window

I think you should try a blackmanharris or some other window that doesn't require parameters - otherwise you will have to start guessing numbers that you don't understand.

Simply using a window at all (other than boxcar) should improve things considerably so that you don't need such long filters - numtaps equivalent to the beep frequency should not be needed - you should be able to get away with well under 100 for numtaps


QPSK Data Rate vs Bandwidth at Passband: Why can it only transmit 1 bit per Hz and not 2 bits per Hz? 1) This is a follow up question to the most up-voted answer to this question.

2) The accepted answer is that "QPSK transmits 1 bit per Hz at passband".

3) The theory found in books seems to agree with this answer.

4) But this defies my logic and my experience.

5) My logic and experience tells me that QPSK transmits 2 bits per Hz at passband.

6) From experience, I know that in order to acquire 1 MHz of complex bandwidth (at baseband), I need to set the IQ rate of my receiver to 1 Mega samples per second. (i.e. IQ rate = Complex Bandwidth, and keep in mind that the IQ rate usually needs to be slightly higher to account for the roll-off of an anti-aliasing filter. But let's assume that there is no filter for simplicity.)

7) At baseband, this is the equivalent of 500 KHz of bandwidth if you don't consider the negative frequencies.

8) And at passband, the total bandwidth will be 1 MHz because the negative frequencies "appear" as the frequencies on the "lower" side of the carrier frequency and they are now taken into consideration on the amount of bandwidth, effectively "doubling" the bandwidth.

9) With an IQ rate of 1 Mega samples per second, I can send 2 bits on each sample (because it is QPSK).

10) So the effective data rate is 2 Mbps for a bandwidth of 1 MHz at passband (in my experience).

11) Refer to image below, describing relationship between data rate and bandwidth in my experience.

12) Other people seem to have obtained similar results to mine:

A) Refer to the last video on these student experiments: They show a transmission of 2 Mbps with a 1 MHz QPSK signal.

B) In this video a 4 KHz data rate signal results in a 2 KHz QPSK 3dB spectrum.

Let us assume you have a bandwidth $W = 1\,\text{Hz}$ in passband. This means that you have a bandwidth $B = W/2 = 0.5\,\text{Hz}$ in baseband.

In baseband, using sinc pulses, we can transmit at a symbol rate $R_p = 2B = 1\,\text{Bd}$, where $\text{Bd}$ (baud) is the SI unit for "symbols (or pulses) per second".

Using BPSK modulation, we can transmit one bit per symbol, for a bit rate $R_b = 1\,\text{b/s}$.

Now, QPSK is esentially two BPSK signals in parallel. In other words, we have two BPSK baseband signals $s_I(t)$ and $s_Q(t)$, each having a bandwidth of 0.5 herz and transmitting one bit per second.

By using quadrature, we transmit both $s_I(t)$ and $s_Q(t)$ at the same time in passband, using the signal $$ s(t) = s_I(t)\cos(2\pi f_c t) - s_Q(t)\sin(2\pi f_c t). $$

The QPSK signal $s(t)$ has a bandwdith of 1 herz and transmits 2 bits per second (one over the in-phase component and one over the quadrature component).

However, note that the relationship $R_p = 2B$ is valid only for sinc pulses, which are almost never used in practice. If you use other pulses, then the answer will be different. Many introductory textbooks are not explicit about this.

Let us say that you wish to use a pulse $p(t)$. Furthermore, $p(t)$ is a Nyquist pulse. Then, the spectrum of the pulse determines the pulse rate and subsequently the bit rate.

As an example, let us repeat the exercise above for pulses shaped like a half-sine wave. In this case, the pulse rate is $R_p = B/2.5$. Then, if $B=0.5\,\text{Hz}$, $R_p = 0.2\,\text{Bd}$ and, using QPSK, you will only be able to transmit at a bit rate $R_b = 0.4\,\text{b/s}$ per herz in passband.

Note that the $2.5$ figure I used above is my (conservative) personal measurement and reflects what I'm comfortable with. Many introductory textbooks are extremely optimistic about the actual bandwidth of different pulses. In particular, it is common for textbooks to claim that $R_p = B$ for rectangular pulses, which I believe borders on the dishonest (Stallings and Tomasi do this, among others). This is what is behind the claim that "QPSK transmits 1 b/s/Hz" that you (very sensibly) suspect is false.

Note also that the answer you linked to is correct, for the specific pulse shape that is assumed. My answer is more general and allows you to calculate the bit rate for any Nyquist pulse you choose to use.


Python: Resample without equal I would like to resample (downsample) a signal using python in order to get an even spacing and fill gaps.
About the signal:

It consists of a vector for y (amplitude) and x (timestamps)
Very slow; $F_s$ is probably >100 times higher than needed (for the desired content; steps and noise are present)
The samples aren't exactly equally spaced
The signal has gaps with missing samples
The signal is not periodic (FFT resampling should be fine though, I can trim off beginning and end)
It doesn't really matter how the gaps are filled, a linear interpolation would do

In matlab could likely just use y = resample(x,tx,fs) , however, scipy.signal.resample() can take a vector for x but still doesn't work for signals with non-uniform spacing.
It is probably a bad idea to use scipy.interpolate.interp1d() for filling gaps and downsampling in one step because of aliasing. Filtering before interpolating is likely not a good idea either, since the discrete filters wouldn't work with a non-uniform spacing of the samples.
Should I first interpolate() to the approximate sample frequency and then downsample using resample()?

Should I first interpolate() to the approximate sample frequency and then downsample using resample()?

That sounds like the best approach. Use the initial interpolation to "straighten out" the sampling grid, i.e. you should interpolate so that the grid is uniform but has roughly the same amount of samples than the original. I would go with a cubic spline for that.
If you interpolate to a grid that is an integer multiple of your target sample rate, you don't need to call resample() You can just run a lowpass filter and throw away the extra samples.  That lowpass filtering will also remove a lot of the noise and the interpolation artifacts. resample() does the same thing, but by doing the lowpass manually, you have much better control over the process.


Spectral flatness or Wiener Entropy for AR(1) and AR(2) I'm sudiying compressibility of random processes by using Spectral flatness aka Wiener Entropy

I would like to know if there is any reference which derives this quantity, for autoregressive processes AR(1) and AR(2) in terms of their coeficients.

I found the answer:

The AR(1) Gaussian process is defined as:

$$ x(n) = a x(n-1) + \sqrt{1-a^2}\epsilon (n) $$

Where $\epsilon (n)$ is zero-mean stationary white Gaussian noise.

It has PSD:
$$\gamma_{xx}(f)=\frac{\sigma^2_{x}(1-a^2)}{1-2a\cos(2\pi f)+a^2} \ \ \  |f|<1/2$$

Where $$\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}=\frac{1}{1-a^2}$$

The AR(2) is defined as:

$$ x(n) = a_2 x(n-2)+a_1 x(n-1) + \sqrt{G}\epsilon (n) $$

The definition of $\epsilon (n)$ is the same as above.

Its PSD reads:
$$\gamma_{xx}(f)=\frac{\sigma^2_{x}G}{|1-a_1 e^{-j2\pi f}-a_2 e^{-j4\pi f}|^2} \ \ \  |f|<1/2$$

Here:
$$\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}=\frac{1}{G}$$

According to this this Paper for an AR(p) Gaussian, the Spectral Flatness Measure (SFM) is very easy to compute:

$$ \hbox{SFM=}e^{-2\rho}\  \hbox{ where } \rho=\frac{1}{2}\log_2\left(\frac{\sigma^2_{x}}{\sigma^2_{\epsilon}}\right)$$

Their SFM is respectively:

$$ e^{\log_2(1-a^2)} $$
and
$$e^{\log_2(G)} $$


Ideal low-pass filter in the Frequency Domain I try to create a low-pass filter to apply it for an image.

$$H(U, V)=1 \text{, if }D(U,V) ≤ D_0 \text{, and } 0 \text{ otherwise}$$

and $D(U, V) = \sqrt{U^2 + V^2}$

U,V = np.shape(padded)
H = np.zeros((U,V))
D = np.zeros((U,V))
#cut off point
D0 = 5
for u in range(U):
    for v in range(V):
        D[u,v] = np.sqrt(u**2+v**2)
        if D[u,v] <= D0:
            H[u,v] = 1
        else:
            H[u,v] = 0
plt.imshow(H,cmap='gray')
plt.show()

After I Run it I expect this image

But I get this circle at the left top of the image as bellow :

so I think I should get the center of $\frac{U}{2}$ and $\frac{V}{2}$ but I do not know how.

How can I get the expected filter image?!

The rule should be something like this i suppose $$D(|U - U_o/2|,|V - V_o/2|)≤ D0$$ where $U_o$ is the length of image in horizontal dimension and $V_o$ is the length of image in vertical dimension


k-NN match for object recognition The $\frac{d_1}{d_2}<0.6$ gives a large number of false matches and it assumes the nearest neighbour is a correct match.

Anybody knows any good method to reduce number of false matches or filter them ?

Can we use k-NN matches instead of $\frac{d_1}{d_2}<0.6$ ratio for matching the SIFT features and then choose the best of k points by:

Changing the Euclidean distance by considering the spatial neighbourhood, or
Spatial neighbourhood matching.

Anybody ever tried something like that? Or any paper related?

Generally, one wouldn't gain much from such parameter tuning or workarounds.
The accuracy heavily depends on the image features that you are using, together with KNN algorithm. If you are using an approximate ANN, try using FLANN. In FLANN, increasing the number of trees and the iterations contribute to the accuracy.

To improve the matching accuracy, enhancements over descriptors are also possible. Zisserman proposes RootSift, which is very easy to implement and enhances SIFT with an intuitive modification:

http://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/arandjelovic12.pdf

There are also other good details in this paper.

Other than that, if you are already using BoW, you end up with sparse features. For a more meaningful representation, try VLAD (http://www.vlfeat.org/api/vlad.html). Yet, a quantization through a vocabulary tree as in Nister and Stewinius (http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.61.9520) would help.

Imposing spatial neighborhoods isn't straightforward, but can be done. However, I recommend you to first try the aforementioned methods with post processing (geometric verification). A nice one is found here: http://web.stanford.edu/~sstsai/selected_papers/2010_GeometricReranking.pdf

Cheers,


How do I calculate the bandwidth from a waveform? Suppose I have a waveform $x(t)$
where $x(t) = 5 + 30\mathrm{cos}(2000\pi t) + 10\mathrm{cos}(6000\pi t)$ that is to be uniformly sampled for digital transmission, how do I calculate the bandwidth of this signal?

I reckon that it would $4000 \textrm{Hz}$ because we're subtracting the highest and lowest frequencies $(6000-2000)$. Is this right?

However, the answer in my book says it's $3000 \textrm{Hz}$. I do not understand how did they get this answer.

Thank you for reading

A sinusoidal signal is represented as
$$x(t) = \mathrm{cos}(\omega t) = \mathrm{cos}(2\pi f t)$$

$\omega$ is the angular frequency and $f$ is the frequency. See Frequency definition.

Your signal \begin{align}
x(t) &= 5 + 30\mathrm{cos}(2000\pi t) + 10\mathrm{cos}(6000\pi t)\\
 &= 5\mathrm{cos}(2\pi \times 0 t) + 30\mathrm{cos}(2\pi \times 1000 t) + 10\mathrm{cos}(2\pi \times 3000 t)\end{align}

There are three frequencies $0, 1000, \textrm{ and } 3000\textrm{Hz}$ thus the bandwidth is $3000\textrm{Hz}$.


Conjugate symmetry Problem I'm seeing something very curious in the Matlab code below.

X = fft(x);
Xt = X;
Xt( (end/2+1):end ) = 0;
xr = 2 * real( ifft( Xt ) );
Xr = fft(xr);

where x is a real zero mean sinus function sampled at Fs Hertz according to Nyquist.

The spectra X and Xr should be the same. However, Xr(Fs/2) is not equal to X(Fs/2). Also, xr and x are slightly deviating.

Am I missing somehting here?

You are not handling DC and Nyquist correctly. These are both real and don't have a conjugate "partner" so you need to multiply them with .5 instead

Try

x = randn(1024,1);
X = fft(x);
Xt = X;
Xt( (end/2+2):end ) = 0;
Xt( (end/2+1)) = 0.5*Xt( (end/2+1));
Xt( 1) = 0.5*Xt( 1);

xr = 2 * real( ifft( Xt ) );
Xr = fft(xr);
delta = abs(Xr-X).^2;
fprintf('Diff = %6.2f dB\n', 10*log10(mean(delta)/mean(X.*conj(X))));


Deriving a Kalman Filter Equation for a Linear Gaussian Filtering Model with Non Zero Mean Noise I am trying to answer an exercise question from the book Simo Sarkka - Bayesian Filtering and Smoothing.
The question is:

Does anyone know if there is a resource that has the solutions for this book?

Since the model is linear taking care of known mean is pretty simple.
Measurement
Since the "Bias" of the measurement is known all you need is to define $ \hat{y}_{k} = {y}_{k} - {m}_{r} $. Then use it as you'd use a measurement with no bias.
Model
The model is used in the prediction step. What you should do is use the regular prediction step and add to it $ {m}_{q} $ after the prediction.
All other steps only include the covariance matrices and are the same.


What does imaginary Component of a bin's phase in the Fourier Transform mean? I have been recently working with the FFT and realized that each frequency within a bin corresponds to a phase. This concept I have no problem with because if you wanted to construct the summation of sinusoidal waves back into the original signal, you would need the phases of those sinusoidal waves.

My Questions:

Why must there be an imaginary component to each of the bins corresponding phases? It makes no sense to me why a phase would need an imaginary component.
In reference to DSP with the FFT, how does this imagery component come into play when characterizing a signal?

It simple trigonometry.  A vector (FFT result bin) can have a magnitude and a phase.  Or you can specify the same vector identically with X,Y coordinates.  The complex representation simply uses X for the real part and Y for the imaginary component.  Without an imaginary component, there’s no place to put the Y, and your vector would be under specified.

The phase ends up being atan2(I,Q) == atan2(Y,X) == atan2(im(z),re(z))

As for why use a complex representation rather than 2 strictly real quantities?  The complex math is more elegant, and requires less writing on the chalkboard to demonstrate similar things.


Baseband vs Bandpass Signal I cannot find any reliable or reasonable differences between the two even after searching for it! Please tell me the possible differences with reasons and reference links!

Starting with a Bandpass signal - this is a signal that has a bandwidth $BW$ about a center frequency $f_c$, so the signal content goes from $f_c-1/2BW$ to $f_c+1/2BW$. For a real signal (i.e. not complex) the their is also the corresponding negative frequency component i.e. the signal also exists at $-f_c$.

Basebanding often refers to shifting this signal down to DC. This is done by mixing the signal (modulating) it with $\exp(-j2\pi f_c t)$ and low pass filtering to remove the component that was originally at $-f_c$ and gets moved to $-2f_c$ due to the shifting. The component that was at $f_c$ and gets moved to DC is kept by the low pass filter. This signal is now complex due to the loss of conjugate symmetry. This resulting signal is referred to as the Baseband signal.

Note that the mixing signal $\exp(-j2\pi f_c t) = \cos(j2\pi f_c t)-j\sin(j2\pi f_c t)$, which is often how you see it implemented in communication systems.

You can achieve the same result by using the Hilbert transform which essentially removes the negative frequency component and then shifts the signal down to baseband.

References:
Haykin - Communication Systems (but almost any communications textbook), also Crochiere - Multirate Digital Signal Processing (somewhere in the first couple of chapters they explain this).


the phase spectrum of a zero-phase wavelet is not zero! I am going to plot the phase spectrum of a zero-phase ricker wavelet. I calculate the fft of the wavelet and use the phase or angle to calculate the phase spectrum but the result is not zero!

I use the following code:

Fs = 250;
t = 0:1/Fs:0.5*(Fs-1)/Fs;
f1 = 30;
w = (1-2*pi^2*f1^2.*t.^2) .* exp(-(pi*f1.*t).^2);
X = fft(w);
plot(angle(X));

Your vector w is not symmetric, so don't expect a zero-phase FFT result. Just do plot(t,w) to see what the FFT sees as its input. The reason is that you computed w only for positive t. What you need to do is the following:

Choose an odd number of equidistant time domain points $2N+1$
Compute N+1 points of the Ricker wavelet for non-negative $t$ starting with $t_0=0$
Since the wavelet is symmetric, the $N$ points for $t<0$ are identical to the $N$ points for $t>0$
append these values for negative $t$ at the end of the time domain vector to be transformed: w = [r_0,r_1,...,r_N,r_N,r_{N-1},...,r_1]. The reason for doing this is that the DFT/FFT assumes the first vector entry to be the value for $t=0$. Appending the values for negative $t$ at the end simply reflects the implied periodicity of all signals transformed by the DFT.

The result will be a real-valued FFT up to numerical errors.

In Matlab this would look like this:

Fs = 250;
N = 50;
t = (0:50)/Fs;
f1 = 30;
w = (1-2*pi^2*f1^2.*t.^2) .* exp(-(pi*f1.*t).^2);
w = [w, w(N+1:-1:2)];
X = fft(w);
max(abs(imag(X)))     % 2.2682e-16 on my machine


Problem identifying the analytic expression of such determined signal I came across this problem

I am supposed to find the Fourier transform of $g(t)$, but I am not able to find the analytical expression of such signal.

The teacher suggests that I should consider the inclined part in $\left[-T, 0\right]$ as the integral of a determined function.

It turns out it is the integral of some $\rm{rect}()$ function. However, the teacher noted that I should subtract from it the (???) step function.

Can you explain to me how to find the fourier transform of $g(t)$ more comprehensively?

EDIT 1:

I tried to solve it, and I came to this solution. Is it correct?

The simple answer is “Linearity”

The Fourier Transform is linear.  $g(t)$ can be expressed as a sum of signals that have known, as in there are tables of transform pairs, in most books, probably also in your course notes, along with the shift and scale properties. Look at the tables you have and find entries that you can add to get $g(t))$.  You already found rect().  Look for a triangle piece.

A full triangle is the convolution of 2 rects.  A half triangle is a full triangle multiplied by a step.


How can I move a signal in frequency so that it ends at the frequency zero? I am trying to move a signal to 0 in the frequency domain but I am not sure how to do it.

One way of moving a frequency response of signal is multiplying it with a sinusoid. Assume x(t) is your signal.Let the sinusoid be $$cos(w_0t)$$
If I multiply this to x(t) as $$y(t)=x(t)cos(w_0t)$$
Then Fourier transform of y(t) is $$Y(w)=0.5(X(w-w_0)+X(w+w_0))$$
This is because multiplying in time domain of two signals is equivalent to its convolution in frequency domain.Since fourier transform of the above sinusoid is$$C(w)=\frac{\delta(w-w_0)+\delta(w+w_0)}{2}$$Due to the sifting property of delta function we get the above Y(w) after multiplying.
Although here we get two signal the required one can be obtained by passing it through band pass filter


Understanding correlation in image processing using examples I have been trying to understand the intuition behind the correlation. I understand from this guide that correlation is a way to detect a known waveform in a noisy background. I have seen from this SO answer that correlation can be used for template matching. So I decided to try it and hence the reason for my confusion.

Consider a 2D array,

+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+
| 0 | 0 | 1 | 1 |
+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+

I did correlation without padding with the below array,

+---+---+
| 0 | 0 |
+---+---+
| 1 | 1 |
+---+---+

The output I got is,

+---+---+---+
| 2 | 2 | 2 |
+---+---+---+
| 0 | 1 | 2 |
+---+---+---+
| 2 | 2 | 2 |
+---+---+---+

I was expecting a maximum value where the smaller array was present at the big array. But I got multiple peak values. Can you help me understand why this happened and how to overcome it?

Thank you.

As it was already posted multiple times: The problem comes from an inaccurate definition of correlation in your application.

The Pearson correlation coefficient does require the data to be

centered, ie the mean must be subtracted
normalized, ie the data must be divided by the standard deviation

This centering and normalization must be done for the mask as well for each sub-matrix of your larger matrix.

In your example, you would end up with a correlation matrix as:

$\left(
\begin{matrix}
0.0&0.0&0.0\\
-1.0&-0.577&0.0\\
1.0&0.577&0.0
\end{matrix}
\right)
$

Lets have a look at the matrix entries:

0.0 They are there, since there are sub-matrices with all ones in them. The mean of such a matrix is 1 and therefore it is centered to the corresponding matrix with only zeros. Hence the correlation is zero.
1.0 is the correlation that you would like to see. The corresponding sub-matrix fits the mask perfectly, hence full correlation, i.e. 1
-1.0 is the correlation of the sub-matrix that fits the mask, but the sub-matrix entries are switched (i.e. the sub-matrix is 0 where the mask is 1 and vice versa). This is a perfect anti-correlation, i.e. -1
The 0.577 cases: They occur where the sub-matrix deviates from the mask by one number (i.e. there is a 1 where the mask has a 0 and vice versa). It also happens to be a symmetric case in your example, therefore there is a pair that correlates positively with one switch and one that correlates negatively with a corresponding switch.

Here is a pretty ugly brute-force code to illustrate it. It is neither optimized in any way, nor tested for any other case than yours. I hope it is somewhat correct ;)

import numpy as np
In [2]:

img = np.ones([4,4])
img[2,0:2] = 0.0
print(img)
[[1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [0. 0. 1. 1.]
 [1. 1. 1. 1.]]
In [3]:

mask = np.ones([2,2])
mask[0,0:2] = 0.0
print(mask)
[[0. 0.]
 [1. 1.]]
In [4]:

def simple_convolve(img, mask):
    max_row    = img.shape[0] - mask.shape[0] + 1
    max_col    = img.shape[1] - mask.shape[1] + 1

    output = np.zeros([max_row, max_col])

    for curr_row in range(0, max_row):
        for curr_col in range(0, max_col):
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    output[curr_row, curr_col] += img[curr_row + curr_mask_row, curr_col + curr_mask_col] * mask[curr_mask_row, curr_mask_col]
    return output

def better_convolve(img, mask):
    max_row    = img.shape[0] - mask.shape[0] + 1
    max_col    = img.shape[1] - mask.shape[1] + 1

    output = np.zeros([max_row, max_col])

    mask_mean  = mask.mean()
    mask_sigma = (((mask - mask_mean)**2).sum())**0.5
    mask = mask - mask_mean
    mask = mask/mask_sigma

    curr_img_mean  = 0.0
    curr_img_sigma = 0.0

    for curr_row in range(0, max_row):
        for curr_col in range(0, max_col):

            curr_img_mean  = 0.0
            curr_img_sigma = 0.0

            # Compute the mean value of the sub-image of img
            ctr = 0.0
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_mean += img[curr_row + curr_mask_row, curr_col + curr_mask_col]
            curr_img_mean = curr_img_mean / ctr

            # Compute the std value of the sub-image of img
            ctr = 0.0
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_sigma += (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)**2.0
            curr_img_sigma = curr_img_sigma**0.5

            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    x = (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)
                    if curr_img_sigma == 0.0:
                        x = 0.0
                    else:
                        x = x/curr_img_sigma

                    output[curr_row, curr_col] += x * mask[curr_mask_row, curr_mask_col]
    return output
In [5]:

simple_convolve(img, mask)
Out[5]:
array([[2., 2., 2.],
       [0., 1., 2.],
       [2., 2., 2.]])
In [6]:

better_convolve(img, mask)
Out[6]:
array([[ 0.        ,  0.        ,  0.        ],
       [-1.        , -0.57735027,  0.        ],
       [ 1.        ,  0.57735027,  0.        ]])


upsampling and downsampling operator I am trying to prove that for all $x,y \in l^{2}(\mathbb{R})$, it holds that $$\left\langle\left(\downarrow M\right)\left[x\right],y \right\rangle = \left\langle x,\left(\uparrow M\right)\left[y\right]  \right\rangle \text{.}$$
How should I proceed to prove this?

The left-hand side is

$$
\sum_{i \% M=0}\sum_{j\%M=0} x(i/M,j/M)y(i,j)
$$

where the sum can only be done over all i,j which are dividable by M (otherwise x(i/M)=0). Hence, by substitution $i=Ma, j=Mb$ we get

$$
=\sum_{a\in\mathbb{Z}}\sum_{b\in\mathbb{Z}}x(a,b)y(aM,bM)
$$

which equals the right-hand side.


What does librosa.load return? I'm new to signal processing and I'm working with librosa library quite a bit.

I would like to know what is the information that is returned by librosa.load function when I read a audio(.wav) file using it. Is it the instantaneous sound pressure in pa? or the just the instantaneous amplitude of the sound signal with no unit?
How can I use this information to calculate the SPL(no-weigting) of the signal ?

Is it the instantaneous sound pressure in pa?

How would a piece of software know about sound pressure? What it gets is the digitized voltage from a microphone, with zero information about the physical sensitivity of that microphone.

  or the just the instantaneous amplitude of the sound signal with no unit?

Exactly! :)

  How can I use this information to calculate the SPL(no-weigting) of the signal?

Calibrate your recording setup with a sound source of known power. That's the only way short of buying an expensive calibrated measurement microphone, amplifier, digitizer. But then, there will be little need for librosa.


Calculating SNR from Frequency Domain As I am newbie to signals, I have collected an appliance signal in frequency domain as follows:

(PSD- Power Spectral Density). So, I need to calculate the SNR (Signal to Noise Ratio) and BER (Bit Error Rate) caused by this signal in Powerline Calculation.
So, I have made these codes in Matlab, but what I am stuck with is calculating the SNR and BER.

bits=10000; %number of bit
b=randi([0,1],1,bits); % generate random [0,1]
t=0:1/30:1-1/30; % Time period allocated for the signal

%ASK Carrier Signals
carrier_signa_l= sin(2*pi*t);
E1=sum(carrier_signa_l.^2);
carrier_signa_l=carrier_signa_l/sqrt(E1); %unit energy
carrier_signal_0 =0 * sin(2*pi*t); % zeros for 0 bits in the carrier signal

%MODULATION
ask=[];
for i=1:bits
    if b(i)==1 % If bit = 1
        ask=[ask carrier_signa_l];
    else
        ask=[ask carrier_signal_0];
    end
end

If it is a sinusoidal signal, there will be peak (among the frequency bins) in the frequency spectrum corresponding to the tone's frequency.

Ratio of the magnitude of this peak to the sum of the magnitudes of all other bins (which are noise) correspond to Signal to Noise Ratio.

But when its a non sinusoidal signal (like the one in your plot) you have to consider the relevant band of the signal instead of a single peak. you can quantify it by specifying a frequency bin corresponding to the mojor frquency component in the signal plus some leakage (related to bandwidth) into the nearby bins.

Please take a look in the following matlab code.

N               = 8192; % FFT length
leak            = 50;
% considering a leakage of signal energy to 50 bins on either side of major freq component

fft_s           = fft(inptSignal,N); % analysing freq spectrum
abs_fft_s       = abs(fft_s);

[~,p]           = max(abs_fft_s(1:N/2));
% Finding the peak in the freq spectrum

sigpos          = [p-leak:p+leak N-p-leak:N-p+leak];
% finding the bins corresponding to the signal around the major peak
% including leakage

sig_pow         = sum(abs_fft_s(sigpos)); % signal power = sum of magnitudes of bins conrresponding to signal
abs_fft_s([sigpos]) = 0; % making all bins corresponding to signal zero:==> what ever that remains is noise
noise_pow       = sum(abs_fft_s); % sum of rest of componenents == noise power

SNR             = 10*log10(sig_pow/noise_pow);


Pronunciation research for app Hello does any one know or at least can lead me to the right way?

What does not change pronunciation of words when someone is speaking. The idea is to create a app that checks if word/sentence is pronounced correctly. To do so I need something in pronunciation that does not change at all or changes in small range.

For example someone reads a text. His voice is recorded. Now someone else, who can not read so well reads it too. Then the recordings are compared to eacother. Then software displays on the computer screen how wrongly the second person pronounced.

Please do not suggest any software that does it!

You need the right audio features, and the modern way to do this is with deep neural networks, which learn them automatically. I would use one built for speech recognition (it's basically the same problem), then model the users by the output of the top layers. Accents should manifest as clusters of the features, and languages as superclusters. I don't have much experience here but I'm making an educated guess.

Here is a relevant paper (Deep Neural Networks for Acoustic Modeling in Speech Recognition), and here is a relevant student project (Accent Recognition with Neural Network).


What if the Down sampling factor is not an integer? I have developed a function which down samples an input signal.Say I have an input signal with a sampling rate of 512 samples/sec and would like to down sample it 128 samples/sec, then the down sampling factor is discrete and it is 4. In this case, I believe that the input signal will successfully be down sampled to 128 samples/sec without any loss of information.But, what if we have to down sample a signal from a sampling rate of 511 samples/sec to 127 samples/sec,then the down sampling factor 4.0236. When I run my code for a down sampling non integer factor, I still get the output down sampled signal,but I am sure that, since I down sampled the signal by an non-integer factor, the output may not be correct and would have lost information.Any methods to overcome this issue or Is there a way to approximate it without damaging the quality of the down sampled signal?

My Understanding - I found a method after a while of researching. In this case we have to initially up sample the input sampling rate followed by down sampling, known as "Multirate conversion". So here is my understanding, Input Sampling rate = 511 b/s and Desired output sampling rate = 127 b/s. So, 127*5 = 635 and 635 -511 = 124. Therefore, i need to up sample 511 b/s by a factor of L = 5 to get 635 and then down sample this 635 by a factor M = 124 to get the desired down sampled output of 127 b/s. What do you have to say?. Please correct me if its wrong or illogical.

Resampling (by rational or even irrational ratios) can be done by low-pass filtering in conjunction with high quality interpolation of all the samples needed for the new rate, directly.  No two-step up-sampling following by downsampling is required (although that is one possible implementation for simple ratios).

A Sinc kernel, being a reconstruction formula, is the ideal interpolator for band-limited sampled waveforms.  In practice, a windowed Sinc interpolator is usually as good or better than the FIR or IIR filters used inside other resampling methods.

Also, information may be lost in going to any lower sample rate (even in your original 4:1 resampler) if the original sample data contains higher frequency spectrum that has to be removed by the anti-alias low pass filter (required to prevent aliasing at the lower sample rate).