Difference between acoustic emissions and ultrasonic waves I need some clarification about the difference between acoustic emissions and ultrasonic waves. I have the idea that both are high-frequency signals and informative frequency bands more or less overlap.  Is there any difference between the waveform or in other characteristics that may the propagation speed through a medium. Moreover, will the velocity of both the signals through a medium (carbon steel) will be the same or will there be any difference.

The propagation of a mechanical wave within a medium depends on the characteristics of that medium (e.g. density).

An easy way to experience this is to knock a piece of material with a sharp hit and listen to it trying to dissipate the energy you gave it. Consider for example a bell. When you hit it, it goes bonnnnnnnnnnng at some tone. If you were to weld a few rings around the bell and try to hit it again, it would go bonng. Different density / stifness of the material, different sound.

Acoustic emissions are akin to the sound of the bell and the cause of the "ringing" are the internal cracks. Consider for example the sound of ice cracking (similar sounds are produced even as ice cracks on its own, without forcing it externally) or even closer to your question, glass cracking under thermal shock.

Within a block of ice, the acoustic emissions are low enough in frequency that can be heard. But, if that was a piece of kevlar breaking, the disturbance would propagate with different characteristics through the material. You can see that we are already at hundreds of kHz in this case and we would be needing special equipment to sense the material vibrating.

Hope this helps.


Initial rest condition applied on $x(t)$ vs $h(t)$ Define the LTI system $\mathcal{H} : x\mapsto y$
Define the convolution for continuous-time system :
$$
(x*h)(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)\;\text{d}\tau
$$
The initial rest condition states that :

An LTI system with an input signal $x(t)$ is causal if and only if $x(t)=0$, $\forall t<0$

Now I have noticed that in some textbooks, the definition of convolution would have $h(\tau)$ instead of $x(\tau)$ and $x(t-\tau)$ instead of $h(t-\tau)$ (since $x*h=h*x$) but the initial rest condition now becomes :

An LTI system with an impulse response $h(t)$ is causal if and only if $h(t)=0$, $\forall t<0$

My question is :

Why does commutativity allow that for a LTI system to be causal then in one definition $x(t)=0$, $\forall t<0$ while in another definition $h(t)=0$, $\forall t<0$ knowing that $x$ and $h$ are two different terms?

It should be clear that a property of a system, such as causality, cannot be determined by looking at its input signals. For a linear time-invariant system, it is its impulse response $h(t)$ from which properties such as causality or stability can be determined.
Only the second definition in the question is correct: a causal LTI system has an impulse response $h(t)$ which equals zero for $t<0$. This implies that the system's output at time $t_0$ only depends on values of the input signal $x(t)$ for $t\le t_0$, and not on future values $x(t)$, $t>t_0$.
The output signal of a causal LTI system can be written as
$$y(t)=\int_{0}^{\infty}h(\tau)x(t-\tau)d\tau=\int_{-\infty}^tx(\tau)h(t-\tau)d\tau$$
An initial-rest condition just means that if the input $x(t)$ is zero for $t<t_0$ then the output $y(t)$ must also be zero for $t<t_0$.


Tuning a PID Controller I have the following transfer function
$$10000 \cdot\frac{(s-0.012)(s+1.05)(s-18)}{(s+0.22)^2(s+45)(s+1000)}$$
For which I am trying to tune a PID controller for. I'm using the pidTuner in MATLAB to do so. I wish to have the settle time be less than 3 seconds if possible, however, the best I have been able to do is 77 seconds with the settings $K_p = 0.049955$, $K_i = 0.0087525$ and $K_d=0$. At these values, going in any direction seems to increase the settling time more. Any help would be really great!

You have 2 right-half plane zeroes : 0.012 and 18. The zero at +18 is "fast" and will not affect the performance much, but your slow zero at 0.012 will severely limit your performance. You can't cancel this right-half plane zero with a right-half plane pole in your controller, your controller output will be unbounded.
Is this homework or a real life problem?


DOA estimation; Music (Multiple Signal Classification) algorithm Why do we say that signal subspace and noise subspace are orthogonal in e.g. MUSIC algorithm? precisely, Suppose we have 2 sources (signals) which are WSS, in WGN in a uniform linear array. After computing the covariance matrix of received data "Y=AS+N", and decomposing R (covariance matrix) into 2 subspaces, we say that these 2 subspaces are orthogonal. Why?

