Frequency response  Impulse response I am trying to calculate an actual frequency response from a desired frequency response.

fftSize = 128;
fs = 8000;

I have 2 input arrays: f[] and aDesired[], both of which are Real and length fftSize+1 (DC + Nyquist included). The f[] array is linearly spaced between 0 and fs/2 (0:4000). The aDesired[] array is either gains in dB or linear values.

I can easily change the length of f[] and aDesired[] as I have an interpolation function which gives any size outputs, but i thought this was a good number.

I would like to take an IFFT of ADesired[] to generate coefficients (with a configurable length), window and truncate the coeffs, and then FFT to show what the actual response would be, aActual[].

How can I take an IFFT/FFT of this data when the inputs are both Real? Or how do I prepare the data more?

Finally, it doesnt matter about the phase, if its linear or not it doesnt matter.

Is there maybe a library that could be of use with my input data?

Summary:

Input: f[129] : linearly spaced fc's
       aDesired[129] : Real values of linearly spaced gains at fc's (either in lin or log),

Processing:
       coeffs = IFFT(f, aDesired)
       coeff_new Truncate(window * coeff)
       aActual = FFT(f, coeff_new)

Output: aActual[]

Many thanks

As Jason R pointed out in his comment, your way of obtaining an actual frequency response is just one way of doing it. There are many other ways of determining an actual filter response. Any other frequency domain FIR filter design method could do the job, and with each method you would get a slightly different result. Having said that, let's have a look at the way you want to do it:

You have a real-valued desired response. This is the magnitude of the desired frequency response, because for causal (realizable) systems, the frequency response is always complex-valued. However, this is no problem since any finite length impulse response can always be made causal by adding some delay. First, we need to make the desired response symmetric, as is necessary for real-valued systems. Let the column vector D be the vector of desired magnitude values on the equidistant frequency grid between 0 and $f_s/2$, then the symmetric desired response is (using Matlab/Octave syntax):

Dsym=[D;D(N/2:-1:2)];

where $N$ is the FFT length. We can obtain the (non-causal) impulse response by doing an inverse FFT:

h=real(ifft(Dsym));

Note that taking the real part is only necessary to remove the tiny imaginary part due to numerical errors. To make the impulse response causal, we need to shift the second half of it to the beginning:

h_c=[h(N/2+2:N); h(1:N/2+1)];

The impulse response can now be multiplied with a window function:

h_w = h_c.*w;

Your actual frequency response is then

H=fft(h_w);

It is complex valued, so in order to compare it to the desired response, you need to plot its magnitude:

plot(abs(H))

The phase of the actual response is linear, i.e. the system just adds some delay but does not introduce any phase distortions.


Can anybody recommend a good source for material of Gaussian and Complex Gaussian wavelet? Can anybody recommend a good source for material of Gaussian and Complex Gaussian wavelet?
Please share the link if possible.

The notion of "Gaussian wavelet", whether real or complex, is not clear from the literature. Indeed, the Gaussian has a positive average, while wavelets are (or should be) zero-mean (to satisfy the admissibility condition).

Two standard wavelets may suit your need:

the Morlet wavelet, or Gabor wavelet, a complex exponential  multiplied by a Gaussian window envelope (or its real part),
derivatives of Gaussians (including  differences of Gaussians or DoG), sometimes called Ricker wavelet in 1D in geosciences or mexican hat in 2D.

Their complex extensions, or real restrictions, are relatively natural.


Is it required to quantize wavelet coefficients before one derives features from? I am using wavelet coefficients for feature extraction in classification problem. As the wavelet coefficient values are real, positive and negative, is it required to quantize them before feature extraction?

The global answer is no, it is not mandatory. Yet weak forms of quantization (of values and locations) may help, like absolute values and range compression.

First, wavelet coefficient values can be complex. Wavelet  transform are usually linear, and are often shift-variant.  If you want a classification that does not depend on data sign or phase, and on (some) offsets, ie if you want the same class for:

black object on white or white object on black.
slightly offset signals

at least you need to take an absolute value, one of the weakest form of quantization.

Then, depending on what classification you are used, it is not mandatory. You can for instance fit a distribution to wavelet coefficients in each subband, and use the distribution parameters for clustering, without further quantization. For instance, generalized Gaussian-Laplacian distributions were once considered a good fit for wavelet coefficients, see for instance Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback–Leibler Distance. The following plot is taken from a paper related to adaptive filtering, however this generalized Gaussian modeling from wavelet frame decomposition with different power laws may apply to classification too.

More generally, range compression, like applying a power-law (Box-Cox, Anscombe transformations) to absolute wavelet coefficients, is often used to set thresholds, that are another form of classification.

The last quantization we often do not think about is location. We previously talked about values only. But you can gather wavelet coefficients in vectors or blocks, packed into a "Reference Observation Vector", and compute statistics on them, that will be used for further processing.


Does the Fourier series coefficient of AC components remains same if DC component is subtracted form the given signal? Suppose a signal is defined by
$
x(t)= \begin{cases}
t & 0\leq t \leq 1 \\
2-t & 1\leq t\leq 2 \\
\end{cases}
$

Since $x(t)$ has even symmetry, I can calculate fourier coefficient as
$$
a_n = \frac{4}{T} \int_0^1 x(t).\cos{n\pi t}.{dx}
$$
I have calculated $$\begin{equation}
a_n = 2\big[\frac{\cos{n\pi} - 1}{n^{2}{\pi}^{2}}\big]\tag{1}
\end{equation}$$

The DC value of $x(t)$ i.e $a_0 = 0.5$. If we subtract DC value we get,

From this we can see that given signal has hidden half wave symmetry in addition to Even symmetry. So we can find fourier coefficient as

$$
a^{'}_n=\frac{8}{T}\int_{0}^{\frac{1}{2}}(t-\frac{1}{2})\cos{n\pi t}.dt
$$

I have calculated $$\begin{equation}
a^{'}_n = 4\big[\frac{\cos{\frac{n\pi}{2}} - 1}{n^{2}{\pi}^{2}}\big]\tag{2}
\end{equation}$$

My question is, shouldn't $a_n$ and $a^{'}_n$ be equal for $n\neq0$ ?

PART-I: I would like to provide the general proof considering the title of the question and imposing no specific properties on the signal $x(t)$ other than having a CTFS representation.

The following is a simple analysis to conclude that the CTFS coefficients of any signal $x(t)$ and that of the DC removed signal are equivalent. (except $a_0$ of course).

Consider a continuous-time periodic signal $x(t)$ with period $T$ divided into two components: $x_{dc}$ and $x_{ac}$, with periods $T$ also, where $x_{dc}$ is the pure DC component of $x(t)$ and $x_{ac}$ is the pure AC component of $x(t)$, then we have:

$$ x(t) = x_{dc}(t) + x_{ac}(t)$$

Computing the CTFS coefficient $a_k$ of $x(t)$ yields:
$$
\begin{align}
a_k &= \frac{1}{T} \int_{<T>} (x_{dc} + x_{ac}) e^{-j k \frac{2\pi}{T} t } dt  \\
&= \frac{1}{T} \int_{<T>}  x_{dc} e^{-j k \frac{2\pi}{T} t } + \frac{1}{T} \int_{<T>}  x_{ac} e^{-j k \frac{2\pi}{T} t }  \\
a_k &= b_k + c_k \\
\end{align}
$$

where $b_k$ and the $c_k$ are the CTFS coefficients of DC and AC parts of $x(t)$.

By definition of any DC signal, it's known that $b_k = 0$ for all $k \neq 0$ and by defition of any AC signal it's known that $c_0 = 0$. Then using the relation $a_k = b_k + c_k$ we get the following:

$$a_0 = b_0 + c_0 = b_0$$

and
$$a_k = 0 + c_k = c_k ~~~,~~~ \text{ for all } k \neq 0 $$

From which we define :

$$
b_k = \begin{cases} a_0 ~~~&, ~~~\text{ for } k = 0 \\ 0 ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases}
$$

and

$$
c_k = \begin{cases} 0 ~~~&, ~~~\text{ for } k = 0 \\ a_k ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases}
$$

Hence we conclude that the CTFS coefficients, $a_k$, of any periodic signal $x(t)$ and the CTFS coefficients $c_k$ of DC-removed part, $x_{ac}$, are the same for all $k \neq 0$.

PART-II: Based on OP comments, the relation for an even and real signal, is the following.

For a signal $x(t)$ which is real and even we have
$$x(t) = x(t)^{*} = x(-t) = x(-t)^{*}$$
and the associated CTFS coefficients has the property of
$$ a_k = a_{-k}^{*} = a_{-k} = a_k^{*} $$
which indicates that the coefficients $a_k$ are also real and even.

Using this, we can obtain the trigonometric (cosine) Fourier series coefficients as.
$$
\begin{align}
a_k &= \frac{1}{T} \int_{<T>} x(t) e^{-j k \frac{2\pi}{T} t } dt  \\
&= \frac{1}{T} \int_{<T>}  x(t) \left( \cos( k \frac{2\pi}{T} t) + j \sin(k \frac{2\pi}{T} t) \right) dt \\
&= \frac{1}{T} \int_{<T>}  x(t)  \cos( k \frac{2\pi}{T} t) dt + j \frac{1}{T} \int_{<T>}  x(t)  \sin( k \frac{2\pi}{T} t) dt\\
a_k &= \mathcal{Re}\{a_k\} + j ~~ \mathcal{Im}\{a_k\} \\
\end{align}
$$

Now since the property states tat $a_k$ are real, then the imaginary part is zero and we have:

$$
a_k= \frac{1}{T} \int_{-T/2}^{T/2}  x(t)  \cos( k \frac{2\pi}{T} t) dt
$$

Furtermore since $x(t)$ is also even;  $x(t)=x(-t)$, then we also have

$$
\boxed{ a_k= \frac{2}{T} \int_{0}^{T/2}  x(t)  \cos( k \frac{2\pi}{T} t) dt }
$$

As the trigonometric cosine series coefficients of the real and even signals.

In addition to this, for a real & even signal, $x(t)$ of period $T$, which has no DC part, the following is also observed:
$$ x(t-\frac{T}{2}) = -x(t) $$

And based on the time-shift property of CTFS we can conclude that
$$ a_k ~~e^{-j\frac{2\pi}{T}k \frac{T}{2} } = - a_k $$

$$ a_k ~ e^{-j\pi k } = - a_k \implies a_k = \begin{cases} -a_k &, k=0,\pm 2, \pm4,... \\ a_k &, k=\pm 1,\pm 3,...\\ \end{cases}$$

Which indicates that the CTFS coefficients $a_k = 0$ for $k=2m$ (even) for a real, even, (and having no DC) signal $x(t)$. Indeed we can get rid of the DC removed adjective and state for all real & even signals, as the DC will only affect $a_0$ being non-zero.

PART-III: Finally apply these to the example signal to see that it works.

The signal  defined as:
$$
x(t) =
\begin{cases} t &, 0<t<1 \\ 2-t &,1<t<2 \\ \end{cases}
$$

Then for the CTFS coefficeints (in the trigonometric form) we have:

$$a_k = \frac{2}{T} \int_{0}^{1} t \cos(\frac{2\pi}{T} k t) dt = \int_{0}^{1} t \cos(k \pi t) dt $$