The receive-signals covariance matrix $\mathbf R$ is by its definition
$$\DeclareMathOperator{\expect}{E}\newcommand{\E}[1]{\expect\left[#1\right]}\DeclareMathOperator{\diag}{diag}$$
$$\mathbf R= \E{\mathbf{YY}^*}=\mathbf R^*$$ hermitian.

We don't just find any decomposition of $\mathbf R$ into subspaces – we pick the Eigenvalue Decomposition, which gives us

$$\mathbf R= \mathbf U \mathbf \Lambda \mathbf U^* \text,$$

where

\begin{align}
\mathbf\Lambda &= \diag(\lambda_1,\ldots,\lambda_N) &\text{diagonal matrix of Eigenvalues,}\\
\mathbf U &= ( \mathbf g_1,\ldots, \mathbf g_N)&\text{matrix: Eigenvectors as columns}
\end{align}

We can without loss of generality have $\mathbf \Lambda$ sorted such that $\lambda_1,\ldots\lambda_M$ and thus $\mathbf g_1,\ldots \mathbf g_M$ correspond to transmit signals, and $\lambda_{M+1},\ldots,\lambda_N$ and $\mathbf g_{M+1},\ldots,\mathbf g_N$ come from to noise.

Assuming our receive array is well-formed and the signal is not incident from an angle that leads to ambiguities, $\lambda_i\ne\lambda_{l\ne i}$.
Now, in the presence of uncorrelated noise, we know that

$\mathbf R$ has full rank, and
since $\mathbf R$ is hermitian, different Eigenvalues have orthogonal Eigenvectors.

So, if we have $M<N$ independent receive signal components, we get $M$ different, orthogonal vectors.

The $M+1,\ldots, N$th eigenvalues of $\mathbf R$ are the noise variance $\sigma^2\ne\lambda_i,\,\forall i\in\left\{1,\ldots,M\right\}$, we see that the space spanned by the signal-caused Eigenvector must be orthogonal to the noise subspace.


Optimal image histogram for Human eye I need to visualize some data as gray-scale images after contrast enhancement. The processing algorithm must be built on few samples, and then it will be applied on unseen images (suppose, the [min, max] range is given).

I plan to:

Use sample images to estimate actual histogram H1.
Define desired histogram H2.
Build a transformation T (look-up table) that transforms H1 into H2.
Apply T on unseen images.

Question: How to build H2? Are there studies in regards to the optimal gray-scale histogram for Human eye?

Below are two examples of different transformations, applied on the same image:

Human eye perception is a complicated thing, the human perception tend to fail in different situations even in large changes of gray tones, also is capable of note minimal changes in the gray levels but all depends on the shape of the figure rather than its histogram.

In some cases you can enhance an image by applying histogram equalization techniques but this is only to emphasize some gray tones in order to observe more details in an image. This is related to gray images.


Sobel and Prewitt operators decomposition I was reading in the wikipedia page of Sobel operator and Prewitt operator that is possible to decompose these two operators (quote form the Formulation paragraph):

  "as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing."

I know that this means that I can rewrite and simplify a Prewitt mask (for example) in this way:

\begin{align*}
\begin{bmatrix}
+1 & 0 & -1 \\
+1 & 0 & -1 \\
+1 & 0 & -1
\end{bmatrix} = \begin{bmatrix}
1\\
1\\
1
\end{bmatrix}
\begin{bmatrix}
+1 & 0 & -1
\end{bmatrix}
\end{align*}

But I don't understand what are in these cases the averaging and the differentiation kernel? And, why it is written that they can compute the gradient with smoothing?

The part:
\begin{bmatrix}
1\\
1\\
1
\end{bmatrix}
acts on vertical pixels
\begin{bmatrix}
x_{1,\cdot}\\
x_{2,\cdot}\\
x_{3,\cdot}
\end{bmatrix}
as a weighted sum with equal weights, and if divided by $3$ like a 3-point average $(x_{1,\cdot}+x_{2,\cdot}+x_{3,\cdot})/3$.

Then,
\begin{bmatrix}
+1 & 0 & -1
\end{bmatrix}
acts on a row of horizontal pixels
\begin{bmatrix}
x_{\cdot,1}&
x_{\cdot,2}&
x_{\cdot,3}
\end{bmatrix}
like the discrete 3-point derivative: $(x_{\cdot,1}-x_{\cdot,3})/2$ up to the $2$ factor. Here, up to a global $6$ factor, it acts like a smoother vertically combined with a gradient horizontally. The same reasoning works when you switch directions. This notion relates to the separability of a 2D mask into the product of two 1D operators.

Additional sources:

What's logic behind the constuction of Sobel's filter in image processing?
How to find out if a transform matrix is separable?


Fourier transform of frequency division I need a mathematical description of what happens to the spectrum of a signal when that signal is processed to divide its frequency.

As an example, construct an input signal as the sum of two sine waves of different frequencies (as commonly used in a "two tone" test). Processing this signal with a buffer that changes state at the rising zero crossings of the input signal produces a squared wave whose frequency can be divided by a flip-flop:

The spectrum of the the input signal comprises tones at 102kHz and 105kHz:

"Squaring" the input preserves the two equal tones, albeit with added side tones (spaced at the difference between the two original tones) that result from removal of the amplitude envelope and added harmonics from the squaring:

Passing the squared signal through a flip-flop produces a second square wave which in which the two-tone nature is no longer present:

The single dominant tone is at 1/2 the frequency of the higher of the original two tones, while the spacing between the peaks continues to be the spacing between the two original tones.

Passing the squared result through an LC low-pass filter produces a sine wave with cycles that are half as frequent as the input:

with the same general spectral distribution as the squared output:

Can someone, please, explain what happened to the other tone?

A sum of two sine waves can be expressed equivalently as the product of two sine waves.  This is a phenomenon commonly referred to in musical circles as a beat frequency, which can be used to tune instruments by ear.  For your example, your input can be expressed as a 103.5kHz frequency multiplied by a 1.5Hz frequency.

I wouldn’t use the term “squaring” for what you are doing as much as “clipping the ever living hell out of”.  This would probably be expressed best using concepts from intermodulation distortion.  However, in your application, the square nature of the output of the waveform can be analyzed more simply anecdotally.  I would expect to see a square wave at the 103.5kHz with glitches at the 1.5Hz zero crossings.  But notably, it would look a lot like a 103.5kHz waveform, with some side band nonsense caused by glitches.

I do not believe that the flip-flop can be expressed using any common signal processing techniques, but again it can be analyzed anecdotally.  As you’ve hypothesized, I would expect the time between edges to double, because you’re skipping every other edge.  However, this is a grosse generalization.  It would be exactly true for a sinusoid, less true for the input you’ve proposed, and untrue for a generic input.  So for you’re example, mostly energy at 51.75kHz.

At this point, you still have a square wave, roughly at 51.75kHz, with some side band garbage cause by the glitches.  One could apply an arbititary amount of filtering to extract the sine wave at this frequency, resulting in the output waveform you described.

Something else to keep in mind.  A frequency transform has practical limitations.  A proper transform requires all time domain points of the input signal.  Using less than infinite samples is equivalent to windowing, which causes side lobe patterns about dominant frequencies, which is a topic unto itself.  Additionally, a signal sampled at discrete time intervals (like on a computer) has a limited bandwidth.  When you “clip the ever loving hell” out of a signal, you are widening the bandwidth of the input signal significantly.  If care isn’t taken to ensure that the discrete time sample rate is high enough, a significant amount of aliasing distortion will occur, hindering proper analysis.

Sorry for being a bit of a rant, but your question touches on a lot of points.  Cheers.


why does quantisation noise sound like white noise? To record sound and transform it to digital form, audio signals are recorded and its pressure levels are sampled and quantised.

Why does quantisation noise sound like white noise?

It often doesn't sound like white noise. The quantization noise of a sine wave is fairly harmonic with a clearly detectable pitch. The quantization of a pulse train sounds like a pulse train.

If the quantization has enough bits, and if the signal is reasonably stationary in time and reasonably dense in the frequency domain, the quantization noise will be more or less white. For music quantized at 16 bits, that's often the case. However the quantization noise will always be correlated with the original signal which can produce undesirable artifacts. Hence commercial audio production uses dithering and noise shaping and decrease the audibility of quantization.


confused in difference between calculations of convolution and multiplication? How are two signals multiplied? And how is it different from convolving two signals??

I have read the above question on dsp stack exchange and i also tried it in my matlab , code is as follow:

clc
clear all
close all
x=[2,4,1]
p=[5,1,8]
y1=x.*p   % result of multiplicaton
y2=conv(x,p) % result of convolution

But Matlab gives answer as given in that link but i am confused,about convolution steps, why only 5(first element of 2nd vector ) is multiplied with all there 3 element of first vector making up 3 different terms of final answer(as shown in attached photo) but 2nd element(blue under line) & 3rd element(red under line) of 2nd matrix do not repeat the same steps as followed for first element?why they are not multiplied by all elements of first vector as was the case with 5

?

Just write down the convolution sum to see what's going on:

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]p[n-k]\tag{1}$$

where we define the elements of the sequences $x[n]$ and $p[n]$ as equal to zero for the index $n$ outside the range of non-zero values, i.e., outside $n\in[0,2]$. Taking into account those zero values, we can rewrite $(1)$ with finite summation limits:

$$y[n]=\sum_{k=0}^{n}x[k]p[n-k]\tag{2}$$

because for $k<0$ all $x[k]$ are zero, and for $k>n$ the elements $p[n-k]$ are zero. Now you see that as $n$ increases from $0$ to  $2$ (outside that range the term $x[k]$ will be zero), you'll get more terms in the sum $(2)$:

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


duty cycle of a WiFi Packet I am doing a simple over-the-air(WiFi 802.11a) test in a shielded chamber involving one Tx and one Rx (around 5 feet).. I need to find the duty cycle given the length of the payload (in bytes), MCS etc.
As far as I understand , duty cycle is time a device spent on channel / total time .
I was going through this paper https://www.academia.edu/15142561/Determination_of_the_duty_cycle_of_WLAN_for_realistic_radio_frequency_electromagnetic_field_exposure_assessment
Where the paper calculates the theoretical upper limit of duty cycle for a given MCS

The payload is 1500 bytes
Max data rate is 54 Mbps

How is the data 248 micro secs in calculations? I am thinking,it should be 222.22 micro seconds (1500*8/54 micro seconds)

can any one please let me know  about it

And If anyone has done such analysis on duty cycles , please do guide me thanks

Wifi uses OFDM, which uses symbols that can transport many bytes at once. That means you can't have arbitrary long packets, but always need to use the next multiple of a symbol duration.


Given a continuous time signal, does the minimum Nyquist sampling rate depend on the choice of the set of basis functions? This is my first question on a StackExchange.  When the basis functions to represent a signal are chosen as  $e^{j\omega t}$ such as in a continuous-time Fourier transform then the sample rate $f_\text{s}$ must be more than twice the maximum frequency taken by the chosen continuous time signal. Two questions:

Is it possible to get a Nyquist sampling rate less (I doubt) or more than this by choosing a different basis set? I understand that "frequency" itself may no longer be well defined, so basically: would it be possible to choose a set where one could afford to sample at a higher/lower rate without loss.
If yes, is there some reasoning behind how the set was chosen? Any theory/references would be appreciated.

Would be glad to clarify any doubts/inconsistencies.

In the most general case, if you want to sample a continuous-time signal without loss of information, the minimum sampling rate is independent of any choice basis functions. The faster the signal changes with respect to the independent variable (which doesn't need to be time), the faster you have to sample. And if the signal is not band-limited or if it can't be parameterized in some other way, sampling will be lossy.

Note, however, that the (Nyquist-Shannon) sampling theorem is only a sufficient criterion for sampling a continuous-time signal without loss of information. In general you have to choose the sampling frequency greater than twice the highest frequency component in the signal. However, under certain conditions you can get away with sampling at a lower rate.

One example are bandpass signals, which can often be sampled at rates much lower than their highest frequency, especially if their bandwidth is small compared to their center frequency. But note that it is generally not sufficient to sample at a rate equal to twice the signal's bandwidth. See this answer for more details on bandpass sampling.

Another special case are sparse signals. The sparsity can be exploited to represent these signals with fewer samples than required by the sampling theorem. The corresponding technique is called compressed sensing.

In general, if a given signal has a "finite rate of innovation (FRI)" (such signals are called "FRI-signals"), they can be sampled at or above their rate of innovation. This subject treated in this important paper.

I also recommend reading this overview paper: Sampling - 50 Years After Shannon.


Is Step the same as an Impulse + 1/s? Isnt a Step (as a Transfer function -> $ \frac{1}{s} $) and Impulse (trans fcn -> 1) therefore Impulse + $ \frac{1}{s} $ should be equal to a Step in Simulink right?
I made a Simulink model:

But the outputted plots are not the same ?
Why is that? Am I wrong? Should those two blocks be essentially the same?

The problem is that you can't use a Dirac delta impulse as an input, which would be the theoretically correct thing to do. Since you only have a discrete impulse, which has some finite value, and since Simulink probably uses something like a zero-order hold to make a continuous-time signal out of that impulse you actually use two different input signals: a step and a smoothed step.


Spatial domain vs. Frequency domain filtering of an image. Which one is better? The question might seem naive but I'm new to this and a little confused.
So, I recently came to know that gaussian smoothing can be applied to an image both in spatial domain and in frequency domain. But I am confused on how should the results be compared?
So, suppose I apply gaussian to an image A in spatial domain and get output GS(A), and similarly get GF(A) when applied in frequency domain. Which one of both these will be better? I know this is dependent on filter size and all but I am confused on how can the comparison be made.
Also, is there a standard practice where we use frequency domain filtering for certain tasks while spatial domain filtering for certain other tasks? If yes, can you please share your views and/or reading material for the same?
EDIT 1 (CLARIFICATION):
I know that "you can implement spatial filtering either by direct convolution or indirectly by multiplying the image's Fourier transform by the filter's Fourier transform, then taking the inverse Fourier transform to get the whole thing back into the spatial domain".
What I want to know is the difference in output of these two methods.
Also, which method should I opt for and why?
There might be speed differences obviously, but are there any quality changes?
When should I use which method? (if it depends on the use case)

OK.  In theory, convolution in the time domain is the same as multiplication in the frequency domain.  And filtering a signal with a linear filter is mathematically the same as convolving by its impulse response.
So in the time domain, if you have an image $x_{n, m}$, and a filter with an impulse response $h_{n, m}$, then after filtering $x$ with $h$, you get a signal $y = x * h$.  (Note that especially in image processing, a filter's impulse response is often called a kernel.  Until you're smarter about this than I usually am, just treat "impulse response" and "kernel" as synonyms).
If you let $X = \mathcal F \left \lbrace x \right \rbrace$ denote the Fourier transform, then if $y = x * h$, $Y = H X$.  This turns out to be a happy result from 1-dimensional Fourier analysis, and an even happier result is that if you're using Cartesian coordinates, it applies to 2-dimensional (on 3- or 4-, or 42-dimensional) Fourier analysis, too.
That's the theory.  Basically, it's what you've found in your research*.
So if you have $h$ then you can compute $y$ as $y = x * h$ (or $y = h * x$ -- convolution commutes).  But you can also compute $y = \mathcal F^{-1} \left \lbrace \mathcal F \left \lbrace \mathcal h \right \rbrace F \left \lbrace \mathcal x \right \rbrace \right \rbrace $.  This seems to be dreadfully counter-productive, but for two things.
First, you can pre-compute $H = F \left \lbrace \mathcal h \right \rbrace$.  If it's settled enough, you can even compute it during design time and just have a pre-defined array in code.
Second, there's this thing called the fast Fourier transform (FFT).  From a mathematical perspective, a fast Fourier transform would be just an overly complicated way of doing a plain old discrete Fourier transform except that for large images it uses far fewer operations (hence -- Fast Fourier Transform).
Convolution of an $N \times M$ image by a $N_h \times M_h$ kernel takes -- if I'm getting my math right -- $N M N_h M_h$ multiply-add operations.  Basically, it's the area of the kernel times the area of the image.  Taking the FFT and inverse FFT of that same image area takes $k N M \log(N M)$ multiply-add operations**.  So for a large enough filter kernel, using the FFT can be faster -- particularly because you can use a technique called "overlap-add" to speed things up.
So, theoretically there is no difference between doing your filtering using convolution in the spatial-domain or by multiplying in the frequency domain.
Practically, doing it in the frequency domain is always more complicated, but it can be faster.  Also practically, when you use the FFT there can be some slight loss of numerical precision -- basically, the speed comes at the cost of a bit of extra numerical scrud due to rounding errors and whatnot.  This numerical effect is almost certainly of no concern to you.
Finally -- practically --

it's rare that you have useful filter kernels that are large enough that doing the computation in the frequency domain will speed things up.
Moreover, if all you're interested in is Gaussian filtering, there are some fast almost-Gaussian filtering technique that apply a small (3x3, I think) kernel to an image repeatedly.  Because of the Central Limit Theorem, if you choose a small kernel correctly and filter by it repeatedly, the resulting filtering will tend to approximate filtering by a Gaussian.

* Things may have been a bit clouded by the fact that there's a result that's specific to Gaussian filters, in that if $h_{n, m}$ is Gaussian, then so is $H(\omega_1, \omega_2)$ -- it's just that in general it has a different scaling.  This is not always the case, but if you're looking at some page that doesn't fully explain it, the notion of multiplying by a Gaussian in the frequency domain being the same as convolving by a Gaussian in the spatial domain would look extra-weird.
** I don't remember the exact number, but I think for a naïve FFT implementation, $k=6$, and there's more clever implementations that get it smaller yet.


Which transform most closely mimics the human auditory system? The Fourier transform is commonly used for frequency analysis of sounds. However, it has some disadvantages when it comes to analyzing the human perception of sound. For example, its frequency bins are linear, whereas the human ear responds to frequency logarithmically, not linearly.

Wavelet transforms can modify the resolution for different frequency ranges, unlike the Fourier transform. The wavelet transform’s properties allow large temporal supports for lower frequencies while maintaining short temporal widths for higher frequencies.

The Morlet wavelet is closely related to human perception of hearing. It can be applied to music transcription and produces very accurate results that are not possible using Fourier transform techniques. It is capable of capturing short bursts of repeating and alternating music notes with a clear start and end time for each note.

The constant-Q transform (closely related to the Morlet wavelet transform) is also well suited to musical data. As the output of the transform is effectively amplitude/phase against log frequency, fewer spectral bins are required to cover a given range effectively, and this proves useful when frequencies span several octaves.

The transform exhibits a reduction in frequency resolution with higher frequency bins, which is desirable for auditory applications. It mirrors the human auditory system, whereby at lower-frequencies spectral resolution is better, whereas temporal resolution improves at higher frequencies.

My question is this: Are there other transforms which closely mimic the human auditory system? Has anyone attempted to design a transform that anatomically/neurologically matches the human auditory system as closely as possible?

For example, it is known that human ears have a logarithmic response to sound intensity. It is also known that equal-loudness contours vary not only with intensity, but with the spacing in frequency of spectral components. Sounds containing spectral components in many critical bands are perceived as louder even if the total sound pressure remains constant.

Finally, the human ear has a frequency-dependent limited temporal resolution. Perhaps this could be taken into account as well.

In designing such transformations, one should take into account competing interests:

fidelity to the human auditory system (that varies with people), including non-linear or even chaotic aspects (tinnitus)
easiness of the mathematical formulation for the analysis part
possibility to discretize it or allow fast implementations
existence of a  suitable stable inverse

Two recents designs have catch my ears recently:
Auditory-motivated Gammatone wavelet transform, Signal Processing, 2014

  The ability of the continuous wavelet transform (CWT) to provide good
  time and frequency localization has made it a popular tool in
  time–frequency analysis of signals. Wavelets exhibit constant-Q
  property, which is also possessed by the basilar membrane filters in
  the peripheral auditory system. The basilar membrane filters or
  auditory filters are often modeled by a Gammatone function, which
  provides a good approximation to experimentally determined responses.
  The filterbank derived from these filters is referred to as a
  Gammatone filterbank. In general, wavelet analysis can be likened to a
  filterbank analysis and hence the interesting link between standard
  wavelet analysis and Gammatone filterbank. However, the Gammatone
  function does not exactly qualify as a wavelet because its time
  average is not zero. We show how bona fide wavelets can be constructed
  out of Gammatone functions. We analyze properties such as
  admissibility, time-bandwidth product, vanishing moments, which are
  particularly relevant in the context of wavelets. We also show how the
  proposed auditory wavelets are produced as the impulse response of a
  linear, shift-invariant system governed by a linear differential
  equation with constant coefficients. We propose analog circuit
  implementations of the proposed CWT. We also show how the
  Gammatone-derived wavelets can be used for singularity detection and
  time–frequency analysis of transient signals.

The ERBlet transform: An auditory-based time-frequency representation with perfect reconstruction, ICASSP 2013

  This paper describes a method for obtaining a perceptually motivated
  and perfectly invertible time-frequency representation of a sound
  signal. Based on frame theory and the recent non-stationary Gabor
  transform, a linear representation with resolution evolving across
  frequency is formulated and implemented as a non-uniform filterbank.
  To match the human auditory time-frequency resolution, the transform
  uses Gaussian windows equidistantly spaced on the psychoacoustic “ERB”
  frequency scale. Additionally, the transform features adaptable
  resolution and redundancy. Simulations showed that perfect
  reconstruction can be achieved using fast iterative methods and
  preconditioning even using one filter per ERB and a very low
  redundancy (1.08). Comparison with a linear gammatone filterbank
  showed that the ERBlet approximates well the auditory time-frequency
  resolution.

And I shall mention also:

An Auditory-Based Transform For Audio Signal Processing, WASPAA 2009

  An auditory-based transform is presented in this paper. Through an
  analysis process, the transform coverts time-domain signals into a set
  of filter bank output. The frequency responses and distributions of
  the filter bank are similar to those in the basilar membrane of the
  cochlea. Signal processing can be conducted in the decomposed signal
  domain. Through a synthesis process, the decomposed signals can be
  synthesized back to the original signal through a simple computation.
  Also, fast algorithms for discrete-time signals are presented for both
  the forward and inverse transforms. The transform has been approved in
  theory and validated in experiments. An example on noise reduction
  application is presented. The proposed transform is robust to
  background and computational noises and is free from pitch harmonics.
  The derived fast algorithm can also be used to compute continuous
  wavelet transform


Similarity between two Signals: cycle by cycle I need to find the similarity between two signals, i.e., measured signal and the reference signal(e.g. a simple sinusoidal signal). What I mean by similarity is, if there is any decentralized glitch or any change in freq for one particular cycle, the similarity measure should give 'False' or a value in 0 - 1 scale for that cycle. Phase difference also should be considered. Could this be done using correlation? Or is it possible to use wavelets? Or any other method?

Many Thanks...

so i'll throw the cross-correlation and/or matched-filter thing at you.  a little bit similar to Karthik but maybe a little more rigor.

so you have a reference signal, $y(t)$, and a measured signal $x(t)$.  suppose they have something in common and some differences and may one of the differences is a misalignment in time.  but if you investigate what these two signals have in common (or what they don't) for a variety of different offsets in time, that dissimilarity due to time misalignment can be made non-salient (or maybe you want to know the time offset).

so first consider measuring the dissimilarity (at a variety of different offsets or "lags") as:

$$ Q_{xy}(\tau) \triangleq \int\limits_{-\infty}^{+\infty} \left( x(t) - y(t+\tau) \right)^2 dt $$

so you're looking at how much $x(t)$ and $y(t+\tau)$ are different, for different lags, $\tau$, squaring that difference or "error" so that it's always positive (this is analytically easier than if you use the absolute value), adding up the total squared-difference over all of time and you have an overall measure in how dissimilar $x(t)$ and $y(t+\tau)$ are.  then picking the lag $\tau$ so that $Q_{xy}(\tau) \ge 0$ is minimum might be considered to be a measure on how dissimilar $x(t)$ and $y(t+\tau)$ are.

the smaller this measure (which is always non-negative), the more similar $x(t)$ and $y(t+\tau)$ are.  this is the Average Squared Difference Function (ASDF) and has a similar motivation in its definition as does the more well-known and historical AMDF function.

now the cool thing you get that you don't if it's AMDF is

$$ \begin{align}
Q_{xy}(\tau) & = \int\limits_{-\infty}^{+\infty} \left( x(t) - y(t+\tau) \right)^2 \ dt  \\
& = \int\limits_{-\infty}^{+\infty} x^2(t) + y^2(t+\tau) - 2x(t)y(t+\tau) \ dt \\
& = \ \int\limits_{-\infty}^{+\infty} x^2(t) dt \ + \ \int\limits_{-\infty}^{+\infty}y^2(t+\tau) dt \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
& = \ \int\limits_{-\infty}^{+\infty} x^2(t) dt \ + \ \int\limits_{-\infty}^{+\infty}y^2(t) dt \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
& = \ E_x \ + \ E_y \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
\end{align} $$

$E_x$ and $E_y$ are the energies of $x(t)$ and $y(t)$ and not functions of $\tau$.  so $Q_{xy}(\tau)$ will be minimum with a $\tau$ that makes the integral on the right maximum.  we'll give that integral a name:

$$ R_{xy}(\tau) \triangleq \int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt $$

and we see that

$$  Q_{xy}(\tau) = E_x + E_y - 2  R_{xy}(\tau) $$

so so $Q_{xy}(\tau)$ is minimum when $R_{xy}(\tau)$ is maximum.  in fact, with a scaling constant of $2$ tossed in there $Q_{xy}(\tau)$ and $R_{xy}(\tau)$ are upside-down mirror images of each other.

we call $R_{xy}(\tau)$, the cross-correlation function of $x(t)$ and $y(t)$ and is the inner product (if we're pretending we're living in a Hilbert space) of $x(t)$ and $y(t+\tau)$.  if we're dealing with complex $x(t)$ and $y(t)$ (we're not, i hope) then $y(t)$ should be complex-conjugated.  where $R_{xy}(\tau)$ is maximum is the time offset that $x(t)$ and $y(t)$ are most similar to each other.

the above can be rearranged:

$$ R_{xy}(\tau) = \frac{1}{2} \left( E_x + E_y - Q_{xy}(\tau) \right) $$

or

$$ \frac{R_{xy}(\tau)}{\frac{1}{2}(E_x + E_y)} =  1 - \frac{Q_{xy}(\tau)}{E_x + E_y} $$

so the cross-correlation is getting normalized by the average energy of $x(t)$ and $y(t)$.  we're always subtracting a non-negative number from 1 so that normalized cross-correlation is never larger than 1

and pick the $\tau_0 \ne 0$ where $R_{xy}(\tau_0)>R_{xy}(\tau)$ for all other $\tau$ and calculate that normalized cross-correlation and you have a measure of similarity with

$$ -1 \le \frac{R_{xy}(\tau_0)}{\frac{1}{2}(E_x + E_y)} \le 1 $$

if it's -1, they are identical but with opposite polarity.  if it's 1, they're totally similar (at that offset of $\tau_0$).


STFT - Plotting is not a true representation I might be going crazy here. Following my question: here

I did some thinking and probably understood that my STFT algorithm was not performing correctly.

I followed this algorithm:

-> Get 256 of data

-> Apply windowing function

-> Get the FFT

-> get the first 128 points

Here is my result. I have used NFFT = 256 which should give me 128 points. I firstly calculated the FFT (after the windowing function was applied) and only counted the first 128 values of this. I.e. FFT_RESULT / 2

Whereas my last output was the following:

I don't quite understand why the values along the X axis are particular that high?

EDIT:

Ok, so following the advice, I did the following:

1) Created a 2D vector for each of the "256 resulting FFT values"

for(size_t i=0; (i < signal.size()); i++)
{
    std::vector<double> temp_vector;
    // store blocks of 256 inside temp_vector

    // Compute Hanning Window

    // Compute FFT

    std::vector<std::vector<complex> stft;

    stft.push_back(temp_vector);
}

This gave me a total of 74 blocks within the 2D vector.

I then took the absolute value of the REAL values and stored these inside a double, read these in and with pylab library I used imshow() and got the following result:

However, If I plot normally, I get the following:

I think my calculations are correct, I just don't think I'm plotting them right using the spectrogram!!

Regarding your first 2 plots: for FFT results, the X axis is usually labeled with the FFT_bin_index * sample_rate / FFT_length, and the Y axis with either the magnitude (absolute value of the complex bin value), or the log magnitude of the FFT result bin.


Cross-correlation and auto-correlation of an image using $\tt fft$ I'm working on a project where I have to find the auto-correlation and cross- correlation of types of pixels in an image:

I am using MATLAB and my code for the autocorrelation is:

fftimage = fft2(image);
congfft = conj(fftimage);%conjugate
ans1 = fftimage*congfft;
ans2 = ifft2(ans1);
ans3 = abs(ans2);
ans3 = fftshift(ans3);
image(ans3,'CDataMapping','scaled')

for the cross-correlation

 ans1 = fftimage*congfft2; %congfft2: conjugate the other type of pixel.

I don't understand clearly the outputs of the code. Here is the autocorrelation output of the yellow pixels:

For 1D signals, the autocorrelation of a rectangular pulse will be a triangular pulse.
For your image, the autocorrelation will be a triangular pulse in both directions.
The image below shows:

The original image
The autocorrelation without zero padding
The autocorrelation with zero padding
The autocorrelation with zero padding, viewed as a mesh.

For your example image, the zero padding is not required because the support of the non-zero pixels is way less than the image size. For non-zero pixels covering the whole image, you'll need to zero pad when using the FFT.
Your image does not appear to show the gradual change. I suspect you will need to normalize the pixel values manually. It's not clear to me that 'CDataMapping','scaled' does the right thing (it may, I'm just not familiar with it).