A simple by-parts integration yields the following result:
$$\begin{align}
a_k &=  \int_{0}^{1} t \cos(k \pi t) dt \\
&=  \frac{t ~\sin(\pi k t) }{\pi k} |_0^1 -\int_{0}^{1} \frac{\sin(k \pi t)}{\pi k} dt \\
&=  \frac{\sin(\pi k) }{\pi k}  + \frac{1}{\pi k} \left( \frac{ \cos(\pi k) - 1}{\pi k} \right)\\
\end{align}
$$
we conclude that

$$\boxed{ a_k = \frac{ \pi k ~\sin(\pi k) + \cos(\pi k) - 1}{ \pi^2 k^2 } }$$

Form which it can also be seen that

$$a_k = \begin{cases} 0.5 &, k=0 \\ \frac{-2}{\pi^2k^2} &, k=\pm 1, \pm 3,... \\ 0 &, k=\pm 2, \pm 4,...\\ \end{cases} $$

Note that the term $\sin(\pi k)$ can be ignored (as it's all zero) except for $k=0$. Also note that those zero values $a_k$ turns out to be the case without assuming half wave symmetry.

Finally, we shall compute the CTFS $c_k$ coefficeints of the DC removed signal $x_{ac}(t)$ to see if they are equivalent. From the definition of the signal we see that

$$
x_{ac}(t) = x(t) - 0.5 =
\begin{cases} t-0.5 &, 0<t<1 \\ 1.5-t &,1<t<2 \\ \end{cases}
$$

then the $c_k$ become:
$$c_k = \frac{2}{T} \int_{0}^{1} (t-0.5) \cos(\frac{2\pi}{T} k t) dt = \int_{0}^{1} t \cos(k \pi t) dt - \int_{0}^{1} 0.5 \cos(k \pi t) dt $$

Note that this integral is the same for the case of $a_k$. The only diference is in the last term which is $1$ for $k=0$ and $0$ for all $k\neq 0$ and can be ignored for $k \neq 0$. Then the equation relating to $c_k$ will be identical to that of $a_k$ except at $k=0$ which yields:

$$c_k = \begin{cases} 0 &, k=0 \\ \frac{-2}{\pi^2k^2} &, k=\pm 1, \pm 3,... \\ 0 &, k=\pm 2, \pm 4,...\\ \end{cases} $$

hence we again concluded that
$$
\boxed{ c_k = \begin{cases} 0 ~~~&, ~~~\text{ for } k = 0 \\ a_k ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases} }
$$


Codec for forensic video Security camera and recorder vendors tout HVEC/H265 codec support in their systems.  As far as I can tell, this is the industry top-of-the-line.  I’ve been using such a system for a few years now and I am continually frustrated by the (I believe) characteristics of the HVEC for this application of forensic video limited by the available settings in the (also terrible) management software.  It seems to me that camera sensors are capable of VERY crisp image capture, but that the lossy encoding makes most video recordings useless for moving faces and plates at reasonable distances.  In security camera footage, almost all content is redundant and uninteresting.  When the scene is still, it seems I could zoom in and count the blades of grass across the street.  I can clearly read distant street signs or the license plates of parked cars.  But when a vandal of foot smashes a car window, or a drunk driver sideswipes a parked car, all I can recover is a useless lossy blur.  It’s very frustrating to have clear boring footage and blurry action footage.

The camera sensors MIGHT be to blame if the issue is motion blur because of shutter speed, but it doesn’t seem to me to be the main problem.  I’ve tried a variety of CBR/VBR, resolutions, framerates to try to get better results.  The software enforces bandwidth limits per camera feed and among all feeds, but I wonder if they are too conservative or if the real limit is the realtime processing power needed to encode.

HVEC was developed for entertainment applications mostly, wasn’t it?  Where subjects of interest are much closer to the camera, we get familiar with the characters, and lossy-related motion blur is more acceptable.

Is there is a better codec choice out there for forensic video?  Is there a codec in which still regions are de-prioritized and motion or even prerecognized areas of interest, such as faces, figures, letters or areas of just high contrast edges, are prioritized in the encoding process?  Is this all possible in HVEC but just apparently not implemented in the commercial-professional recommended systems that I’ve seen?

Don't discount the contribution of physics to the blur in your pictures: Especially in low light, you need to have long exposure times, to get a picture that isn't completely dominated by grain. Conversely, things become less sharp if they're moving.

In fact, security cameras tend to have low frame rates (that might actually be for historical reasons – because nobody actually cared about the lawn growing in 25 fps, or about getting 100 grainy shots of someone smashing a window in 4 seconds, instead of 10 shots that have better exposure). So, that might actually be the limiting factor here.

  HVEC was developed for entertainment applications mostly, wasn’t it?

Certainly, the highest money flow would be from the movies industry, but they have the same problem, namely that action scenes are harder to compress than still live shots of landscape. So, HEVC tries very hard to be better at accomodating different types of video content, so, I'd assume it's the best standard solution.

All in all, HEVC does a lot of things to make action scenes work best – motion prediction on a sub-pixel level, high dynamic range, you name it.

Especially, it compresses still-standing parts very efficiently, leaving more headroom on storage space when there's movement. Of course, the encoder has to be configured to actually spend the bitrate exactly when necessary. Maybe your vendors misconfigure it. Impossible to tell. Maybe look at the instantaneous bitrate over time; if it doesn't change much when there's something appearing in the picture, the encoder was configured to aim for a constant bitrate rather than a constant quality.

Of course, there are numerous approaches to making surveillance video capture more efficient. (I wouldn't be surprised if a lot of that wasn't available for civilian purposes, though. Automated video surveillance is very much something that states do.) I know that several "smart" home cameras draw rectangles around areas they consider noteworthy based on some more or less artificially intelligent algorithm.

Source: Pyimagesearch/Adrian Rosebrock

Of course, if you attach enough neural network processing power to your camera, and train accordingly, you can make sure your system reliably doesn't even store background, but recognizes e.g. cats and people in the picture.

Source: Micah Scott (Scanlime): Tucoflyer Smart DNN-based cat tracking indoor camera gimbal system

If you know what you don't care about and that you care about the rest, and only record the latter, you can encode more efficiently by not even trying to record your grass growing. But: this takes a solid amount of processing power and hence electricity.

Taking that thought further:

Video compression is literally nothing but taking a high-dimensional input and mapping it to a lower-dimensional description and quantizing that cleverly.

That is also a pretty good description of what a video object classifier neural network does. Thus, there's plenty of approaches to encode video through neural networks.

You might then store video regions of interest not even as "collection of picture information", but simply as coefficients necessary to recreate e.g. the same face from a neural network optimized at representing faces by feature vectors with little error. This works pretty well. In fact, so well, that by synthesizing such feature vectors, we can create convincingly real-looking faces. For example This Person Does Not Exist:

It seems at first that the quality as evidence of a face picture recreated from a camera neural network, which you cannot prove to be right, is very limited. Or at least, should be, in a court of law.

  My algorithm says this is the visual representation of the lowest-error feature vector of a face that would match what the camera saw at 3 a.m. when someone broke into my car. I cannot formally prove an upper bound for the error, but it looks just like the defendant. So, it's the word of my trained-on-thousands-of-pictures neural network against the word of the defendant! If that doesn't prove him guilty, what will?

shouldn't be an all too successful argument when in dubito pro res is a principle in criminal justice.

Then again, all video compression claims the same; got a JPEG picture of the thief breaking into your store done with your iphone? These JPEG coefficients just are the lucky coincidence of image content, shot noise and the JPEG encoder.


Preprocessing to improve Harris corner tracking between video frames? In our unmanned aerial vehicle grayscale video image stabilization application, we're having difficulty finding the "good" Harris corners in frame N+1 selected from frame N. The source of the difficulty appears to be radical nonuniform pixel contrast changes between frames. Perhaps being airborne and using a slow frame rate (~3 fps) is the root cause of the pixel contrast shift.

We've tried various histogram equalization techniques to try to smooth out the pixel contrasts between frames in order to improve tracking "good" Harris corners in frame N+1. Results continue to be poor.

Does anyone have any suggestions on how to improve Harris corner tracking between video frames in a slow frame rate daytime aerial environment? Thank you kindly in advance.

Edit: 30 Jan 2012, added test case (not actual frame size) images

Summary Update: 8 Feb 2012. People suggest Harris corners are not so useful in grayscale video feature tracking. Answers below suggest and provide links to various alternatives. We are evaluating these alternatives and I will report results when we get to that point. Thank you, all, for your comments and answers.

Here's the previous frame N with 35 "good" 5x5 harris corners selected. The original frame is 8 bpp raw pixels.

Here's an excellent 5x5 Harris corner located at row 59 col 266:

The current frame N+1 with a few tracked 5x5 Harris corners, only one of which is valid:

The previous frame 5x5 Harris corner appearing in frame N+1 at r47 c145:

Note how the pixel intensities in the selected 5x5 have all changed in a nonuniform manner from the previous frame to the current frame. Contrast equalization techniques between frames do not help in detecting previous frame selected 5x5 pixels in the current frame. All suggestions welcome.

Can you try a different feature detector? FAST may be, erm, faster, and a higher frame rate will make matching easier (assuming your features are moving a lot between frames)

Looks like you are trying to use the grayscale region around the identified feature point to match from frame to frame.  This is likely to be poor, especially if there is lots of movement between frames.

You may get better performance by using what is often called a "descriptor" of the region around the feature point on which to match.

Some possible descriptors are used within the SURF and SIFT algorithms - they also have techniques for identifying regions to track, but you've bypassed the requirement for that.  A simpler descriptor designed for use in stereo-vision matching is the Census transform although that may not work as well for matching from frame to frame, for much the same reason as your current method (which is also widely deployed in stereo-vision)

The book to read on this is Multiple View Geometry in Computer Vision.


Performance evaluation in terms of channel MSE and BER The first figure on the left (in this paper) provides the channel NMSE performance and the the second figure from left provides detector's BER performance.

I am curious, when the channel MSE overlapped in the first figure, why there is a difference between the BER performance.

I think, that there should be no difference, Moreover, even if there is a gap between MSE the use of quantizer of quantizer can bridge the gap in BER performance.

Thanks

From the simulation section of the paper:

  The SER for each Monte Carlo simulation is calculated as follows: only
  when the activity state and the transmitted symbol are both detected
  correctly for a certain user, the corresponding user’s information is
  considered to be correct, and all symbols transmitted by the user
  whose activity state is mistakenly detected are considered to be
  incorrect. To provide a benchmark for performance comparison, we
  consider the oracle LS by assuming the actual active user support set
  known at the BS.

  However, the SER performance gap of the proposed joint CE and MUD
  framework and the oracle LS always remains. This is mainly due to the
  CS-based sparse signal reconstruction algorithm with non-orthogonal
  training matrix suffers significantly, so that the data detection
  error is the main cost of user activity detection.

The non-oracle receiver has to figure out who is transmitting and hence can do additional errors. The oracle knows this information already. That's why you have the difference in the SER.


How to put specific frequency in a specific bin in frequency domain? I have eight tones as represented by f, and I want to the distance between every two tones is different as shown in the figure,
I have the frequency spacing, but I need to put the tones properly, and I need the resolution and the number of samples to be suitable for FFT consideration.