what does the output value in each pixel mean for the auto and the cross-correlation?

As calculated, not much.  It would mean a little more if the mean was subtracted from the image before calculation. As it is, the  most you are seeing is the auto-correlation of the mean pixel value (128?) over a limited range.
If the mean were subtracted, then it would tell you about how each pixel influences the pixels around it: if the autocorrelation values are close to zero, then that means "not much". If the autocorrelation values are high, then it means "greatly".

R Code Below
# 32079

par(mfrow=c(2,2))
image <- array(0,c(20,20))
image[8:13,8:13] <- 1
image(image, col= grey(seq(0, 1, length = 256)))
title("Original")

fftimage <- fft(image)
congfft <- Conj(fftimage)
ans1 <- fftimage * congfft
ans2 <- Re(fft(ans1, inverse=TRUE))/20/20
image(ans2, col= grey(seq(0, 1, length = 256)))
title("Without zero padding")

imagePadded <- array(0,c(40,40))
imagePadded[11:30,11:30] <- image
fftimagePadded <- fft(imagePadded)
congfftPadded <- Conj(fftimagePadded)
ans1a <- fftimagePadded * congfftPadded
ans2a <- Re(fft(ans1a, inverse=TRUE))/40/40
image(ans2a/max(ans2a), col= grey(seq(0, 1, length = 256)))
title("WITH zero padding")