N=8;
r=1;
f_d=1e3;
M=4;
m=0:M-1;
fc=2.45e9;
%% freq spacing
for j=1:M
f_spac(j)=(1+m(j)*r)*f_d;      %freq spacing
end
%% Generate the frequency vector
for i=1:M
    f(:,i)=[fc+((1+m(i)*r)/2)*f_d; fc-((1+m(i)*r)/2)*f_d];
end
B = reshape(f,[1,N]);
fs=20*f(N/2); % Sampling frequency
Ts=1/fs; % Sampling time
t=0:Ts:1e3*Ts; % Time vector

I ran your code, there were some errors. f_spac variable is not used after it is filled inside the for-loop. What you can do is to define a variable f_delta=1000. This will make sure the 'f' variable is filled properly with spacing of 1000 and its multiples.

N=8;
r=1;
f_d=1e3;
M=4;
m=0:M-1;
fc=0;
%% freq spacing
for j=1:M
   f_spac(j)=(1+m(j)*r)*f_d;      %freq spacing
end
f_delta = 1000;
%% Generate the frequency vector
for i=1:M
    f(:,i)=[fc+((1+m(i)*r)/2)*f_delta; fc-((1+m(i)*r)/2)*f_delta];
end
B = reshape(f,[1,N]);
fs=20*f(N/2); % Sampling frequency
Ts=1/fs; % Sampling time
t=0:Ts:1e3*Ts; % Time vector

At the bare minimum, you need FFT size to be 8 with Fs = 4kHz. But before that you need to down-convert the signal to baseband.


How to compare vocal words in matlab We as a group are new to signal processing and speech recognition. We have stored .wav sounds of three different words and we are trying to record a new .wav sound and have the system compare the sound to these three pre-recorded sounds to see which one correlates the most to the newly generated sound.

We have altered the .wav sounds in various ways trying to figure out how we should tackle this problem. Thus far, we have been stuck on this issue for a while now and are only able to achieve sound comparisons with just two different sounds.

How should we tackle this problem? We would be very grateful for some general directions, thanks in advance.

Dan Ellis has some neat Matlab scripts that allow you to pull audio features from files . . .

http://labrosa.ee.columbia.edu/matlab/

The dynamic time warp might be a good starting point for the task you describe.


python; audio; detecting silence in audio signal I have real-time audio application written in python, and I have question about detecting silence.
What's the most efficient way to distinguish silence from non-silence, in an input audio signal? By non-silence I mean any sound, even noise. I need to know when silence-period is in the process, as accurately as possible. I don't necessarily ask for some code-algorithm; theoretical tips and explanations would be also useful. I know that silence have low level of oscillations (signal is flat); in numpy.fft.fft output, silence has highest amplitude in [1] bin, and so on, there are several significant things when silence is there.
But using any of these conditions in code,  cannot give me absolutely accurate result to  detect silence and completely split it from sound. Any advanced techniques for this task?

i'm not gonna deal with python nor numpy nor any other parochial computational platform.

first, since "silence" is a perceptual property, you need to apply a weighting filter, such as A-weighting to boost the frequency components of the audio that our ears are more sensitive to and attenuate the portions we're less sensitive to.  there are other weighting curves besides A-weighting.  i like Robert Wannamaker's E-weighting the best.
then square the filtered signal (to get power or energy).
then low-pass filter the squared and filtered signal (to get an envelope).
then compare this envelope to a threshold that will correspond to the threshold of hearing.
lastly, you may need to put in a little hysteresis in the threshold comparison so that the silence/no-silence result does not jitter as your envelope crosses the threshold.

this is similar to the "gate" function in a compressor/limiter/gate process that sometimes finds expression in audio applications.  the "gate" was to make sure that something that's perceived as silence (or sufficiently low-level that it should be perceived as silence) becomes zero during the "silent" periods of time.


complex multiplier in divide and combine FFT I am studying radix 2 algorith from Proakis' book.

But I'm a bit confusied why 1st DFT $G_1$ is not multiplied by complex entity while 2nd DFT $G_2$ is being multiplied by complex entity $W$ as shown highlighted in attached figure.

Now as the document says, let $x[n]$ be a length $N$ (even) sequence whose even and odd indexed samples are denoted as $g_1[n]$ and $g_2[n]$ of length $N/2$ each.

Then the $N$-point DFT $X[k]$ of $x[n]$ can be written as folows:

$$\begin{align}
X[k]  &=  \sum_{n=0}^{N-1} x[n] e^{ -j \frac{2\pi}{N} n k} ~~~,~~ k = 0,1,...,N-1 \\\\
&=  \sum_{n=0,2,4}^{N-1} x[n] e^{ -j \frac{2\pi}{N} n k} + \sum_{n=1,3,5}^{N-1} x[n] e^{ -j \frac{2\pi}{N} n k}\\\\
&=  \sum_{m=0}^{N/2-1} x[2m] e^{ -j \frac{2\pi}{N} 2m k} + \sum_{m=0}^{N/2-1} x[2m+1] e^{ -j \frac{2\pi}{N} (2m+1) k}\\\\
&=  \sum_{m=0}^{N/2-1} g_1[m] e^{ -j \frac{2\pi}{N/2} m k} + e^{ -j \frac{2\pi}{N} k} \sum_{m=1}^{N/2-1} g_2[m] e^{ -j \frac{2\pi}{N/2} m k}  \\\\
&=  \left( \sum_{m=0}^{N/2-1} g_1[m] e^{ -j \frac{2\pi}{N/2} m k} \right) + e^{ -j \frac{2\pi}{N} k} \left( \sum_{m=1}^{N/2-1} g_2[m] e^{ -j \frac{2\pi}{N/2} m k} \right) \\\\
&=  G_1[k] + W_N^k G_2[k] ~~~,~~ k = 0,1,2...,N-1 \\
\end{align}$$

Where we have recognized the summations inside the parenthesis as the $N/2$ point DFTs of the sequences $g_1[n]$ and $g_2[n]$ respectively. Note that when $k$ spans the range $k=0,1,...,N-1$, the $N/2$ point DFTs will repeat twice.


Verlet integration first step I'm trying to implement color-to-grayscale method from this paper. And they use Verlet's integration as:

$$L^*(t+\Delta t)=\frac{F(t)}{m}\Delta t^2+2L^*(t)-L^*(t-\Delta t),$$

for computing new value of luminance. But they not say how to compute new value of luminance for $t=0$. According to this Wikipedia link I think it could be as:

$$L^*(\Delta t)=\frac{F(0)}{2m}\Delta t^2+L^*(0)$$

Is it correct?

The function varies as you increase the time step ∆t , t will be zero from the start, So at ∆t =0, the initial condition L*(0) is the gray scale value obtained by converting the color image into L* a* b*  color space.


Digital Filters deal only with samples, right? I'm new to digital filters. So I'm trying to get things right and I can't find an explicit answer to my question on the internet.

Question: Digital filters only accept samples as input? I mean the input can not be zeros nor ones. It can only be samples. If that's the case, then the typical position for a digital filter would be after the receiver (I mean before converting to bits).

So receiver would be like:

$$
\mbox{antenna} \rightarrow \mbox{amplifier} \rightarrow \mbox{Sampler} \rightarrow \mbox{Digital Filter} \rightarrow  \underbrace{\mbox{DAC}}_{\tiny \mbox{here I should have filtered samples}}
$$

Please correct me and kindly ignore any coding or LNAs and stuff like that.

You are almost right: digital filters do deal with samples, but a sample can be any numerical representation of a given signal value at a given instant (so in general, they may accept zeros or ones).

Moreover, a sample is usually represented by a binary word (e.g. 0001), so a digital filter actually deals with 0s and 1s.


What bandpass filter design will yield the shortest impulse response? Designing a simple 2nd order IIR bandpass Buterworth filter with a centre frequency of 500 Hz and a 1 octave bandwidth gives me the following frequency response ...

Now, if I take the impulse response and normalise it and convert it to dB, we can observe the decay of the impulse response.

The decay of the impulse response is approximately linear with time when plotted on this scale, allowing us to define a decay time statistic (just like in room acoustics where you can define reverb time). For the impulse response of this filter to drop below 30 dB, it takes about 11 ms.

We are trying to minimize this decay time keeping the following constant:

-3 dB bandwidth
Filter order

I am happy to accept (within limits) passband and stopband ripple, and/or a compromise on the steepness of the transition band to achieve this. Can anyone suggest a method for filtering with the shortest possible impulse response duration as defined above?

The Butterworth filter's frequency response is the result of specific formulas and its characteristic is the flat passband frequency response. Consequently, if the coefficients of the IIR filter are modified in any way, the filter might not maintain the "Butterworth" characteristics.

In addition to the responses by "Hilmar" and "Jason R", maybe you could treat this as an optimisation problem of a fitness function that captures your specifications.

You could for example start with some design (e.g. a Butterworth filter) and then use an optimisation technique to move the zeros and poles about their locations (or modifying the design by adding / removing poles and zeros) trying to achieve your specifications (a sharper time-domain roll-off maintaining bandwidth and filter order).

Along this line, a large amount of work has been performed on designing filters with Genetic Algorithms (and here) and Simulated Annealing (and here) which you might find useful.


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

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

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

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

Thanks for your answers!

Too many things to comment here:

1- Decibels and thousands never go together. -1000dB means $10^{100}$ power attenuation (or $10^{50}$ amplitude attenuation), which is A LOT (an astronomical lot!): the noise level of any electronics, the DAC quantization noise, or the noise introduced by digital arithmetic will make your 1000 dB attenuation useless. Furthermore, designing a filter for such attenuation will normally result in a order so high it is impossible to implement.

Consider maximum gains or attenuations in the 100s of dB. As a reference, signal to noise ratio of CD encoding is around 96 dB. Or a satellite radio link can have a free-space attenuation of hundreds of dB.

2- The $20 N$ dB/decade attenuation of (low-pass) filters is an asymptotic slope. The actual frequency response depends on the type of filter (Butterworth or Bessel have the smoothest response). Which to choose depends on your application.

3- How does your system work?

analog signal -> ADC -> digital LPF -> downsampling -> average

or

analog signal -> analog LPF -> ADC -> average

Why don't you place a big RC filter and sample at a lower frequency? You'll still get your average.

In any case, you always need a low-pass filter before the sampler to avoid aliasing.

4- for your (c) question: if you consider amplitude response of the filter, then it is $20 \log$. If you consider magnitude squared (or power response) then $10 \log$.


How to construct ideal low pass filter in matlab? Firstly, I was asked to construct a signal $x(t)= 2 \cos(2 \pi f_1 t) + 1.5 \sin(2 \pi f_2 t) $ with $f_1= 100, f_2=200, t= 50 \textrm{ms}, Fs= 1\textrm{kHz}$ and number of FFT points $N = 1024$ and my code until this part is below.

f1=100;
f2=200;
fs= 1000;
N= 1024;
t=0:1/fs: 0.05;

x= 2*cos(2*pi*f1*t) + 1.5*sin(2*pi*f2*t)

I don't know that if this is the ideal low pass filter with frequency $f_c$, how to construct it in MATLAB?

The ideal lowpass filter is an infinitely long sinc function. It's Fourier transform is a rectangular shape as shown in your frequency spectrum diagram. In practice you have to window (truncate) it to a certain number of samples. The periodic width of the lobes of the sinc will correspond to the width of your frequency rectangle (lowpass cutoff frequency)


NRZ vs. Manchester pulse shaping filter (Modulator) I'm trying to visualize the differences on the output waveform of a modulator when using NRZ or a Manchester pulse shaping filter. To keep it simple, let's think of a binary PAM signal with an input bit sequence 101011.

My thoughts are that with a Half Sine or a Raised Cosine pulse shaping filter, I will see the transitions between the 1's and 0's smooth out, but I'm struggling to think about the differences between the NRZ and the Manchester.

Manchester coding effectively creates a bit sequence that is twice as long as the given bit sequence and then applies NRZ pulses to the resulting bit stream.
There are two types of Manchester coding. In one form, each $0$ in the original
sequence is replaced by $01$ and each $1$ by $10$ while the other form complements
these bit patterns.  Using the first form,
$$101011 \to 10\,01\,10\,01\,10\,10$$ while the second form gives the complementary
bit sequence $01\,10\,01\,10\,0101$. If each data bit is of duration $T$, the NRZ
waveform with Half Sine pulses has sines of period $2T$ so that each bit
pulse has duration half a period. Thus, that leading $1$ gets
modulated to $\sin(\pi t/T)$ lasting from $t = 0$ till $t = T$.
With Manchester coding, each channel bit
(of duration $T/2$) gets a Half Sine of a sine of period $T$. Put another
way,