persp(ans2a, theta=30)
title('Mesh plot WITH zero padding')


What does "x(n) is complex" mean? (DSP) I am new to EE area and I met a question in my homework:

Can someone help me a little bit? I think it is quite easy!

Thanks!!

It means that $x(n)$ and $y(n)$ can be complex numbers rather than just real numbers. As a result there is no conjugate symmetry in the Frequency domain for these signals. Note that part (c) the asterisk in $x^*(n)$ usually means to take the complex conjugate. If $x(n)$ and $y(n)$ were real, i.e. not complex, then (b) and (c) would be equivalent.


Transfer function intuition What is the meaning of the transfer function of a filter? Please explain intuitively with an example if possible.

Normally, in electrical enginnering, we apply the term "transfer function" and "filter" to an operation that belongs in the class we call Linear Time-Invariant systems (LTI).

Sometimes you might read something from someone where they apply either terms to a nonlinear operation.  A 'filter" in statistics usually means a different thing than what EEs mean.  So the use of "filter" is deprecated in this EE context.

Now either a linear or nonlinear operator will take an input signal $x(t)$ or $x[n]$ and will transform it into an output signal $y(t)$ or $y[n]$. The algorithm or rules or math of that operation is what defines what the operator is, linear or nonlinear, time-variant or time-invariant, deterministic or stochastic.

If the operator is LTI, then I would call it either an "LTI system" or a "filter".  For all practical purposes, what defines an operator as linear is

$$\begin{align}
\text{if } \quad y_1(t) &= \operatorname{LTI}\Big\{ x_1(t) \Big\} \\
\text{and } \quad  y_2(t) &= \operatorname{LTI}\Big\{ x_2(t) \Big\} \\
\\
\text{then } \quad  y_1(t)+y_2(t) &= \operatorname{LTI}\Big\{ x_1(t)+x_2(t) \Big\} \\
\end{align}$$

and what defines it as time-invariant is

$$\begin{align}
\text{if } \quad y(t) &= \operatorname{LTI}\Big\{ x(t) \Big\} \\
\\
\text{then } \quad  y(t-\tau) &= \operatorname{LTI}\Big\{ x(t-\tau) \Big\} \quad \text{for any real } \tau \\
\end{align}$$

Now, if your operator satisfies those two conditions, I will call it a "filter".

When a system is LTI, then we can apply all sorts of math on the operator and the signals going in and coming out of the operator.  This includes the Fourier Transform and the Laplace Transform (for discrete-time signals, $x[n]$ and $y[n]$, this would be the Discrete-Time Fourier Transform and the Z Transform, instead).

Probably, the first result we get from Linear System Theory (also called "Signals and Systems" for some textbooks), is that the output $y(t)$ can always be related to the input $x(t)$ by the operation of convolution:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) x(t-u) \ \mathrm{d}u \\
      &= \int\limits_{-\infty}^{\infty} h(t-u) x(u) \ \mathrm{d}u \\
\end{align}$$

$h(t)$ is the impulse response, and in a sense, fully defines the operation of the filter or LTI system.  I.e. if you know how your LTI system responds to a Dirac impulse as an input, you know how your LTI system will respond to any input.

Then, whether we use the Fourier Transform

$$\begin{align}
 X(f) &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ \mathrm{d}t \end{align}$$

or the (bilateral) Laplace Transform

$$\begin{align}
 X(s) &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-s t} \ \mathrm{d}t \end{align}$$

then the above convolution operation in the time domain becomes a simpler multiplication operation in the frequency domain (or sometimes the "transform domain"):

$$\begin{align}
 Y(f) &= H(f) \, X(f) \\
\\
 Y(s) &= H(s) \, X(s) \\
\end{align}$$

where $X(\cdot)$ and $Y(\cdot)$ and $H(\cdot)$ are all defined the same way from $x(t)$, $y(t)$ and $h(t)$.

The thing that multiplies the (transform of the) input to become the (transform of the) output is the transfer function.  That means $H(f)$ or $H(s)$ transfers the input $X(f)$ or $X(s)$ to the output $Y(f)$ or $Y(s)$.  But that is only in the frequency domain.

In the time domain, the convolution operator (convolving the input $x(t)$ with the impulse response $h(t)$) is what your filter does.

Either the transfer function ($H(f)$ or $H(s)$) or the impulse response ($h(t)$) are sufficient to fully define your LTI system (or filter) from the POV of the input to output relationship.  I.e. two different systems with identical impulse responses $h_1(t) = h_2(t)$ are indistinguishable from each other if both are placed in black boxes.  The same can be said about their transfer functions $H_1(s) = H_2(s)$.  If their transfer functions are identical, from just looking at the inputs and outputs, we cannot tell any difference between the two LTI systems.


Does the DTFT of $\frac{u[n-1]}{n}$ exist? I have started learning DSP on my own and I have this doubt. I have done some googling but haven't found an answer. I hope that someone here would give the answer. It will be of great help.

Note that the sequence

$$x[n]=\frac{u[n-1]}{n}\tag{1}$$

is in $\ell^2(\mathbb Z)$ because

$$\sum_{n\in\mathbb{Z}}|x[n]|^2=\sum_{n=1}^{\infty}\frac{1}{n^2}=\frac{\pi^2}{6}<\infty\tag{2}$$

but it is not in $\ell^1(\mathbb Z)$ since

$$\sum_{n\in\mathbb{Z}}|x[n]|=\sum_{n=1}^{\infty}\frac{1}{n}=\infty\tag{3}$$

We know that for sequences in $\ell^1(\mathbb{Z})$ the DTFT $X(e^{j\omega})$ exists and it is a continuous function of $\omega$. However, it is not necessary for a sequence to be in $\ell^1(\mathbb{Z})$ for the DTFT to exist. We can extend the definition of the DTFT to sequences in $\ell^2(\mathbb{Z})$, which implies that the DTFT converges in a mean square sense (and not uniformly), and that $X(e^{j\omega})$ may have discontinuities.

In that extended sense, the DTFT of $(1)$ exists, just as it exists for the ideal discrete-time differentiator or the discrete-time Hilbert transformer, the impulse responses of which also only decay as $1/n$.

The DTFT of $x[n]$ is

$$X(e^{j\omega})=\sum_{n=1}^{\infty}\frac{e^{-jn\omega}}{n}\tag{4}$$

The series $(4)$ is the complex Mercator series

$$\sum_{n=1}^{\infty}\frac{z^{n}}{n}=-\log(1-z),\qquad |z|\le 1,\quad z\neq 1\tag{5}$$

with $z=e^{-j\omega}$. It is the Taylor series for $-\log(1-z)$ and it converges for all $|z|\le 1$, $z\neq 1$. From $(5)$ we can immediately write down the $\mathcal{Z}$-transform of $x[n]$:

$$X(z)=\sum_{n=1}^{\infty}\frac{z^{-n}}{n}=-\log(1-z^{-1}),\qquad |z|\ge 1,\quad z\neq 1\tag{6}$$

Since $(6)$ also converges for $|z|=1$ except for $z=1$, we can obtain the DTFT directly from $(6)$ by setting $z=e^{j\omega}$. For $e^{j\omega}=1$, i.e., $\omega=0$, the DTFT has a singularity.

  $$X(e^{j\omega})=-\log(1-e^{-j\omega})\tag{7}$$

Equation $(7)$ agrees with WolframAlpha's result.
Also take a look at this related answer.

Let's visually check if it is plausible that the series $(4)$ really converges to $(7)$. I've computed the sum over the first $100$ terms of $(4)$, and I've plotted the real and imaginary parts of that partial sum together with the real and imaginary parts of $(7)$:

From the figure, it appears that the series $(4)$ indeed converges to $(7)$ in a mean square sense. This is of course no proof but a simple sanity check.


Combine FFT's of shorter length than sample data to get spectrum of all data I have a set of sample data, a logarithmic swept sine as it happens, that is 1 second in duration sampled at 100kHz giving 100k samples. I want to calculate the frequency spectrum of the complete data however the FFT library function available only supports 4096 samples maximum.

Is it valid to sequentially FFT 4096 samples and sum the resulting data. Thereby ending up with 4096/2 frequency points that show the frequency content of the entire 100k samples?

I guess it is necessary to window each block of 4096 samples to be FFT'd but this would mean that the portions of the frequency content would be reduce in amplitude artificially.

How does one handle the phase information?

Yes you can do it, but you will loose spectral resolution. The method you mention is a variant of the Welch's spectogram estimation algorithm which relies on the block based estimation of power spectrum of a long sequence in terms of shorter blocks within it.

Given that you have a data of $100$k in length which was sampled at $F_s = 100$k, then a full length DFT after rectangular windowing would provide you a spectral resolution of about $1$ Hz to $3$ Hz, where as a window of $4096$ samples would yield about $25$ Hz to $75$ Hz of spectral resolution.

If this resolution is enough for your purposes, then you can (and indeed you should better) perform the block based averaged spectral analysis in order to trade-off between the reduced spectral resolution and an increased estimation accuracy (decreased estimation variance)


Proof of Bedrosian's theorem I was looking for a straightforward proof for Bedrosian's theorem which says the Hilbert transform of the baseband signal  times the passband signal is the original baseband signal times the Hilbert transform of the passband signal in which by assuming $x(t)$ is a passband signal and and $m(t)$ is a baseband signal, then we have :
$$c(t)=\widehat{x(t)m(t)}=\widehat {x(t)}m(t)$$
Now for the sake of proof since we know Fourier transform of Hilbert  equals to  :
$$\hat X(f)=-j\operatorname{sgn}(f)X(f)$$
Then we say :
$$\hat c(f)=-j\operatorname{sgn}(f) [x(f) \star m(f)]$$
Now can we say( because $x(t)$ is a passband signal) $X(f)$ is only in positive f part and $\operatorname{sgn}(f)$ value is constant and bring it in? and then we have:
$$[-j\operatorname{sgn}(f)x(f)] \star m(f)=\widehat {x(f)}\star m(f)$$
And so we can say
$$\hat c(f)=\widehat {x(f)}\star m(f)$$
$$\hat c(t)=\widehat{x(t)m(t)}=\widehat {x(t)}m(t)$$