with plain NRZ, that first $1$ in your bit sequence will have
pulse $\sin(\pi t/T)$ lasting from $t = 0$ to $t = T$. This is a positive
going pulse rising from $0$ at $t=0$, peaking at $1$ at $t=T/2$, and
falling to $0$ at $t = T$
with Manchester coding, that data $1$ will become channel $10$ giving a positive pulse $\sin(2\pi t/T)$ from $0$ to $T/2$. The next channel bit
$0$ has the delayed half-sine pulse $\sin(2\pi (t-T/2)/T$ from $t=T/2$
to $t = T$, but with a negative amplitude (the channel bit
during this interval is $0$ instead of a $1$) which
works out to be
$$\begin{align*}
-\sin(2\pi(t-T/2)/T) &= -\sin(2\pi t/T - \pi)\\
&= -\sin(2\pi t/T)\cos(\pi) + \cos(2\pi t/T)\sin(\pi)\\
&= \sin(2\pi t/T)
\end{align*}$$
lasting from $t = T/2$ to $t = T$. In short, during the entire data
bit interval $[0,T]$ the signal can be expressed as $\sin(2\pi t/T)$
and is one cycle of the sine wave of period T. If the data
bit were a $0$ instead of a $1$, the channel bits would be $01$
and we would get $-\sin(2\pi t/T)$ from $t=0$ to $t=T$. In short,
Manchester coding with half-sine pulses
effectively creates a phase-modulated signal
from the data sequence with a carrier frequency of $1/T$ Hz and
each data bit interval contains exactly one period of the sinusoid.
If you look at the figure for Thomas's Manchester coding
in Wikipedia,
what we have is the $\pm$ rectangular pulses comprising each
data bit interval, which look like a hard-limited sinusoid
of period $T$, being replaced by the sinusoid.

  The OP says, "My thoughts are that with a Half Sine ....  pulse shaping
  filter, I will see the transitions between the 1's and 0's smooth out,..."
  but this does not happen, there will be abrupt reversals of carrier
  phase at those data bit boundaries where there is a transition
  from $0$ to $1$ or from $1$ to $0$, pretty much as happens
  in plain PSK.


How to discover the caracteristic frequency of the use of an app by someone Let's say I have an user that uses the application one time every monday exactly... That means once a week. If I do the dfft I can see a peak in this frequency

How can I analyse the general case for many users and with apparently random frequencies of use?

Is the fft the best choice?

I can use mat lab or excel to try something

This is the fft of an user of the app by day

I am assuming by it that I can say he uses the app mostly each 2.3 days (1/0.4219)... It is a coeherent conclusion? How can I get the information of all the other peaks?

EDIT
This is the signal I am making the dft

I think a Fourier based frequency spectrum analysis may not be the best (most suitable) approach in providing you the necessary and useful information about that random service usage data.

Instead I would suggest you to use statistical data analysis methods. These methods are the primary tools that are used to analyse user behaviour by collecting usage data and computing the resulting statistics upon them.

Creating the related random variables and estimating their associated probability density functions will allow you to formulate the average quantities such as the expected number of days (or hours) between two consequetive logins, expected day of week (or hour of day) of maximum usage, minimum usage etc, and therefore enable you to make optimizations of resource utilizations based upon them.


Finding the four corners of each quadrilateral shapes in an image I am working on some image processing application using Matlab. The image after segmentation will look (more or less) like the figure below.

Now I want to find the four corners of each shape to crop it (I can access each using shape using its pixels as they have the same label after segmentation). I though of using a corner detector like harris, but it gives me many additional corners (as the sides of the shapes are not really straight due to discretization. for example, a straight line in an image will look like the image below, and hence will introduce many corners.

In even worse cases, due to the imperfection of the image segmentation, the corners may actually be lost and look like a curved line as shown below, even though the image before segmentation had these clear 4 corners.

Another approach that I thought of is first to crop each shape, and then to use the line information (using a line detection algorithm) of the 4 sides and find the intersection of each adjacent sides. However, the line information, again due to imperfect segmentation, may not exist (as a single straight line per side, a side could be detected as multiple non-connected lines). Moreover, again due to imperfect segmentation, the 4 sides may have some artifacts which will hinder proper line detection.

Any suggestions?

I though of using a corner detector like harris, but it gives me many additional corners (as the sides of the shapes are not really straight due to discretization.

In that case you should keep only the strongest feature points.

  the corners may actually be lost and look like a curved line as shown
  below, even though the image before segmentation had these clear 4
  corners.

In that case, you can use scale space, first blur the images excessively, then use Harris, or directly use SIFT detector and keep only the features points of the scale, that match your desired points.


Global vector descriptor Usually, algorithms as SIFT, SURF and many others provides a set of $k$ keypoints and the associated descriptor in $d$ dimensions (for example, in SIFT each descriptor has $d$=128 dimensions).

So, in order to describe an image we need a matrix $k \times d$ ($k$ descriptor vectors, each one in $d$ dimensions). So far so good.

My question is: how can we describe an image through a single vector?

This could be really useful since we could save a lot of space and because certain algorithms (like LSH) requires a vector as input/query.

In some papers (for example this, section 6.5) this approach is described as "global descriptors".

Up to know, I found only this paper but it doesn't seem so accurate (and it's from 2009, not so new).

APPLICATION:

I'm trying to implement a memoization optimization to the following application: given an image (taken with a smartphone) containing a movie poster, return the internet movie database's score of that movie.

The optimization (so what I'm trying to implement) is about "remembering" (memoize) the input image and the relative result (IMD score) and for each new input image, search "the most similar one" and return the already computed result (avoid the database query) if a certain threshold is reached.

Two "similar images" in this context are two images pictures containing the same poster but from different prospective and/or places.

The actual approach is solving an approximate nearest neighbor through LSH in order to find the most similar image, but in order to do that I need to represent a query (an image) through a vector, and not through a matrix.

There are lots of different ways. You could simply put all the image pixel values into a vector with dimension = width x height. :) What you really want to know is how can we describe an image using a lower-dimensional representation and still caption all the important information.

The same approaches that describe an image patch, e.g. SIFT, apply. Just think of the entire image as the patch. However, the method very much depends on the information in the image and the purpose of the vector. For example, you could use some histograms of the gradients, subsample the frequency spectrum, take a certain number of wavelet coefficients, perform edge detection, use moments of the image, divide it into super-pixels and so on. You might need the representation to have different invariances such brightness, rotation, contrast, phase or scale. Basically any method that gives you a condensed version of the information in an image can be put into a vector.

Research into the right method to use for your particular application is probably in the literature somewhere. As for saving a lot of space, this is what compression is all about. e.g. JPEG discards the smallest coefficients of the DCT of blocks in the image to reduce space.


Unambiguously detecting (not correcting) between $t+1$ and $2t$ errors in a Reed-Solomon message I know that Reed-Solomon codes of the form $(N,N-2t)$ can correct up to $t$ errors, and detect up to $2t$ errors, but I've only seen the correction procedure: Berlekamp-Massey or Euclidean to get the error-locating polynomial $\Lambda(x)$, then the Chien search to find its roots, then the Forney algorithm to find the error coefficients.

If I am told that there are at most $2t$ errors, is there a way to unambiguously tell how many errors $e$ there are? Or at least distinguish between the case of $e \le t$ (correctable errors) and $t < e \le 2t$? Is it just a matter of determining the degree of $\Lambda(x)$? Somehow that doesn't seem right, because then finding the roots of $\Lambda(x)$ would tell me where they were, and then I could correct them.

I'm confused, because both the Berlekamp-Massey and Euclidean algorithms appear to find polynomials of degree $2t$ or less -- for Berlekamp-Massey it's mentioned in Blahut's book, and for the Euclidean algorithm it's $\Lambda(x)S(x) + Q(x)x^{2t} = \Omega(x)$ for known syndrome polynomial $S(x)$.

The Berlekamp-Massey algorithm and the extended Euclidean algorithm (both further extended with the Chien search and the Forney calculation of error values) for decoding BCH and RS codes (hereinafter referred to as the decoder) have the following characteristic that is not well understood by many people:

  If there is a codeword $\hat{\mathbf C}$ that differs from the received word $\mathbf R$ in $t$ or fewer locations, then the decoder finds $\hat{\mathbf C}$.

Note the complete absence of any mention of $\hat{\mathbf C}$ being the correct (i.e. transmitted) codeword $\mathbf C$.  Now, if $\mathbf R$ does happen to differ from $\mathbf C$ in $t$ or fewer locations, then the   decoder does indeed find $\mathbf C$. Note that since the minimum distance is $2t+1$, there is no competing codeword $\mathbf C^\prime$ that is also within distance $t$ of $\mathbf R$ to lead the decoder astray: there is at most one codeword within distance $t$ from any given $\mathbf R$ and when such a codeword exists, the decoder finds it. But just because the decoder finds a codeword, there is no guarantee that it is the correct codeword $\mathbf C$; it might be some $\hat{\mathbf C} \neq \mathbf C$ that happens to be at distance $\leq t$ from $\mathbf R$. (The decoder is following the sailor's philosophy of "When I'm not near the girl I love, I love the girl I'm near" with the addendum that on the high seas, no girls are near).  What it is correct to say is that (with the usual assumption that fewer errors are more likely to have occurred than more errors)  with high probability the $\hat{\mathbf C}$ that the decoder finds is indeed the transmitted codeword $\mathbf C$.

But the question is what happens when there have been more than $t$ errors. Can such an occurrence be unambiguously detected? The answer is NO for the reasons described above: the decoder may well find a $\hat{\mathbf C}$ that is different from the transmitted codeword $\mathbf C$ exactly as the decoder is designed to do because $\hat{\mathbf C}$ is at distance $\leq t$ from $\mathbf R$ and the poor schmuck of a system designer will have no idea that s/he has just bought a pig in a poke.  Once again, probability comes to the rescue. In the (fairly unlikely) event that more than $t$ errors have occurred, it is very much more often the case that the decoder fails to decode in several possible ways (which may or may not be detected by the system designer depending on how many safeguards have been built into the system) than that the decoder produces an incorrect codeword.

The Berlekamp-Massey algorithm or the extended Euclidean algorithm find an error-locator polynomial $\Lambda(z)$ of degree $e$ and an error-evaluator polynomial $\Omega(z)$ where the inverse roots of $\Lambda(z)$ (as found by the Chien search) are the $e$ error locations and the Forney error-evaluation formula gives the error values at these $e$ error locations. When more than $t$ errors have occurred, it is usually the case that $\Lambda(z)$ has fewer than $e$ roots and/or that the Forney error-evaluation formula returns a zero error value at a purported error location that the Chien search has found. Both of these conditions indicate that something is awry but testing for these conditions is expensive in hardware implementations and time-consuming in software implementations. The simplest test for detecting a decoder failure is to let the Chien search and Forney formula do their jobs, and the recompute the syndrome of whatever the decoder spits out.  If the recomputed syndrome is zero, then the decoder output is a valid codeword - whether it is $\mathbf C$ or an impostor $\hat{\mathbf C}$ is unknown - but if the recomputed syndrome is nonzero, then the system designer knows that the decoder output is not a valid codeword and that the decoder has actually failed to decode.  But again, such a test is expensive to implement and very often skipped by system designers who know about it, and, of course, not even thought about by those who don't know of the possibility of decoder failure that can be detected.

For a well-designed RS decoder system, it is the case that
$$P_{\text{decoder output correct}} \gg P_{\text{decoder failure}} \gg P_{\text{decoder output incorrect}}.$$
The decoding spheres (the set of vectors or possible received words $\mathbf R$ that are at distance at most $t$ from some codeword and thus lead to a successful decoding, whether correct or incorrect) have very small total volume and so it is far more likely than not that when more than $t$ errors have occurred, $\mathbf R$ is not in any of the decoding spheres, and the decoder will fail to decode and spit out garbage. Whether it is worthwhile detecting such failures is up to the system designer.
It is sometimes the case that $1-P_{\text{decoder output correct}}$ is a perfectly acceptable word error rate and so ignoring decoder failure is an option, but in high-performance systems where very low error rates are required, detecting decoder failure and asking for a
re-transmission of a codeword that was badly munged up in transmission allows the achievement of the desired low error rates at the expense of increased delay.


Extreme audio stretching This tool (Paul's Extreme Sound Stretch) is able to produce extreme timestrech (50x, or even 100000x). Any song can become an interesting sound texture or ambient music (listen to this example).

The technique used is mentionned as:

  "spectral smoothing" the sounds.

How does an extreme time stretch algorithm work?

(it is opensource, I have studied it a bit, but I haven't been able to see a big picture of the algorithm yet)

I thought about :

sound input ====>  STFT ===> duplicate the frames 50x ===> ISTFT ===> output

but I'm really not sure about this.

Here is how the algorithm works:

So the steps are:

Input -> apply some window (hann in the frame) -> FFT -> Get amplitude > randomize the phase -> IFFT -> overlap and add -> out

Do it in a loop based in the strech Factor.

Paul wrote a Python code for all steps in this picture.


Use of the Sobel Filter for Image Gradient [MATLAB] I am using the Sobel filter for an RGB image. I have found two different ways to do that and the results look a little different. What is the difference between these two methods?

Method 1

 sbl=fspecial('Sobel');
 sbl_red=imfilter(Image(:,:,1));

Method 2

 [magred,dirred]=imgradient(Image(:,:,1),'Sobel')

Is there an advantage over one of the two methods? I am trying to find a structure of the image that is parallel to the x axis.

Well, if you go through the documentation of fspecial() you'd see it returns the following filter:

h = fspecial('sobel') returns a 3-by-3 filter h that emphasizes horizontal edges using the smoothing effect by approximating a vertical gradient. To emphasize vertical edges, transpose the filter h'.

[ 1  2  1
  0  0  0
 -1 -2 -1 ]

Namely in the way you applied it above you only got the vertical element of the gradient.

While imgradient() returns the image of the magnitude of the Gradient (Which has Vertical and Horizontal elements) and its direction.


$\mathcal Z$-transform of an equation [Exam question]: Verifying the solution I'm studying for exams at the moment and I'm trying to reproduce a solution from my professor (I have the solutions). The following signal is given:

The excercise says:

  Calculate the Fourier transforms for i = 13 to 15 and N = 64.

The signal in the time domain looks like this (this is where I'm still on par with the solution):
$$
f(k) = \left( -\frac{1}{N}k + 1 \right)\left( \sigma(k) - \sigma(k-N) \right)
$$

The problem starts when trying to transform the whole thing to the $z$-domain. This is my take:
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}k\sigma(k-N) + \sigma(k) - \sigma(k-N)
$$
Corrected for delay:
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}(k-N)\sigma(k-N) + N\sigma(k-N) + \sigma(k) - \sigma(k-N)
$$

Transformed:
$$
F(z) = -\frac{z}{(z-1)^2}\frac{1}{N} + \frac{z}{(z-1)^2}\frac{1}{N}z^{-N} + \frac{z}{z-1}Nz^{-N} + \frac{z}{z-1} - \frac{z}{z-1}z^{-N}
$$

But the solution only gives me three summands and I can't figure out how he is able to simplify the rest of the equation (Being well aware that the solution given by the professor might be wrong). He states:

  \begin{align}
F(z) &= - \frac{1}{N}\frac{z}{(z-1)^2} + \frac{1}{N}\frac{z}{(z-1)^2}z^{-N} + \frac{z}{z-1}\\
F(z) &= \frac{z}{(z-1)^2}\left[-\frac{1}{N} + \frac{1}{N}z^{-N} + {z-1}\right]
\end{align}

What am I doing wrong / am I doing anything wrong?
Can you please sanity check me?

Your mistake is quite simple, originating from a subtle inaccuracy in the reformulation from here
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}k\sigma(k-N) + \sigma(k) - \sigma(k-N)
$$
to here:
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}(k-N)\sigma(k-N) + N\sigma(k-N) + \sigma(k) - \sigma(k-N).
$$

The 3rd term of the second equation is only $\sigma(k-N)$ because the $N$ is cancelled by the $1/N$ from the original term. Then, the (corrected) third term cancels with the last one. And then, transforming this to Z-domain gives exactly the result your professor has provided.


Deriving the Euler-Lagrange Equations for the Chan-Vese Model I am reading the Chan-Vese paper "An Active Contour Model without Edges" and IPOL tutorial "Chan-Vese Segmentation" by Pascal Getreuer.  The energy for this model is given by:
$$
E(\varphi) = {\mu \int_{\Omega} \delta(\varphi(x)) | \nabla \varphi(x)|dx} + \nu\int_{\Omega} H(\varphi(x)) dx \;+$$
$$\lambda_1\int_{\Omega}|f - c_1| H(\varphi(x))+\lambda_2\int_{\Omega}|f - c_2|(1 - H(\varphi(x))dx
$$
Above, $\mu, \nu, \lambda_1, \lambda_2$ are real parameters, $c_1, c_2$ are constants determined for segmentation, $-1 \leq \varphi(x) \leq 1$ is the level-set function in which $\varphi(x) = 0$ specifies the interface,  $f$ is the original image, $H$ is the heavyside function in 1 dimension centered at $0$ and $\delta = H'$.  The Euler-Lagrange equations can be derived by inspecting the following derivative:
$$
\left.\frac{d}{d\epsilon} E(\varphi(x) + \epsilon \eta(x))\right|_{t = 0}
$$
for some test function $\eta(x)$.  Writing the Euler-Lagrange Equations as a gradient descent both papers claim:
$$
\begin{cases}
\frac{d\varphi}{dt} = \delta(\varphi(x))\left[\mathrm{div}\left( \frac{\nabla \varphi}{\mu |\nabla \varphi|}\right) - \nu - \lambda_1(f- c_1)^2 + \lambda_2(f-c_2)^2\right]\\
\frac{\delta(\phi(x))}{|\nabla \varphi|} \nabla \varphi \cdot \mathbf{n} = 0 & x \in \partial \Omega.
\end{cases}
$$
My problem is coming from the differentiation of
$$
{\int_{\Omega} \delta(\varphi(x)) | \nabla \varphi(x)|dx}.
$$
in the variation of the energy.  It should look like:
$$
\int_{\Omega} \underbrace{\left(\frac{d}{d\epsilon} \delta( \varphi (x) + \epsilon \eta(x) )\right)|\nabla \varphi(x)|}_{(*)} + \underbrace{\delta(\varphi(x)) \left(\frac{d}{d\epsilon}|\nabla\varphi(x) + \epsilon \nabla\eta(x)|\right)}_{(**)}.
$$
Now, $(**)$ gives us the curvature/divergence term in the Euler-Lagrange equations as well as the boundary condition after applying Green's first identity.  However, I don't understand how one might argue in this case why $(*)$ goes to 0 or why, for practical considerations, we can just ignore it.

Thanks to Luminita Vese, who responded to this question via email. I will post the answer here.  Let $\varphi_\epsilon(x) = \varphi(x) + \epsilon\eta(x)$ for some test function $\eta(x)$.

\begin{eqnarray}
\int_{\Omega} \frac{d}{d\epsilon} \left[ \delta(\varphi_\epsilon(x)) |\nabla (\varphi_\epsilon(x))|\right]dx &=& \int_{\Omega} \left[\frac d{d\epsilon} \delta(\varphi_\epsilon(x))\right] |\nabla (\varphi(x))|+  \delta(\varphi(x)) \left[\frac{d}{d\epsilon}|\nabla (\varphi_\epsilon(x))|\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| + \delta(\varphi(x))\left[\frac{d}{d\epsilon} \sqrt{\left<\nabla\varphi_\epsilon(x), \nabla\varphi_\epsilon(x)\right>}\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| + \delta(\varphi(x))\left[ \frac{\frac{d}{d\epsilon}{\left<\nabla\varphi_\epsilon(x), \nabla \varphi_\epsilon(x)\right>}}{2|\nabla(\varphi(x)|}\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| \\
&\quad&+ \delta(\varphi(x))\left[ \frac{{\left<\nabla\varphi (x), \nabla \eta(x)\right>}}{|\nabla(\varphi(x)|}\right]\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| \\
&\quad&+ \frac{{\left<\delta(\varphi(x))\nabla\varphi (x), \nabla \eta(x)\right>}}{|\nabla(\varphi(x)|}dx \quad (*)\\
&=&\int_{\Omega}\delta'(\varphi(x))\eta(x)|\nabla \varphi(x)|\\
&\quad&-\delta(\varphi(x))\mathrm{div}\left(\frac{\nabla \varphi}{|\nabla \varphi|}\right) \eta(x) \\
&\quad&- \delta'(\varphi(x))\frac{|\nabla(\varphi(x))|^2}{|\nabla \varphi(x)|}\eta(x)dx\quad(**)\\
&=&-\int_{\Omega}\delta(\varphi(x))\mathrm{div}\left(\frac{\nabla \varphi}{|\nabla \varphi|}\right) \eta(x) dx
\end{eqnarray}
Above, I am writing $\left< \cdot , \cdot \right>$ to denote the vector inner product, From (*) to (**), Green's first identity was applied as well as the product rule for divergence.


What happens when linear phase is critical? I am studying the topic digital filters. And in the slides it said,"if linear phase is not critical, Infinite Impulse Response (IIR) filters yield a much smaller filter order for a given application". So what does it mean ?

Causal IIR (infinite impulse response) filters cannot have a linear phase response, so they can only be used if (exactly) linear phase is not required. Unlike IIR filters, FIR (finite impulse response) filters can have a linear phase response (and, consequently, a constant group delay). So in applications where a linear phase is required, we have to choose FIR filters.

A non-linear phase response causes signal dispersion because different frequency components experience a different delay when being filtered. E.g., in many audio applications a linear phase is not required (although in some it might be), whereas in digital communications, a non-linear phase would cause neighboring symbols to interfere with each other (inter-symbol interference, ISI), so a linear phase is important.

Note that it is possible to design IIR filters with an approximately linear phase in the pass band (as shown in this answer). For many practical applications this might be good enough, and it might save a lot of computational complexity and memory compared to a linear phase FIR filter with a comparable magnitude response.

As mentioned in MBaz's comment, more details on applications where a linear phase may or may not be required can be found in the answers to this question .


Wrong values calculating FFT with EEG Bands using Numpy First of all I have to say I am very new to these matters. I am trying to apply FFT algorithm to some values obtained by EEG bands.
I found the following code but I am unable to make it work correctly.
import numpy as np
fs = 200
#data = np.random.uniform(0, 100, 1024)  # 2 sec of data b/w 0.0-100.0

data = np.array([(-0.27,-0.09,0.16,0.01,0.67,0.65,-1.59,0.11,-0.36,0.69,0.94,0.43,-0.53,-0.57,0.98,-0.97),(-0.74,-0.56,0.76,0.58,1.79,2.22,-4.96,-0.25,-1.19,1.57,2.95,1.47,-1.00,-1.06,2.97,-3.43),
(-0.89,-1.50,1.76,1.72,2.09,4.04,-7.97,-1.30,-2.09,1.11,4.82,2.57,-0.01,0.05,4.64,-6.44),
(-0.84,-2.78,3.00,2.88,2.00,6.36,-11.46,-2.36,-3.22,-0.01,6.96,3.84,1.85,2.11,6.64,-10.31),
(-0.39,-4.27,4.26,4.57,1.09,8.77,-14.78,-4.09,-4.40,-1.87,9.03,5.12,4.86,5.31,8.76,-14.77)])

print("After data")
# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))
print("After fft_vales")
print(fft_vals)
# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)
print("After fft_freq")
print(fft_freq)
# Define EEG bands
eeg_bands = {'Delta': (0, 4),
         'Theta': (4, 8),
         'Alpha': (8, 12),
         'Beta': (12, 30),
         'Gamma': (30, 45)}

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:
   freq_ix = np.where((fft_freq >= eeg_bands[band][0]) &
                   (fft_freq <= eeg_bands[band][1]))[0]
   eeg_band_fft[band] = np.mean(fft_vals[freq_ix])

print(eeg_band_fft)

I am getting the following errors when it goes into the for.
/home/xxxx/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3372: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
/home/xxxx/.local/lib/python3.8/site-packages/numpy/core/_methods.py:170: RuntimeWarning: invalid value encountered in double_scalars
ret = ret.dtype.type(ret / rcount)

What am I missing or what's wrong?
As I said above, I am very new to this matter so please, try to explain everything as simple as possible.
This is the result I get:
band        val
0  Delta  64.642945
1  Theta        NaN
2  Alpha        NaN
3   Beta        NaN
4  Gamma  25.851368

This is how the file I am reading looks like. It is a file  from a 16 electrodes reading with nearly 200k entries.
I have only taken the first entries of the file in order to test this code snippet.

There are a few things going on:
1.) Your data is actually only length 5 (made up of 5 longer components). This means when you call fft_freq = np.fft.rfftfreq(len(data), 1.0/fs) your array containing the frequencies of the FFT is only length 3! I think the main issue I'm having is that I'm not sure what your data actually represents: multiple electrodes?
2.) I think you're using np.where wrong. Basically it's not returning a slice which appears to be what you're trying to get.
3.) Even if a single array in your data represents an FFT of a single channel, the length is not nearly long enough to get anything useful. If we take the frequencies of data[0] we get [  0.   12.5  25.   37.5  50.   62.5  75.   87.5 100. ]. Not even one bin per band so you need to increase the length of your data. This is probably why freq_ix returns an empty array sometimes.
I was basically doing the same thing here in lines 45-68. A little different but more or less the same idea.
Also highly recommend using the MNE package. :)


help to calculate phase shift in certain frequency(identified from FFT) I am designing a 24ghz doppler radar, with 1TX 2RX. Four adc channel are recording at 25khz, applying 256point FFT on adc of both of RX I and Q, I could identify the target in certain doppler frequency(speed) . The next step is to calculate the phase shift to get the angle of the target, that's where I stuck. Since I locate the target on FFT, how to trace it back to the time domain and calculate the phase shift between the different RXs. Any input will be very welcomed. I am a computer engineer with limited knowledge of dsp though, but facinated about radar.

Phase shift in time domain translates to frequency shift in the Fourier transform.

$FT of exp(j2πf0t)x(t) = X(f−f0).$

If there is a phase shift in time, the Fourier transform that you take is already shifted in frequency. Find the shift in frequency f0 and you should get the phase change.


Trouble Generating a Sine Wave in C/Matlab from Basics this is a crosspost from stackoverflow as I figured there may be more who can help me here...

I've been writing a basic C program to generate a sine wave and write it out to STDOUT for piping into a filter which I am also going to write but I'm having issues.

I am able to generate what sounds like a good clean sine wave, but the frequency is off.  When I listen to the output I am generating and compare to a true signal generators output, the frequency of mine is slightly too high.

The C code should generate raw 16 bit signed stereo wav @ 44100Hz for a frequency passed as an argument to the program, or 500Hz as default.

A matlab version of the code is here, (modded very slightly for 1 indexing in Matlab, WHY MATHWORKS WHY!?) minus the passing to STDOUT etc as I know that works fine

CHANNELS = 2; SAMPLING_RATE = 44100; NUM_SAMPLES = 512;

frequency = 1000;

%% BEGIN WHILE(1) LOOP HERE

output_buff = zeros(1,CHANNELS*NUM_SAMPLES);

for i = 1:2:CHANNELS*NUM_SAMPLES
    output_buff(i) = 30000 * sin(frequency * pi * (i-1)/NUM_SAMPLES);
    output_buff(i+1) = output_buff(i);
end

%% OUTPUT TO STDOUT AND REPEAT

I should add that this code runs (in the C version) inside a while true loop, generating a full output_buff of values, then pushing the buffer to STDOUT.

I have written some further testing code to see what's actually being generated as follows:

plot(1:CHANNELS*NUM_SAMPLES, output_buff)

output_buff = output_buff .* hanning(length(output_buff))';
Y = fft(output_buff);
Mag=abs(Y(1:length(output_buff)/2)).^2;
[a,b]=max(Mag);
% Result
SAMPLING_RATE*b/length(output_buff)

When I run this script, I can see at the end the frequency of the generated signal is actually 1.0767e+03Hz... Close but no cigar...

I've tried tweaking some of the parameters but I have no idea what's wrong or how to make the generated frequency more accurate.

The C code is as follows:

#include <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <unistd.h>

#define M_PI 3.1415926535897932384626433832795L
#define NUM_CHANNELS 2
#define NUM_SAMPLES 1024
#define SAMPLE_RATE 44100

double frequency = 500.0;

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);

    while(1){
        for (int i = 0; i < cbBuffer; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 30000 * sin(frequency * M_PI * (double)i /(double)(NUM_SAMPLES));
            samples[i+1] = samples[i];
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}

UPDATE:

I changed the C code as follows and I now get a nice smooth sine wave but the frequency is way too low

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    int counter = 0;

    while(1){
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 3000 * sin(frequency * M_PI * (double)counter /(double)(SAMPLE_RATE));
            samples[i+1] = samples[i];
            counter++;
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}

So I've managed to fix it.  My problem was I wasn't properly taking into account the Sampling Rate of the produced raw wav file, and that by using the loop variable, there is a discontinuity for each iteration of the loop.

The final working code is as follows:

#include <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <unistd.h>

#define M_PI 3.1415926535897932384626433832795L
#define NUM_CHANNELS 2
#define NUM_SAMPLES 512
#define SAMPLE_RATE 44100

double frequency = 500.0;

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    int counter = 0;

    while(1){
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 3000 * sin(frequency * 2 * M_PI * (double)counter /(double)(SAMPLE_RATE));
            samples[i+1] = samples[i];
            counter++;
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}


Fixed point channel model for OFDM system I'm testing an OFDM system implemented in fixed point. The data format is Q11. My system work's fine but I need to test it under some channel for evaluation of the design before field testing. The system will be tested on Rician fading channel but for now I'm working on AWGN testing first. I asked a similar question here:
Apply AWGN noise to QPSK-OFDM symbol

But I realized that the problem is in the fixed point - float point conversion. I had a simulation of my system using 32b and the results agree with the theoretical BERXEbN0 for QPSK systems.

I'm using the following code for apply noise:

def awgn(self,data,noise_power):
    data_dq = data*(2.**(DATA_WIDTH-1))
    Es = sum(abs(data_dq)**2.)
    Eb = Es/(self.OFDM_SIZE*2.)
    EbN0_dB = noise_power + 10*log10(self.t_symbol/(self.t_symbol - self.t_cp))
    EbN0 = (10.)**(EbN0_dB/10.)
    N0 = Eb/EbN0
    noise = sqrt(N0/2.)*(random.standard_normal((data_dq.size,)) + random.standard_normal((data_dq.size,))*1j)
    return around(data_dq + noise)/(2.**(DATA_WIDTH-1))

where DATA_WIDTH = 12.

With this implementation the system agree for some regions of EbN0 curve and goes far from the curve in the other regions. With this number of bits I'm getting $10^{-2}$ BER with an EbN0 of 10dB.

My question is: how to handle the conversion between fixed and float point in this context?

update

I had changed some parts of the above code:

def awgn(self,data,noise_power):
    data_dq = data
    Es = sum(abs(data_dq)**2.)
    Eb = Es/(self.OFDM_SIZE*2.)
    EbN0_dB = noise_power + 10*log10(self.t_symbol/(self.t_symbol - self.t_cp))
    EbN0 = (10.)**(EbN0_dB/10.)
    N0 = Eb/EbN0
    noise = sqrt(N0/2.)*(random.standard_normal((data_dq.size,)) + random.standard_normal((data_dq.size,))*1j)
    return data_dq + around(noise*(2.**(DATA_WIDTH-1)))/(2.**(DATA_WIDTH-1))

I guess it will not affect the system performance but it's here for keep updated.

Below are the curve from my last simulation:

As you can see I have a floor in my present design.

I'm asking if this floor is due to my fixed point system performance or it's an effect of the channel model.

Update[Solution]

Just fro add an update to the problem for future reference. The problem was in the parameters of the fixed-point FFT core used. The effects in the figure above are due to clipping, as stated from the helpful comments in this topic. Thank you for your support.

The effect you're observing is most certainly due to quantization errors. If your calculation of noise power was wrong you would get an offset between the simulated and the theoretical curve but not an error floor. The fact that it works fine with 32 bit resolution supports that hypothesis. You can analyze the output signal of your VHDL transmitter in order to find out whether it's already corrupted. There are several things I would consider:

Calculate the BER without channel (back-to-back). It should be zero but from your diagram I guess it will be something around 0.02.
Plot the constellation diagram. You should see distinct points not "clouds".
Calculate the signal-to-quantization-noise ratio $\gamma_\mathrm{q}$. Let $a_{\nu k}$ be the complex value of subcarrier $\nu$ in OFDM symbol $k$ after mapping and $b_{\nu k}$ be the respective subcarrier value at the transmitter output (can be obtained by applying FFT to the tx output). Then
$$
\gamma_\mathrm{q}=\frac{\sum_{\nu.k}|a_{\nu k}|^2}{\sum_{\nu.k}|b_{\nu k} - a_{\nu k}|^2}
$$


Does a .wav file created from a .mp3 file have better quality than the .mp3 file itself? A .wav file created from a .mp3 file has the same frequency domain plot as the mp3 file itself https://stackoverflow.com/questions/58706674/difference-between-frequency-domain-in-mp3-and-wav.

However, the .wav file created from .mp3 file has larger file size. According the post https://stackoverflow.com/questions/15593806/why-is-a-wav-file-created-from-a-mp3-file-much-larger-in-size, this is because the .wav file is storing much more data than the mp3 file is.

Now, is there any difference in the quality between the .wav file and .mp3 file?

Does the quality of .wav file sound better since it has more data stored?

MP3 is a lossy compression format, so encoding an uncompressed lossless WAV file to MP3 will lose some quality in the encoding process.

Decoding a MP3 file back to WAV file has at most the same quality as the MP3 file has, it is just as good or bad as the contents in the MP3 file was, it is not in any way different to decode it for saving into a file or decoding it for playback. Most importantly, you don't get the lost quality back when decoding, and basically the result actually depends on the quality of the decoder as well.


Distinguishing FIR and IIR from difference equation Find the transfer function of the difference equation $$y_n = x_n + 1.2y_{n-1}$$

I fail to understand how one can distinguish between an FIR filter and an IIR filter by looking at the equation given above.

How to find if an IIR or a FIR filter is stable from the transfer function?

$$ H(z) = \frac{y(z)}{x(z)}$$

$$x_n = y_n-1.2y_{n-1}$$

$$x(z) = y(z)-1.2z^{-1}y(z)$$

$$H(z) = \frac{y(z)}{y(z) -1.2z^{-1}y(z)}$$

$$H(z) = \frac{z}{z-1.2}$$

Thus we arrive at the transfer function.

I am a week into signal processing. Any suggestions and tips appreciated. Please forgive my mistakes.

Welcome to DSP! First of all I would recommend reading chapter 2, specially 2.2 to 2.5 of Discrete Time Signal Processing 3rd edition by Oppenheim and Schafer. Then 3.2 and 5.2 of that very same book.

About your question. Well, in general (see comments on Axel answer for an exception) if you see an equation which depends on past output terms, like $1.2*y [n−1] $ in your equation then it is an IIR system. If it only depends on current, past or future inputs, it is an FIR system.

If it is an FIR system, you may determine stability, causality, linearity and time invariance from the Linear Constant Coefficient Difference Equation. However, if it is an IIR system, unless initial conditions are given to you, there's no way you could prove any of these because you don't know the value of the past output the first time you are going to "use" the system and it could determine whether it behaves one way or another.

For a system for which the input and output satisfy a linear constant coefficient difference equation, if the initial condition is that the system is initially at rest, then the system will be linear, time invariant, and causal. For stability analyze the Region of Convergence of the $Z$ transform and make sure it includes the unit circle


Which Fourier Transform properties should be applied first when they are combined together? A sequence $x[n]$ has Fourier Transform $X(e^{jw})$, and the Time shfiting, Frequency shifting and Time reversal properties are shown below, respectively :
$$x[n-n_d] \Rightarrow e^{-jwn_d}X(e^{jw}) $$
$$e^{jw_On}x[n] \Rightarrow X(e^{j(w-w_O)}) $$
$$x[-n] \Rightarrow X(e^{-jw}) $$

Suppose we were given that the Fourier transform of a signal $Y(e^{jw})$ is $3e^{j4w}X(e^{-{j(w-w_O)}})$, and we want to work backwards to find the original sequence.  If I use Time reversal first then:
$$x[-n] \Rightarrow X(e^{-jw})$$
followed by Frequency shifting:
$$e^{jw_On}x[-n] \Rightarrow X(e^{-j(w-w_0)}) $$
and then finally Time shfiting and Linearity:
$$3e^{jw_O(n+4)}x[-n+4] \Rightarrow 3e^{j4w}X(e^{-j(w-w_0)}) $$

So, I have obtained $y[n] = 3e^{jw_O(n+4)}x[-n+4]$,but if I had changed the order of properties which I applied,I would have gotten a slightly different answer. Should certain properties always be applied first when dealing with such problems?

If you do things right the order doesn't matter. Note that replacing $n$ by $n+4$ in $x[-n]$ results in $x[-(n+4)]=x[-n-4]$ and not in $x[-n+4]$.


Overlap-Add versus Overlap-Save What differences or other criteria can be used to help decide between using overlap-add and overlap-save for filtering?  Both overlap-add and overlap-save are described as algorithms for doing FFT based fast convolution of data streams with FIR filter kernels.  What are the latency, computational efficiency or caching locality (etc.) differences, if any?  Or are they the same?

Essentially, OS is slightly more efficient since it does not require the addition of the overlapping transients.  However, you may want to use OA if you need to reuse the FFTs with zero-padding rather than repeated samples.

Here is a quick overview from an article I wrote a while ago

  Fast convolution refers to the blockwise use of circular convolution
  to accomplish linear convolution. Fast convolution can be accomplished
  by OA or OS methods. OS is also known as “overlap- scrap” . In OA
  filtering, each signal data block contains only as many samples as
  allows circular convolution to be equivalent to linear convolution.
  The signal data block is zero-padded prior to the FFT to prevent the
  filter impulse response from “wrapping around” the end of the
  sequence. OA filtering adds the input-on transient from one block with
  the input-off transient from the previous block. In OS filtering,
  shown in Figure 1, no zero-padding is performed on the input data,
  thus the circular convolution is not equivalent to linear convolution.
  The portions that “wrap around” are useless and discarded. To
  compensate for this, the last part of the previous input block is used
  as the beginning of the next block. OS requires no addition of
  transients, making it faster than OA.


Which pattern (circle pattern or checkerboard pattern) should be used for automotive camera calibration (fisheye/wide webcam)? Calibrating a Fish-eye camera / webcam with large angle mounted on car by using opencv 3

One important issue is to define the calibration pattern and its position !!!!

The affects error that paly an essential role in the camera calibration process; of either localizing the control point or using geometric properties of the pattern directly.
Therefore, the precise localization of control points, or accurate determination of the geometric properties is very important, when a small error lead to imprecise camera calibration .Therefore I decided to use a regular grid of circles as calibration pattern.
But when using the wide angle lens the camera location has a big inclination with respect to the calibration pattern. Therefore, the size of the control points circles varies a lot (circle radius of the most remote circles could be smaller than those close to the lens).

According to that how can determine rows=? X  columns=?  of the calibration sheet, and the distance = ?  between each row an column, and the position=? for the calibration sheet, in order to obtain the accurate calibration for  the camera mounted on the car .

The answer here suggests that checkerboard patterns may yield more accurate (subpixel) calibration results and be more robust.

You may have edited your question because the title asks which pattern but the text asks about rows and columns. In either case, you may consider using checkerboard instead of circle pattern maybe?


How do you choose the optimal number of taps in a windowed image scaling filter? I'm attempting to use a filter to scale an image and for the most part, I have something that appears to work.  But now I want to be able to generate a filter with an optimal number of taps.  I know in general that the more taps, the better.  But there are diminishing returns in how much each additional tap will increase the quality of the image.  So how would one go about picking a number of taps that makes a good tradeoff between implementation complexity and image quality?  Are there standard measurements that would allow different filters to be evaluated to some "ideal" filter?

Thinking about this more, I think this is really two questions.  One is a matter of deciding which windowed sinc to use (in this case, I've already settled on using lanczos2 or 3) -- that is a more general question that has various tradeoffs depending on the application.  But once that decision is made, the number of taps is easily determined.  For lanczos2, the ideal number of taps is 4 and for lanczos3, the ideal number of taps is 6.  The reason it's not 5 or 7 is that one of those taps will always be zero due to the windowing.

In the more general case, I think the easy way to state how many taps to have for a given window is to simply say:

taps = max - min

In the case of lanczos2, the max is 2 and the min is -2.  Therefore, taps = 4.

Images quantized to some bit depth have error due to the quantization. Your filter also introduces error, compared to using an infinitely large ideal filter. Using larger and larger filters, eventually you will arrive at a knee in the total error vs filter size diagram where those two errors are equal. After that the quantization error dominates so it doesn't pay off much to improve the filter further.

A commonly used error metric is peak signal-to-noise ratio.


Filter to add 3dB per octave? If I want to add 3dB per octave to signal (e.g. to flatten out the power spectrum of an exponential sine sweep), is it really just as simple as...

Take the Fourier transform of the time-domain signal, let's call this A(f).
Multiply this by the square root of the frequency  (square root because power goes like amplitude squared)?
Maybe normalize by the lowest frequency, so that there's unity gain for that frequency.
Take the inverse Fourier transform to get back to the time domain.

Mathematically, the filter would be, simply,...

A'(f) := A(f) * sqrt ( f / f_low )

Is this right, and/or is there a better way?
(I've searched around and...it seems this is such a simple matter that people don't post how to do this.)

take any of the pinking filters you can find and swap the zeros with the poles.

e.g. here's my pinking filter with 3 real poles and 3 real zeros (all inside the unit circle):

$$ H(z) = A \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

$p_1$ = 0.99572754,     $q_1$ = 0.98443604

$p_2$ = 0.94790649,     $q_2$ = 0.83392334

$p_3$ = 0.53567505,     $q_3$ = 0.07568359

just swap the poles with the zeros.


kalman filter with time-varying noise? in regular discrete-time (1 dimensional) kalman filter, it is assumed that we have white gaussian noise affecting the transitions and the observations:

$x(t+1) = Ax + w$

$y(t) = Cx(t) + v$

assuming $w$ follows Gaussian(0, $\sigma_1$) and $v$ follows Gaussian(0, $\sigma_2$) distributions.

will the standard kalman filter not work or not be optimal if the gaussian noise variables are time varying, meaning $w$ and $v$ become $w(t)$ and $v(t)$ and somehow change with time, becoming systematically bigger or smaller?
related: the kalman filter generally has the behavior that its variance is decreasing with time. will that be true even if $w(t)$ and $v(t)$ are are time-varying even if they are assumed not to be? will the kalman filter's variance still get smaller even if the data are sampled from time-varying noise process?

1) It depends on what you call the standard Kalman filter -- I will call the equations in the picture below to be the "standard Kalman filter". You can easily derive an expression for the Kalman filter where the covariance matrices of the noise processes are time varying in terms of the covariances of the noise processes (and cross-covariance of the noise processes).

Moon and Stirling's Mathematical Methods and Algorithms for Signal Processing derives the Kalman filter under $w(t)$ has mean $0$, covariance $E[w(t) w(s)^T] = Q_t \delta_{ts}$, $v(t)$ is also mean zero with covariance $E[v(t) v(s)^T] = R_t \delta_{ts}$ and $E[w(t) v(s)^T] = M_t \delta_{ts}$ (if you assume $M(t) = 0$, the equations simplify considerably).

For the system model
$ x_{t+1} = A_t x_t + w_t , y_t = C_t x_t + v_t$  [$w_t = w(t), v_t=v(t)$] with $M_t = 0$, the following are the resultant Kalman equations (from Moon & Stirling) to get the predicted state estimate (13.45) and filtered state estimate (13.47):

These equations can be shown to be optimal using either a Bayesian approach (MAP with Gaussianity assumptions) or Linear MMSE approach for the given system model.

2) Convergence of the Kalman filter is a bit of an annoying subject. You can read these notes for more details or check the book Stochastic Systems: Estimation, Identification and Adaptive Control by Pravin Variya and P.R. Kumar or one of the other texts on stochastic systems. Essentially, if the systems are controllable and reachable and yada yada, you can get some statements on convergence properties.


Counter-example to show that this system is non-linear? I have a discrete-time system y[n] = x[n]*x[n-1]

I need to show that this system is non-linear by using a counter example rather than by disproving with algebra and the properties of additivity and homogeneity.

I think that the system is homogeneous although not additive, but I can't think of an example to show this. Would someone help me out please?

OK. Let's try homogeneity again:
$$
y[n] = x[n]\cdot x[n-1]
$$
for input $x[n]$.

For input $a\cdot x[n]$ we get
$$
y'[n] = a\cdot x[n] \cdot a x[n-1] = a^2 \cdot x[n]\cdot x[n-1] \not = ay[n]
$$
so homogeneity doesn't apply.

Let's try that with an example, as per MBaz's suggestion.

Let's try:
$$
x[n] = 2
$$
so that
$$
y[n] = 4
$$
Then let's try twice that:
$$
x'[n] = 2\cdot x[n] = 4
$$
So that
$$
y'[n] = 16 \not= 2\cdot y[n]
$$


Help designing Butterworth filter I am trying to bandpass filter a signal using a butterworth filter, but I am getting the following result (click to enlarge):

The original signal is the blue and the filtered one is the green.

Why does the first part of the filtered signal (from 0 to 50) is changing so much from the original? Can you help me understand why this is happening?

Here goes some details:

Signal length: 256
Frequency: 22Hz
Cutoff frequencies: 0.7Hz and 3.0Hz
Filter order: 6

Python code:

nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
y = lfilter(b, a, data)

The first part of the response may be a result of the transient response of the Butterworth filter.  Also note that the signal is phase shifted and the basic shape is altered.  This looks like the result of phase distortion.  You may want to try again using a Bessel filter which has better transient and phase characteristics.  If you are using DSP, you can also try a true linear phase FIR approach.  This will yield the best transient and phase characteristics.


Ideal Low Pass Filter Concept in MATLAB Please help me understand the following MATLAB code for Ideal Low pass filter. I am unable to understand the Part2 in the below code. Please explain me why we are doing like this.

I have read the Rafael C. Gonzalez's Digital Image Processing Using Matlab 2E which explains my question but I couldn't understand properly. It states that DFTUV Computes meshgrid frequency matrices. I don't understand this "meshgrid frequency matrices".It will be very much helpful if someone could explain me clearly.

 %IDEAL LOW-PASS FILTER

%Part 1
        function idealfilter(X,P) % X is the input image and P is the cut-off freq
        f=imread(X);  % reading an image X
        [M,N]=size(f); % Saving the the rows of X in M and columns in N
        F=fft2(double(f)); % Taking Fourier transform to the input image
%Part 2 % I don't understand this part
        u=0:(M-1);
        v=0:(N-1);
        idx=find(u>M/2);
        u(idx)=u(idx)-M;
        idy=find(v>N/2);
        v(idy)=v(idy)-N;
        [V,U]=meshgrid(v,u);
        D=sqrt(U.^2+V.^2);

%Part 3
        H=double(D<=P);       % Comparing with the cut-off frequency
        G=H.*F;               % Multiplying the Fourier transformed image with H
        g=real(ifft2(double(G))); % Inverse Fourier transform
        imshow(f),figure,imshow(g,[ ]); % Displaying input and output image
        end

I tried to run each commands in Part2 individually for M= 8 and N=8. I get

u=0:(M-1); ==> u = 0 1 2 3 4 5 6 7

v=0:(N-1); ==> v = 0 1 2 3 4 5 6 7

idx=find(u>M/2); ==> idx = 6 7 8

u(idx)=u(idx)-M; ==> 0 1 2 3 4 -3 -2 -1

idy=find(v>N/2); ==> idy = 6 7 8

v(idy)=v(idy)-N; ==> 0 1 2 3 4 -3 -2 -1

[V,U]=meshgrid(v,u); ==>

V=

     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1
     0     1     2     3     4    -3    -2    -1

U =

     0     0     0     0     0     0     0     0
     1     1     1     1     1     1     1     1
     2     2     2     2     2     2     2     2
     3     3     3     3     3     3     3     3
     4     4     4     4     4     4     4     4
    -3    -3    -3    -3    -3    -3    -3    -3
    -2    -2    -2    -2    -2    -2    -2    -2
    -1    -1    -1    -1    -1    -1    -1    -1

I am unsure why they are doing like this. How it affects the image.Please help me understand. The fftshift is used to avoid the wraparound error. Is the above operation anyway related to fftshift command? My understanding is that they are doing this to avoid using the fftshift command. Thanks in advance. Help me learn.

Reference: Rafael C. Gonzalez's Digital Image Processing Using Matlab 2E

meshgrid is simply a way to create a discrete grid. If you think of a regular MATLAB matrix of size MxN, it defines indices from 1 to N in x and from 1 to M in y, with the step size of 1. But let's say you want to evaluate some function f(x,y) on the x-interval [-3, 1] and the y-interval [-1, 2], and you want your step size to be .2.  The meshgrid helps you create a matrix of all the x-values and a matrix of all the y-values.  If you do

[X,Y] = meshgrid(-3:0.2:1, -1:0.2:2);

You can then easily evaluate your function f(x,y).  For example f(X(1),Y(1)) = f(-3, -1). Of if f is a vectorized MATLAB function, you should be able to evaluate it on your entire grid with one line:

Z = f(X(:), Y(:));

In this particular example, the meshgrids contain the indices of the frequency coefficients in the same order used by fft2 function. Remember, the Fourier transform of an image is periodic, it wraps around. For calculations it is convenient to have the zero frequency at (0,0), then the positive frequencies, and then the negative
frequencies wrapping around.

On the other hand, for visualization it is convenient to have the zero frequency in the center.  This is what the fftshift function does.  It is used purely for visualization.


Tetrahedral microphone array beamforming I have a microphone array of 4 channels taken as channels [6,10,22,26] from Eigenmike spherical microphone array. I wish to do 3D beamforming and create $8\times 4$ beams, 4-elevations with 8 wavefronts per elevation equally distributed.

I can't seem to work out the math as to the delays between the microphones. I also can't find an "off the shelf" code for that. This is my first beamformer so I am not even sure I am not missing anything. Is delay and sum process all I need?

What are the equations for determining the delays between the channels per beam?
Do I need any other equations to create a beamformer?

Specifically, I am trying to create beams in direction from all possible combinations of $\varphi_s \in \{-\frac{3\pi}{8},-\frac{\pi}{8},\frac{\pi}{8},\frac{3\pi}{8}\}$ and $\theta_s\in\{-\frac{7\pi}{8},-\frac{5\pi}{8},-\frac{3\pi}{8},-\frac{\pi}{8},\frac{\pi}{8},\frac{3\pi}{8},\frac{5\pi}{8},\frac{7\pi}{8}\}$.

My microphone array contains 4 microphones $m_i=(r_i,\theta_i,\varphi_i)$, at the following positions:
 \begin{pmatrix}
  0.042 & \frac{\pi}{4}  & \frac{7\pi}{36} \\
  0.042 & -\frac{\pi}{4} & -\frac{7\pi}{36} \\
  0.042 & \frac{3\pi}{4} & -\frac{7\pi}{36} \\
  0.042 & -\frac{3\pi}{4} & \frac{7\pi}{36}
 \end{pmatrix}

$r_i$ is in meters.
$\varphi_i, \theta_i$ are azimuth and elevation respectively, according to the microphone's documentation, given herein radians.

A quick sketch of the configuration from two views. Note that the distance of all mics is 42mm from the origin:

$(x,y,z)$ values in this system is extracted in the following manner to my understanding:
$$x=r\cdot \cos(\varphi)\cdot \cos(\theta)$$
$$y=r\cdot \cos(\varphi)\cdot \sin(\theta)$$
$$z=r\cdot \sin(\varphi)$$

I don't know what the fuss is about.  This is a simple vector problem.

"Well, if it is that easy please elaborate. I have 4 microphones on the following locations (0.042,45,35),(0.042,-45,-35),(0.042,135,-35),(0.042,-1355,35). These are coordinates for the microphones with m1=(ri,θi,φi), where −π≤θ≤π and −π2≤φ≤π2. create a beam towards a sourse at (θs,φs)=(22.5,67.5). – havakok"

You need to define five vectors in Cartesian coordinates, one for each mic and one for the beam.  The beam vector should be normalized to a length of one.

If the sound source is sufficiently far away, the assumption that the wavefront is coming in as a planar wave is a fairly accurate approximation.  Therefore, you simply have to find the projection of each mic location onto the line defined by the beam vector.  The projection for each mic is found by dotting the mic vector with the normalized beam vector.  The value you get is the distance along the line (where one is the beam vector length).  Use the speed of sound to convert this distance to a time shift.  Use your sampling rate to convert the time shift to a sample shift.  Shift and add.  The last mic to receive the signal is not delayed, and the previous mics are delayed relative to their projection's distance to the last mic's projection.

It should be obvious that you only get one beam direction at a time.  Don't be oversold on the "beam" notion either.

Ced

Two Answers:

1) Rounding to the nearest sample should be negligible.

Some ball park figures to consider:  Suppose you have a 1000Hz tone.  Sound travels at roughly 1130 feet per second.  Thus one cycle of the tone is roughly 1.130 feet long, or about 13 and 1/2 inches.  At CD quality sampling, about 44 samples.  So each sample is about 3/10 of an inch (this is frequency independent).  Compared to the length of the cycle, this isn't a lot.  However, if you are concerned about frequencies above 10k, this may matter.  In that case, a simple linear interpolation in your delay should do the trick.

2) The beam effect is formed by "coordinated interference effects", i.e. signals coming from other directions will be added together in an out of phase manner (only the signals in the desired direction are in phase), thus you are relying on interference to destructively interfere.  This is frequency dependent.  Also, signals coming in nearly at the beam angle won't be affected much, so perhaps "cone forming" is more accurate.

You didn't mention whether these were omnidirectional mics or not.  No bearing on the math, but may influence your results.