One easy way to see the validity of the theorem is to consider the corresponding analytic signals. Let a signal $s(t)$ be defined by
$$s(t)=m(t)c(t)\tag{1}$$
where $m(t)$ is a low pass (message) signal, and $c(t)$ is a bandpass (carrier) signal. Note that we don't need to assume that $c(t)$ is sinusoidal. Let $\hat{x}(t)$ denote the Hilbert transform of $x(t)$, and $\hat{c}(t)$ is the Hilbert transform of $c(t)$.
The equality claimed by Bedrosian's theorem
$$\hat{s}(t)=m(t)\hat{c}(t)\tag{2}$$
is equivalent to
$$s_a(t)=m(t)c_a(t)\tag{3}$$
where $s_a(t)$ and $c_a(t)$ denote the analytic signals corresponding to $s(t)$ and $c(t)$, respectively. Analytic signals are complex-valued signals having only positive frequency components. The analytic signal is given by the original signal plus its Hilbert transform as the imaginary part:
$$s_a(t)=s(t)+j\hat{s}(t)\tag{4}$$
In the frequency domain, the signal $m(t)c_a(t)$ is centered at the (positive) center frequency of the bandpass signal $c(t)$. Its bandwidth equals the bandwidth of $c(t)$ plus twice the bandwidth of the lowpass signal $m(t)$. It has only positive frequency components if the spectra of $m(t)$ and $c(t)$ do not overlap. In that case it equals the analytic signal $s_a(t)$, as claimed in Eq. $(3)$. Otherwise, $m(t)c_a(t)$ has also negative frequency components and $(3)$ is not satisfied.
Note that $(2)$ and $(3)$ are also valid if $c(t)$ is a highpass signal, as long as the spectra of $m(t)$ and $c(t)$ don't overlap.
Also take a look at this related answer and the figure therein.


Filter's which requires physical model This can be very general question, but I am trying to under the basics. I want to know the difference between the filters which require a physical model in order to perform the filtering  and the ones which doesn't. I read someone one such type of filter is Luneberg. But I could not find anywhere about them. If anyone could explain difference between these two for understanding purpose, it'd be great.

PS- this is not a homework. I am just trying to understand.

After reading your comments, the real question you're asking is:

"What's the difference between an estimator and filter?"

An estimator is a system that takes in inputs and attempts to estimate a signal (or state) given knowledge of the system. A simple example is system that takes in positional data and calculates the corresponding acceleration. In this case, your input data is position, the system is the second derivative with respect to time, and the output is acceleration. You've just estimated the acceleration by feeding it through a model.

A filter is just a system that takes an input and attempts to modify the frequency content of a signal. A Butterworth filter is a filter, not an estimator.

So the real question you have to ask your is: "What am I trying to do?"

A problem where you're trying to estimate something, you use an estimator (i.e. a system that produces an output using some physical model). A problem where you trying to change the frequency content of a signal, you use a filter. The filter doesn't care about the underlying process or signal, it's just setup to modify the frequency content of the signal.

Mathematically, a filter and estimator are the same thing (put signal in, get signal out), but they typically have different physical meanings and are used to achieve different things.

Edit: It also seems that your question may have been inspired by the 'Kalman Filter'. Even though this contains the word 'filter', its still considered more of an estimator. However, depending on how you set up the problem (i.e. how much uncertainty you have in your model and how much uncertainty you have in your feedback), the system will experience different frequency response characteristics.


Basic difference between Multi-scale and Multi-resolution images While working on Convolutional Neural Network (CNN) I saw many research articles working on the Multi-scale images. In the literature of image processing I saw there is also Multi-resolution analysis of image which is at different level of images. Whereas in Multi-scale it seems that size of images are not changed instead it gets more blur and edges are getting more clear of the objects.
I am not using any wavelet transform. So in terms of general images how can we explain to beginners about the basic difference between Multi-scale and Multi-resolution images?

In scale space the square of width of blurring kernel considered as the scale (variance of Gaussian and not its Standard deviation because if you sequentially blur an image with 2 Gaussian kernel variance of effective blurring kernel will be sum of variance of each kernel) but in multi-resolution analysis the resolution defined by the scaling parameter.


Why is the Kalman gain constant when position and velocity change? I used the measured position to estimate the velocity using a Kalman filter.

But in simulations, the Kalman gain changes quickly and then remains constant when position and velocity continue to change. For example, position and velocity change at times  0 to 0.5(s) and 3 to 4(s), respectively. But the Kalman gain only changes  0 to 0.1(s) and then remain constant.

My question is: Why is the Kalman gain constant when position and velocity change?

The Kalman Gain $K_k$ for a discrete, linear system is computed using the state covariance matrix $P_k$, measurement matrix $H_k$ and measurement noise covariance matrix $R_k$:

$$
K_k = P_kH_k^T(H_kP_kH_k^T+R_k)^{-1}
$$

If the measurement matrix and the measurement noise covariance matrix are constant and the state covariance matrix converges to a steady-state value $P_\infty$, the Kalman Gain also converges to a steady-state value $K_\infty$, regardless of actual measurements or state.


derivation of discrete Fourier Transform I'm learning DSP on my own. Found a free online book. Currently reading a chapter on Discrete Fourier Transform. Right out of the gate I can't understand several things. I'm pretty good with math and have a decent understanding of complex numbers and their representation with the help of $e^{j\phi}$. But this got me stumbled.

Below is the paragraph I can't understand and the questions pertaining to it.

  We now develop a Fourier representation for finite-length signals; to
  do so, we need to find a set of oscillatory signals of length $N$
  which contain a whole number of periods over their support. We start
  by considering a family of finite-length sinusoidal signals (indexed
  by an integer $k$) of the form

  $$ w_k[n]=e^{j \omega_k n}, \space n = 0,...,N-1 $$ where all the $\omega_k$’s are
  distinct frequencies which fulfill our requirements. To determine
  these frequency values, note that, in order for $w_k[n]$ to contain a
  whole number of periods over $N$ samples, it must conform to

  $$ w_k[N]=w_k[0]=1 $$ which translates to $$ (e^{j \omega_k})^N = 1 $$

  The above equation has $N$ distinct solutions which are the $N$ roots
  of unity $e^{j\frac{2\pi m}{N}}$, $m = 0,…,N - 1$

which contain a whole number of periods over their support.
What does "periods of their support" mean?
note that, in order for $w_k[n]$ to contain a whole number of periods over $N$ samples. Why does $w_k[n]$ must contain a whole number of periods over $N$ samples? Where does this requirement come from?
it must conform to $w_k[N]=w_k[0]=1$. Why? Where does this requirement come from? I do understand $w_k[N]=w_k[0]$, that's because it is periodic, what I don't understand is $=1$ portion. Why at sample $0$ and sample $N$ it must be equal to 1?
which translates to $(e^{j \omega_k})^N = 1$. How does $w_k[N]=w_k[0]=1$ translate to $(e^{j \omega_k})^N = 1$?
The above equation has $N$ distinct solutions which are the $N$ roots of unity $e^{j\frac{2\pi m}{N}}$, $m = 0,…,N - 1$. Can someone demonstrate the finding of these roots?

The Discrete Fourier Transform is one-and-the-same as the Discrete Fourier Series.  The DFT (and explicitly DFS) maps one periodic sequence, $x[n]$ of period $N$:

$$ x[n+N] = x[n] \qquad \forall n \in \mathbb{Z} $$

to another periodic sequence, $X[k]$, having the same period $N$:

$$ X[k+N] = X[k] \qquad \forall k \in \mathbb{Z} $$

and the inverse DFT (or "iDFT") maps it back.

I have a pretty rigid Nazi-like position about it (see the other answer), but I don't think that the math cares.  Now if you wanna connect this to the more common (but often misinterpreted) derivation, you can look at a previous answer of mine.

The term "support of $x[n]$" simply means the domain of values of $n$ in which $x[n]$ (presumed to be non-periodic and to have "finite support") is potentially non-zero.  It could also mean the duration of $x[n]$.  But when you periodically extend $x[n]$, then that sequence has no finite support.  $x[n]$ is defined and potentially non-zero everywhere (or all integers).


FIR Filter matlab command i'm trying to understand the fir1 filter but i still don't get it. For example here i got an audio signal that i consider noise, i'm passing it through a low pass filter.

n = 100000
fs = 11025
handles.noise = wavrecord(n, fs, 'double');
nfilt = fir1(11,0.4);
fnoise = filter(nfilt,1,handles.noise);

Why did i need to pass the noise signal to a low pass filter in the first place? and why is the number of coefficients 11? also the 0.4 shouldn't it be a range of values?

I know these are a lot of questions but when i get to understand them, ill be grateful to everyones help.

Matlab FIR-1 command designs FIR filters with at least two given parameters. P1: The order of the filter, P2: cutoff frequency, assuming a default lowpass filter. Order of the filter N determines a number of properties of the resulting filter such as its comptational complexity, transition bandwidth, passband and stop band ripples. In general large N gives better (closer to ideal) frequency response but worse computational load. Also order of filter being even or odd will have effect on zero-locations being at Nyquist frequency which might be important for some applications with high/low pass filters...

The second parameter, cutoff frequency, is basicaly determined by the application.


How to find DC bias of periodic signal from samples that include a non-integer number of periods? I'm using an ADC to sample a periodic signal and I want to calculate the DC bias of the signal from the samples.  The obvious solution is to take the average of the samples.  This works fine if there are an integer number of periods of the signal in the sample set.  But if the sample set contains an additional fraction of a period then those extra samples can throw the average off a bit.

The signal frequency can vary (from one sample set to another) and is independent of the sampling rate.  The sampling is asynchronous to the signal so sampling could begin anywhere in the period of the signal.  Assume there is at least a couple periods of the signal in the sample set although it could be many more.

Do you know of a better way to calculate the DC bias of the signal from a sample set that contains a non-integer number of periods of the signal?

Update:
For my application this is not a real-time calculation.  My device collects a sample set occasionally and I'm trying to do post-processing of the samples for various measurements.

I don't know the period of the signal by any other information.  All I have is the sample set.

This is admittedly an XY problem.  I want to know the DC bias so I can remove it from the samples (and then I plan to integrate the signal).  I'm hoping for a way to improve my original implementation (which is easy for me to understand) before trying something more complicated (for me).

Low-pass filter with DC gain = 1 (or 0 dB).

i s'pose, if you want to be sophisticated and want a rock-solid DC that might vary very little, then the output of your LPF can go into a median filter.


What is an intuitive understanding of the difference between the minimum required bandwidth and actual bandwidth? When using for example, polar signaling of full width rectangular pulses, the bandwidth is equal to the pulse rate $Rb$, which is also the bit rate.

$B=R_b$ for full width rectangular pulses.

But also we can describe the minimum required bandwidth for transmission of this data.

$B_{min}=\frac{R_b}{2}$

If I am sending data over a channel using polar signaling of the full width rectangular pulses, surely I cannot achieve a bandwidth of $B_{min}$?

How would I go about converting my data to achieve this theoretical minimum bandwidth? I'm looking for an intuitive understanding of minimum bandwidth versus the actual bandwidth of my data.

We must be careful in how we define bandwidth. To maintain truly rectangular pulses in the time domain, infinite bandwidth is required, as in the frequency domain the spectrum would be a Sinc function that extends to infinity (as the Fourier Transform of the rectangular pulse). The requirement for infinite bandwidth is intuitively consistent with a rectangular pulse that can change instantaneously from a low level to a high level (infinite bandwidth is required to transition that fast!). Note for now that the Sinc function has nulls at 1/T and integer multiples of that, where T is the width of the rectangular pulse.

As we constrain the bandwidth in frequency, the edges of the signal will likewise slow in transition from low to high. A good rule of thumb related to this (specific to first order systems but generally a good estimate) is:

$$t_r = 0.35/BW$$

Where $t_r$ is the 10% to 90% rise time in seconds, and BW is the bandwidth in Hz.

Ultimately for a digital signal the fastest bandwidth is required when we send the sequence 1 0 1 0 1 0 1 0 ....   We note that for rectangular pulses this would be a 50% duty cycle square wave which according to the Fourier Series Expansion is decomposed in frequency to the odd harmonics of the fundamental sinusoidal frequency that fits in the 1 0 1 0 pattern.  This is still the Sinc function we described above, just when the time domain signal repeats the only frequency components that can exist are at integer multiples of the repetition rate, and therefore those frequencies become samples of the Sinc function.  As described above, the nulls would exist at 1/T, where T in this case is $1/R_b$ which is also the location of the 2nd harmonic that we already reasoned is nulled (only odd harmonics exist for the 50% square wave). The first harmonic is at $R_b/(2)$, and if we were to represent it as a double sideband signal (positive and negative frequencies) then the minimum bandwidth required to contain both would be $R_b$. This would result in a sinusoidal waveform (sine or cosine) but would be sufficient to represent our 1 0 1 0 pattern without data loss (we recover with a simple threshold at y=0).

However, we can also represent the same sequence as a single sided signal (positive or negative frequencies only) which would then only require half that bandwidth, which is how you came to $R_b/2$. This would be consistent with a use of Single Sided Modulation.  (Note we don't require any bandwidth if the pattern was always 1 0 1 0 as it would be represented by a single frequency, but this is the case of the highest frequency needed, and we would also need all the lower frequencies to represent any other patterns).

Depending on your background this may be a lot to swallow, and if so I also point you to the links below that also help explain these fundamental concepts:

How do I modulate complex time domain signal to a carrier frequency?

Frequency shifting of a quadrature mixed signal


Interpolation algorithms in computer graphics I'm looking for an image-processing algorithm that uses interpolation (polynomial, spline, doesn't matter). I have to implement such algorithm to process a bitmap for my labs. It can't be image resizing, it must be something else. Could you give me some proposals?
I've tried on wikipedia and google, but it's hard to filter algorithms that make use of interpolation

PS. Sorry for the tags, I can't create new ones and couldn't find existing ones that would be suitable for my question.

Some ideas:

rotate the image in 2D
compute the Radon transform.  This is used in tomography applications.
Other 2D effects, like swirling or adding ripples
Map the 2D image onto the surface a 3D shape, such as a sphere
Simulate the effect of a lens, or refraction by some other transparent 3D object in front of your bitmap
Simulate gravitational lensing by a black hole

2D image rotation is what comes to mind first as the simplest application of 2D interpolation after rescaling.  The obvious approach is to simply apply the well-known coordinate transformation; an alternative trick is to decompose the rotation into a sequence of three sheers so that only 1D interpolation is needed.

The 3D examples are more involved, since they will probably require implementation of raytracing.


discrete Haar wavelet transform, fast and efficient method? I'm working on my own implementation of the discrete Haar wavelet transform, I understand the wavelet theory and how to construct the Haar matrix of size N to perform the transform, but obviously there is a problem using the Haar matrix in application - it's simply too big.

I am working on an application that applies the Haar transform to an audio signal. If I sample the signal at 44.1 kHz, even a one second recording would require a 2^16 by 2^16 Haar matrix to do the transform in one step, which is obviously impractical and a hardware constrained machine such as a phone wouldn't have the capability to hold such a matrix in memory.

The other method I've seen used is that a 2X2 Haar matrix is applied to the entire signal iteratively, and the results are stored in two arrays - one array holding the "average" Haar coefficients (first element of the output vector) and the other holding the "difference" coefficients (second element of the output vector). The process is then repeated over the "average" coefficients - as these coefficients are essentially the result of a lowpass filter.  Each time the process is repeated the number of elements needed to process in the following step is halved, until only one lowpass coefficient is left. This seems fine, and definitely works but when I started to think about it, it just seems that it would be really slow.

My main question is, what's a good way to implement a fast and efficient haar transform? Or a practical way to apply one of these two methods?

PS: I've been learning about this completely on my own, so if I made some mistakes in my explanation or way of thinking about this stuff let me know.

PSS: I've never done any kind of audio processing before, so if you know about some other filters that I should apply to the raw signal before doing a wavelet transform let me know!

PSSS: I know that the Haar wavelet may not be the best for this type of signal processing, but when I tried to learn about using other DWTs such as the Daubechies wavelets, the literature seemed very confusing, or at least was directed at more advanced readers. If anyone could point me in the direction of how to implement other DWTs, that would be great.

The Haar wavelet is actually a part of the Daubechies wavelet, for the case D=2. There's some example code on wikipedia that shows the Daubachies transform.

The Haar transform is just a low pass filter combined with a high pass filter, with the coefficients being placed in the first and second halves of the signal. Then, this iteratively (or recursively) keep going.

The low pass filter is just the sum of two components while the high pass filter is their difference.

For the Haar wavelet case, you can check out my vectorized code on Github. Here, idwt2 stands for "Inverse Discrete Wavelet Transform 2D."


Finding the output of a system where the input is a sum of complex exponentials So, I have to find $H\{ x(t)\})$ (which is an LTI system), where $$x(t) = \sum_{k=0}^{\infty} a_ke^{ \ jw_kt}$$ and where the impulse response of the system is given by: $$h(t) = \frac{\delta(t+\tau)-\delta(t)}{\tau}$$

Here's is my attempt. I know that $$H\{x(t)\} = \int_{-\infty}^{+\infty}x(t-u)h(u)\ du \\ =\frac{1}{\tau}\left(\int_{-\infty}^{+\infty}x(t-u)\delta(u+\tau)\ du\ \  - \ \ \int_{-\infty}^{+\infty}x(t-u)\delta(u)\ du\right)$$
but I don't know where to go from there. If you could help, that would be really nice, thank you.

Okay... I managed to do it on my own. Here's the solution if you're interested.

We know that $$H\{x(t)\} = \sum_{k=0}^{\infty}a_k \  H\{e^{jw_kt}\}$$ since $H$ is a linear system.

We also know that $H\{e^{jw_kt}\}  = F(w_k) \ e^{jw_kt}$, where $F(w_k)$ is the frequency response of the system.

By definition, we have that $$F(w_k) = \int_{-\infty}^{+\infty}h(t) \ e^{-jw_kt} \ dt \\ = \frac{1}{\tau} \left(\int_{-\infty}^{+\infty}\delta(t+\tau)e^{-jw_kt} dt \ - \ \int_{-\infty}^{+\infty}\delta(t) e^{-jw_kt} dt\right)\\ =\frac{1}{\tau} \left( \int_{-\infty}^{+\infty} \delta(u) e^{\ jw_k(\tau-u)} \ du \ - \
\int_{-\infty}^{+\infty} \delta(t)e^{\ jw_k(0-t)}\ dt \right) \\ = \frac{e^{ \ jw_k\tau} - e^0}{\tau}$$

Therefore we can simply write $$H\{x(t)\} = \sum_{k=0}^{\infty} \frac{e^{\ jw_k\tau}-1}{\tau} \ a_ke^{\ jw_kt}$$


How to take the FFT in a limited range of frequencies? When using the fft function in Matlab to take the Fourier transform, one can specify the number of FFT points, nfft, and the command is written fft(x,nfft). In this case, what is the frequency vector?

How can the Fourier transform be calculated with the same nfft, but over a narrower range of frequencies? How would the frequency vector be calculated?

If your sampling frequency is $f_s$, and you calculate the FFT with $N$ samples, then the FFT calculates the signal's spectrum at $N$ frequencies from $-f_s/2$ to $f_s/2-f_s/N$, inclusive. In Matlab notation, this would be written as:

freq_vector = (-N/2:N/2-1)./(f_s/N);

(Note that you have to shift the spectrum returned by the FFT operation to correspond to this vector; in Matlab, you would use fftshift.)

I don't know of any FFT algorithm that lets you restrict the frequency range. I would suggest using Goertzel's Algorithm to find the spectrum at the particular frequencies you're interested in. You may also use a bandpass filter to eliminate the unwanted frequencies before performing the FFT.


How to figure out number of multiplications in rational sampling I am trying to understand the following:
Consider a system implementing a rational sampling rate change by $\frac{5}{7}$: for this, we cascade upsampler by 5, a lowpass filter with cutoff frequency $\frac{\pi}{7}$ and a downsampler by 7. The lowpass filter is a 4th order Butterworth filter with transfer function
$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 - a_1 z^{-1} - a_2 z^{-2} - a_3 z^{-3} - a_4 z^{-4}} $$
Assume that the input works at a rate of 1000 samples per second. What is the number of multiplications per second required by the system? Assume that multiplications by zero do not count and round the number of operations to the nearest integer.
I would like to know how to approach this problem and get a general idea. I am taking an online course on my own and there is very little reference material, so I was trying to figure out a good way to approach this problem, but I am stuck with the logic behind it.

That depends a fair bit on how exactly you implement the conversion but if you follow the textbook approach it would probably look like this:

The generic order of sampling rate change is "upsample", "filter", "downsample".
That means that you need to run the filter at the upsampled rate, i.e 5000
Your IIR filter has 9 non-trivial coefficients, i.e. you need to do 9 multiplies per input/output sample
In total that's 9*5000 = 45000 multiplications per second.

That's a lot and that's whey rational sample rate converters often implemented as polyphase FIR filters. Using an IIR requires you to calculate all output samples of the filter although you will be throwing 6 out of 7 away.


Understanding Integration Scale in Multi-Scale Harris Corner Detector I would like to improve a simple Harris corner detector by adding more scales.

I don't need scale-invariant interest points, just to detect more points by processing more scales.

The Wikipedia says there are two sigmas involved: a derivative scale $\sigma_{D}$ and an integration scale $\sigma_{I}$.

So after computing derivative images on a given level, I smooth them with Gaussian filter of size $\sigma_{D}$, then for each location $(x,y)$ in the image I have the matrix:

$H(x,y)=\begin{pmatrix}I_{x}^{2} & I_{x}I_{y} \\ I_{x}I_{y} & I_{y}^{2} \end{pmatrix}$

where $I_{x}$ and $I_{y}$ are the smoothed images (with $\sigma_{D})$ on the respective location.

How do I apply the integration scale $\sigma_{I}$ now?

Do I have to store separate images $I_{x}^{2}$, $I_{y}^{2}$ and $I_{x}I_{y}$ and smooth these with $\sigma_{I}$ before computing the matrix and Harris corner response from it?

I think the first smoothing, by $\sigma_D$, is only done to get more stable derivatives whereas in the second step the convolution by a Gaussian with $\sigma_I$ is done to establish the 'scale-space' in which the operator is applied.

Ignoring $\sigma_D$, this looks like this in Matlab:

dx = [-1 0 1; -1 0 1; -1 0 1];  % Simple mask for derivative
Ix = conv2(im, dx, 'same');     % Convolve against image
Iy = conv2(im, dx', 'same');    % Again for y-direction with transposed mask

% Gaussian filter, as defined by Kovesi
% www.csse.uwa.edu.au/~pk/Research/MatlabFns/index.html
sigma = 1;
g = fspecial('gaussian', max(1, fix(6*sigma)), sigma);

% Determine smoothed squared image derivatives via convolution
Ix2 = conv2(Ix.^2, g, 'same');
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');

% Metric to define the corner score
cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;    % Original measure
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);    % Kovesi measure


Are reconstruction filter always needed? in mixed (digital and analogue application) where you have to convert a signal from continuous to digital time domain and then back again to continuous time an anti-aliasing LPF is needed and in general also a LP reconstruction filter.

Are there some applications (if so what are) in which the reconstruction filter is not needed or not strictly necessary?

If the system naturally has only a low enough frequency response (say a very heavy object relative to the forces involved), and the noise tolerance is high enough, then the effects of any aliasing, due to the lack of any additional low-pass and anti-aliasing filtering, may already be rolled-off below your desired or required noise floor.


Vector length output of discrete time convolution Suppose that the impulse response of a discrete time filter is $h[i]$ where $i=0,1,2,...,N-1$ and the input sequence to the filter is $x[i]$ for $i=0,1,2,.., M-1,$ what would be the length of output vector length?

I have concluded it is $N$  because  according to discrete time convolution

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

Is this correct?

The output $y$ is convolution of $x[n]$ and $h[n]$:

$$y[n]=\sum_{l=-\infty}^{l=+\infty} x[l]h[n-l]$$
As you said, let's assume only $h[0],\cdots,h[N-1]$ and $x[0],\cdots,x[M-1]$ are non-zero.

When $0\le l\le N-1$, then we have a non-zero $x[l]$.
Now consider $h[n-l]$. As long as $0 \le n-l \le M-1$, then we have a non-zero $h[n-l]$.

add up the sides of the two inequalities to get

$$0\le n\le N+M-2$$

as the set of indices where  $y[n]$ is nonzero. There are $N+M-2+1=N+M-1$ terms.


Is there any way to make a sampling rate reduction in the frequency domain? Is there any way to take advantage of the frequency domain to make a clean/nice sampling rate reduction effect ? and how ?

Basically, deducted from my tests:

if I set to 0 the magnitude uppon a certain threshold I only get a filter(quite obvious...)
if I set to 0 one sample over 2 before the buffer transfrom it's not at all some freq domain processing ( if I do the same on the transformed buffers it's also a filter)

Is it possible to make SR reduction in the freq domain ? (actually I have some overlapped phase and magnitude arrays...)

Yes, you can do what you're looking for. Think about what happens in the frequency domain when you downsample a signal in the time domain: the spectrum "folds over" on top of itself and sums, the well-known aliasing effect. If your signal is already sufficiently bandlimited, then you can obtain the downsampled time-domain signal without computing the samples that you would subsequently throw out.

Given a signal $x[n], n = 0, 1, \ldots , N-1 $ that has DFT $X[k], k = 0, 1, \ldots , N-1$, a new signal $x_d[n]$ that is $x[n]$ decimated by a factor $D$ can be generated by:

$$
X_d[k'] = \sum_{d=0}^{D-1} X[k' + d\frac{N}{D}], k' = 0, 1, \ldots , \frac{N}{D} - 1
$$

$$
x_d[n] = IDFT(X_d[k'])
$$

The inverse DFT (and subsequently $x_d[n]$) is of length $\frac{N}{D}$. This method is sometimes referred to as "stack and add".


Am I violating Nyquist? I am using a DAC to transmit 2 different voltage levels. High and low. At 1Gsps. There is no carrier wave or anything. These are just raw 0 and 1 time samples going across a wire.

On my ADC I am also sampling at 1Gsps and I'm able to fully reconstruct the transmitted signal. In other words I can tell which sample is a 0 or 1 level.

Does this not violate the Nyquist sampling theorem? I assume if I am transmitting at 1Gsps my receiver would need to sample at 2Gsps.

Where is my understanding wrong here?

Yes you are violating! But not the Nyquist rate as a matter of fact :))

First of all Nyquist (Shannon) sampling theorem is based on the analog signal's bandwidth, and states that you must sample the signal at a rate twice the bandwidth of it. Now, what's the bandwidth of your signal? You must compute this one first.

Given a DAC (digital to analog converter) modeled as D/C block; i.e., an ideal signal reconstructor with an ideal interpolation (image rejection) filter with $h_r(t)$ whose cutoff frequency is at a half of the DAC sampling rate, then the output of the DAC given at a sampling rate of $1$G samples per second (1 GHz) will be bandlimited to $500$ MHz.

Now given a continuous-time (analog) bandlimited baseband signal $x(t)$ with a bandwidth of $500$ MHz, the required Nyquist sampling rate will be twice the bandwidth of the signal which is $1$ GHz, or stated as $1$ G samples per second.

Hence in principle (ideally) there's nothing wrong with (Nyquist is not violated by) using a DAC-ADC communication pair with the same sampling frequency. Of course due to several nonidealities (please refer to @MarcusMüller 's answer for those) you will find it quite impractical to achieve. Anyway, Nyquist is ok ;-)

But, then, you may be confusing the sampling theorem with channel capacity (information rate) for AWGN channel. So given your DAC-ADC AWGN channel system with a bandwidth of $B = F_s/2$, and given an SNR of the line in between; your maximum rate of information transfer is given by Shannon-hartley law which states that $C = B \times \log_2( 1 + SNR)$ b/s...

Assuming an SNR of $90$ dB (for about $16$ bits of an ADC,DAC with LSB noise) this means in principle you should be able to reach a data rate of $C = 500 \times 10^6 \times 15$ bits per second. This is about $7.5$ Gbps.


How do I find the peak frequency of a wav file in python? I am fairly new to python and signal processing and I was given a task to record audio for 'x' seconds and then find the peak frequency in the audio file.

So far I have successfully implemented the recording part (records as a .wav file, sample rate = 44.1 kHz) but I am unable to correctly find and output the peak frequency in that file.

The goal of this task is that if a frequency above a certain threshold is found, other actions will take place. Could anyone please assist me on how to find the peak frequency.

So far I have attempted this using the scipy fft package and some pieces of code taken from research.

from scipy import signal
from scipy.io import wavfile
from scipy.fftpack import fft, ifft,fftfreq
import matplotlib.pyplot as plt
import wave
import numpy as np
import sys
import struct

frate,data = wavfile.read('output.wav')
print(frate)

w = np.fft.fft(data)
freqs = np.fft.fftfreq(len(w))

# Find the peak in the coefficients
idx = np.argmax(np.abs(w))
print(idx)
freq = freqs[idx]
freq_in_hertz = abs(freq * frate)

print(freq_in_hertz)
print("HZ")

When audio is recorded, the above code sometimes works but however, an incorrect peak frequency is detected (for example I record something that I know has frequencies of 10kHz but it says peak frequency detected is 6kHz)

Other times when I record audio and try to find the peak frequency, i get the following error message: Exception has occurred: IndexError index 547388 is out of bounds for axis 0 with size 480000 File "C:\Users\aviso\Documents\VS Code Projects\FinalDesign\HighFreqTest.py", line 22, in freq = freqs[idx]

import sounddevice as sd
import scipy
import numpy as np
from scipy.io import wavfile as wav
from scipy.io.wavfile import write ,read
from scipy import fftpack as scfft
from matplotlib import pyplot as plt

fs = 44100
seconds = 3
myrecording = sd.rec(int(seconds * fs), samplerate=fs, channels=1)
sd.wait()
write('output.wav', fs, myrecording)

fs_rate, signal = read("output.wav")
l_audio = len(signal.shape)
N = signal.shape[0]
secs = N / float(fs_rate)
Ts = 1.0/fs_rate
t = scipy.arange(0, secs, Ts)
FFT = abs(scipy.fft(signal))
FFT_side = FFT[range(N//2)]
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])
fft_freqs = np.array(freqs)
freqs_side = freqs[range(N//2)]
fft_freqs_side = np.array(freqs_side)

volume=np.array(abs(FFT_side))
audible=np.where(volume>5)

HighestAudibleFrequency=max(freqs_side[audible])
print(HighestAudibleFrequency)

The above code records audio for 3 seconds from the microphone and outputs the highest frequency detected in the 3 seconds. Sorry for the code format, i dont know how to make it look better


Why am I seeing a negative frequency on FFT graph in Gnu Radio Companion? I have a pretty simple flow here in GRC, but for some reason I'm seeing a negative frequency mirroring the 10khz signal that I'm generating. What's causing this negative frequency to show up in the FFT graph?

The stock FFT (in GNU radio?) is a complex-to-complex transform.  Thus any positive frequency peak you see represents a complex signal (phasor) that can include both real and imaginary components.

Since your cosine (or sine) waveform is strictly real (has zero or no imaginary component), the complex FFT result also includes a negative complex conjugated mirror image.  Since complex conjugation inverts the imaginary component, the two imaginary components (of the positive and negative frequency peaks) cancel out, and thus represent your original strictly real signal (with zero imaginary content).  That is why you see a negative frequency peak.

An asymmetric result can occur if you feed the FFT with a complex IQ signal that has non-zero imaginary components (for instance, baseband SSB).  Thus, if you want to see only a positive frequency peak, then you need to feed your FFT with a positive rotating IQ signal or phasor (cos(t) + i sin(t), IIRC).


Speech Signals - Recognition vs Analysis I am looking at working broadly in the area of Speech Signal Processing, and would be grateful if someone can clarify these doubts!

1) Are there any differences between Processing and Analysis of Speech Signals, in terms of the scope they provide for research? Arent both of them concerned in detecting the spoken language from signal?

2) What are the tools that one will be using in these areas?

Thanks!

1) Questions of terminology like that are always debatable! I would say that the term "analysis" describes techniques to convert the speech signal into another representation, closer to the semantic level than a raw list of samples (such as a written transcription in the case of speech recognition, a visualization of the tone curve in the case of a system assisting Mandarin learner, the identity of the speaker in case of speaker identification system, a category of emotion in case of an emotion detection system) ; while "processing" describes problems in which both the input and the output are speech (speech denoising/enhancement, coding and decoding, source separation, 'creative' effects like pitch/timing modifications). Of course, there are contexts in which "processing" describes the whole field, analysis included.

2) Deterministic signal processing (filtering, transforms), statistical signal processing (auto-regressive models, sinusoidal model estimation), statistical modeling and machine learning, psychoacoustics,  knowledge of how humans produce speech... If you venture into recognition: language theory (finite state machines), natural language processing, some elements of distributed computing (for large model training), data structures/algorithms. If you venture into speech coding or enhancement applications: information theory, experiment design (for subjective evaluation campaign), DSPs/embedded systems with hardware acceleration or even in some cases hardware design. Tools and language: mostly matlab or scipy for research, mostly C/C++ for implementation; and a scripting language (perl, python) is always helpful when you have to deal with texts/annotations.


How to design an inverse for a Butterworth low-pass filter? I have a signal that has a spectral shape very similar to the transfer function of a second-order lowpass filter (butterworth design).

I would like to whiten the spectrum of this spectrum in the time domain. For that, I need a filter that reverses the effect of the second-order lowpass filter.

I can calculate the inverse filter of the second-order lowpass by exchanging the A- and B-coefficients, but this will amplify high frequency components to infinity.

How do I design a filter that has the inverse transfer function of a second-order butterworth design but does not amplify to infinity at $f_s/2$? The transfer function should have the same corner frequency, but smoothly rise with 12 dB per octave.

Is this even possible? Any pointers as to how to design such a filter would be helpful!

The simplest and probably best approach is to use a second order IIR filter exchanging numerator and denominator from the original filter - as you suggested - but move the double pole at $z=-1$ (Nyquist) slightly towards $z=0$, i.e. slightly decrease the pole radius to make the filter stable. Just play around with the exact value. The higher the value, the closer to Nyquist will the approximation of the inverse filter be.

This little Matlab/Octave program shows how to do it:

[b,a] = butter(2,.005);     % some 2nd order Butterworth LP filter
r = 0.98;                   % desired pole radius < 1
b2 = a;                     % EQ numerator
a2 = [1,2*r,r*r];           % EQ denominator
k = sum(a2)/sum(b2);
b2 = k*b2;                  % scaling
[H1,w] = freqz(b,a,512);    % Butterworth frequency response
H2 = freqz(b2,a2,512);      % EQ frequency response

% plot EQ'd response
plot(w/pi,20*log10(abs(H1.*H2)))
grid on, xlabel('normalized frequency'),
ylabel('dB'), title('equalized response')


MATLAB Rayleigh fading and white gaussian noise From the MATLAB code below where do these theoretical equations for Rayleigh fading and white Gaussian noise come from? Or how are they derived?

h = 1/sqrt(2)*[randn(nrx,ntx,N/ntx) + 1i*randn(nrx,ntx,N/ntx)]; % Rayleigh channel
n = 1/sqrt(2)*[randn(nrx,N/ntx) + 1i*randn(nrx,N/ntx)]; % white gaussian noise, 0dB variance

The Rayleigh fading channel equation you provided comes from the property that given two independent zero-mean Gaussian random variables with equal variance $X \sim N(0,\sigma^2)$ and $Y \sim N(0,\sigma^2)$, the random variable $R = \sqrt{X^2+Y^2}$ is Rayleigh distributed (see for example wikipedia). In the code you provided, the real and imaginary components (generated by independent calls to randn) generally meet this condition (or at least approximates it quite well for reasonable pseudo-random generator) and the magnitude of h would thus have a Rayleigh distribution.

In addition, it is generally assumed that the signals' power is preserved on average, that is $E\{R^2\} = 1$. Now given how we defined $R$:
$$
\begin{align}
  E\{R^2\} &= E\{X^2+Y^2\} \\
           &= E\{X^2\}+E\{Y^2\} &\mbox{(linearity of expectation)} \\
           &= 2E\{X^2\} &\mbox{($X$ and $Y$ identically distributed)} \\
           &= 2\sigma^2 &\mbox{($X$ and $Y$ zero-mean Gaussian)}
\end{align}
$$

So the signal power is preserved on average when $\sigma = 1/\sqrt{2}$. This is also the scaling factor that must be used if starting with unit variance Gaussian pseudo-random variables (as is the case with MATLAB's randn).

Similarly, given a complex-valued Gaussian noise $n$ defined as $n = n_{\scriptsize \mbox{real}} + i\cdot n_{\scriptsize \mbox{imag}}$, where both real and imaginary components are Gaussian distributed with variance $\sigma_n^2$, the noise power is
$E\{|n|^2\} = E\{n_{\scriptsize \mbox{real}}^2 + n_{\scriptsize \mbox{imag}}^2\}$. By a similar argument as above, it follows that $E\{|n|^2\} = 2\sigma_n^2$. To obtain noise with unitary power, we thus need $\sigma_n = 1/\sqrt{2}$.

Note that the notion of unitary power only makes sense once a measurement unit has been established. Quite often the noise power is defined relative to the signal power. In that case, a noise power of 1 (or 0dB) in signal power's units would mean that the power of the noise is the same as that of the signal.


Echo generation of FIR system using Matlab How can I generate some echo to the system below?

Here is my code. I would appreciate any help.
% Template for the echo cancelling problem
fs = 8000; % sampling frequency (Hz)
f = 400; % sinusoid frequency (Hz)
Tdur = 0.2; % pulse duration (s)
x = pulse(f, Tdur, fs); % Generates a pulse
x = [x zeros(1, 2*length(x))]; % zero pad for processing
n = 0:length(x)-1; % discrete time vector
t = n/fs; % continuous time vector
plot(t, x); % Plot signal
xlabel('Time (s)')
ylabel('Amplitude')
title(sprintf('Original pulse: %.2f Hz', f))
sound(x, fs); % Play the sound
% Zero-pad original pulse
Nzp = floor(fs*Tzp);
nzp = 1:Nzp;
xzp = [x zeros(1,Nzp-length(x))];
figure(2)
plot(nzp/fs,xzp); xlabel('Time (s)'); ylabel('Amplitude');
title(['Original pulse, zero-padded to ' num2str(Tzp,3) ' s'])
pause; sound(xzp);
% Generate echoes
Tdel = 0.1; % echo delay (s)
alpha = 0.5; % echo gain
Ndel = floor(Tdel*fs); % echo delay (samples)
b= [1 zeros(1,Ndel-1) alpha]; a = 1; % filter coefficients
y = filter(b,a,xzp);
figure(3)
plot(nzp/fs,y); xlabel('Time (s)'); ylabel('Amplitude');
title(['Echoed signal, gain = ' num2str(alpha,3)]);
pause; sound(y);
%% Your code goes here
function x = pulse(f, Tdur, fs)
% Sinusoid multiplied by a Hann window function.
% f sinsusoid frequency (Hz)
% Tdur duration (s)
% fs sampling frequency (Hz)
nmax = floor(Tdur*fs); % duration (samples)
n = [1:nmax]; % discrete time index
x = hann(nmax).'.* cos(2*pi*f*n/fs); % waveform samples
end

Generally, an echo of a signal is not just its delayed and weakened copy, but is characterized by an impulse response of the echo path. The basic model of an echo is a linear one.
In order to simulate a linear echo $E$ one will implement the following convolution:
$E = S * H, $
where $S$ is your sampled signal and $H$ is an echo path impulse response, both with the same sampling frequency.
The particular form of the impulse response $H$ depends on the echo path. For example, $H$ may include $N$ zero values at the beginning to simulate an echo delay $N$.
Note, that in case of an acoustic echo, the linear impulse response $H$ may be not enough since an echo path is usually non-linear.


Does this equation have a closed form I am on chapter 6 of "Understanding digital signal processing", 3rd Edition by Richard G. Lyons. But I am stuck on question 6.15. Here is the network in question.

I wrote down the difference  equation as $$y(n)=x(n)+y(n-1)-Qy(n-1),$$ where $x(n)=D$(i.e. constant for all $n$), $Q$ is positive and less than $1$ and network is at rest at $n=0$ (i.e $y(n-1)=0$ when $n=0.$) In order to answer the question, I need to develop the series equation for $y(n)$ and find its closed form.

I have developed the series expression for $$y(n)=D\sum^n_{k=0}(-Q)^k {{(n+1)}!\over{(k+1)!(n-k)!}}.$$

And the terms without simplifying the factorials are $$y(n)=D{(n+1)!\over n!} - DQ{(n+1)!\over 2(n-1)!} +DQ^2{(n+1)!\over 6(n-2)!} -DQ^3{(n+1)!\over 24(n-3)!} . .. +(-Q)^nD$$.

My issue is that I cannot find where I have gone wrong because this series for $y(n)$ does not seem to have a common ratio which makes it hard for me to write it in the closed form as in  $$\sum^{N-1}_{n=p}r^n={r^p-r^N\over{1-r}}.$$

My question is, is my developed $y(n)$ suitable? If it is suitable, does it have a closed form?

Thanks

$$y[n] = x[n] + (1-Q) y[n-1]$$

where $y[-1] = 0$ and otherwise $n \ge 0$.

this is a difference equation with both a "natural response" (when the input is always zero, but you might have initial conditions) and a "forced response".  a.k.a. the "transient response" and the "steady-state response" the math guys might call these the "homogeneous solution" and the "particular solution".

this is LTI so then the exponential function is an eigenfunction for the natural response.

for the natural solution, try
$$ y_\text{n}[n] = y_\text{n}[0] a^n $$

then
$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 y_\text{n}[0] a^n &= 0 + (1-Q) y_\text{n}[0] a^{n-1} \\
               a^n &= (1-Q) a^{n-1} \\
                   &= (1-Q) a^n a^{-1} \\
               1   &= (1-Q) a^{-1} \\
\end{align}$$

so i guess $a = 1-Q$.

for steady state, you posit a solution that is like the input plus all of its finite differences.  since the input is constant and all finite differences of a discrete constant function are zero, then the steady state output is a constant.  you sum the scaled steady state with the scaled transient response for the complete response.

$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 C + y_\text{n}[0] a^n &= D + (1-Q)(C + y_\text{n}[0] a^{n-1}) \\
 C + y_\text{n}[0] (1-Q)^n &= D + (1-Q)(C + y_\text{n}[0] (1-Q)^{n-1}) \\
                           &= D + (1-Q)C + y_\text{n}[0] (1-Q)^n) \\
                      C &= D + (1-Q)C \\
\end{align}$$

i guess that means that $QC = D$ and $C = \tfrac{D}{Q}$.

now you add your natural response to your forced response.

$$ y[n] =  y_\text{n}[n] + C $$

lastly, for your initial condition

$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 y[0] &= x[0] + (1-Q) y[-1] \\
 C + y_\text{n}[0] (1-Q)^0 &= D + (1-Q)\cdot 0 \\
 C + y_\text{n}[0] &= D  \\
\end{align}$$

i guess that means that $ y_\text{n}[0] = D-C = D - \tfrac{D}{Q} $.

$$\begin{align}
  y[n] &= y_\text{n}[n] + C \\
       &= y_\text{n}[0]a^n + C \\
       &= (D - \tfrac{D}{Q})(1-Q)^n + \tfrac{D}{Q} \\
\end{align}$$

for $n \ge 0$.


concept of image rotation i am trying to understand how does the image rotation occur? i studied and came to know that we just multiply the image coordinates with rotation matrix which gives us the value of new coordinates. ok that is fine i get this. But main thing which is confusing me is the pixel values. i read that the pixel coordinates which we get after rotation can be non decimal hence interpolation is needed to get value at that point for this we rotate image back to original coordinate and find value by interpolation. But if we will rotate it back to original coordinate the the pixel will regain its original coordinate and at that coordinates value w are already having the value then why to do interpolation.please somebody clear my confusion.

When you multiply the image matrix with the rotation matrix, you are essentially multiplying an integer matrix with a real number matrix. The result is a real matrix, which has to be quantized back to a discrete matrix.

What is being quantized? The coordinates themselves. The integer coordinates of the previous position of the pixel is now a real valued coordinate, which must be converted back to an integer in order to be placed on the screen.

If you were asked to rotate it back to the original angle you could have been right in questioning the futility of interpolation. However, the user might not want to rotate it back, so you have to perform interpolation or quantization to be able to use the rotated image. If you are concerned about image fidelity, you might want to save a copy of the original real valued matrix in memory or file until or opertions are finished and the user exits program.

Here is a useful link you might want to peruse.