Find Discrete Fourier transform given the inverse I don't know if this question is already answered somewhere else but I'm stuck with that and I need help...

Given an inverse Discrete Fourier Transform (IDFT), with N=16 :

How can I find the F(k) that represent the Discrete Fourier Transform (DTF)?
$$
F(k) = \sum_{n=0}^{N-1} f(n) e^{-i2\pi \frac{kn}{N}}
$$

And how the plot of F(k) will be?

EDIT: for now this is my result but I have no idea if it's correct:

EDIT 2: plus I don't know how to plot the equation since it contains the imaginary number (j in the figure)...

METHOD-I

Plain simple,

The discrete-time sequence you have provided is :

$$
f[n] = \begin{cases} 1 ~~~,~~~\text{ for } n = 0,1,2,14,15 \\ 0 ~~~,~~~ \text{ otherwise } \end{cases}
$$

which is also described as a sum of discrete-time unit impulses located for each $n=0,1,2,14,15$ as

$$f[n] = \delta[n] + \delta[n-1] + \delta[n-2] + \delta[n-14] + \delta[n-15] $$

Using the fact that, the N-point DFT of the unit-impulse located at $n=d$ is: $$ \boxed{ \delta[n-d] \longleftrightarrow X[k] = e^{-j \frac{2\pi}{N} k d} }$$

Then employ linearity property of DFT, to get the N-point DFT of $f[n]$ as:

$$
F[k] = e^{-j \frac{2\pi}{N} k 0} + e^{-j \frac{2\pi}{N} k 1} + e^{-j \frac{2\pi}{N} k 2} + e^{-j \frac{2\pi}{N} k 14} + e^{-j \frac{2\pi}{N} k 15}
$$

which becomes with $N = 16$:
$$
F[k] = 1 + e^{-j \frac{\pi}{8} k} + e^{-j \frac{\pi}{4} k} + e^{-j \frac{\pi}{8} 14 k} + e^{-j \frac{\pi}{8} 15 k}
$$

Now, not always possible but in this case you can merge those complex exponentials, after replacing the frequency of the the latter two according to the following.

$$
e^{-j \frac{\pi}{8} 14 k} = e^{-j \frac{\pi}{8} (16-2) k} = e^{-j 2\pi k} e^{j \frac{\pi}{4} k} = e^{j \frac{\pi}{4} k}
$$

and similarly
$$
e^{-j \frac{\pi}{8} 15 k} = e^{-j \frac{\pi}{8} (16-1) k} = e^{-j 2\pi k} e^{j \frac{\pi}{8} k} = e^{j \frac{\pi}{8} k}
$$

Then the DFT $F[k]$ becomes:
$$
F[k] = 1 + e^{-j \frac{\pi}{8} k} + e^{-j \frac{\pi}{4} k} +e^{j \frac{\pi}{8} k} + e^{j \frac{\pi}{4} k}
$$

Now combine those exponentials, with conjugate angles, into cosine terms according to $2 \cos(\omega n) = e^{j \omega n} + e^{-j \omega n}$ to conlcude:
$$
\boxed { F[k] = 1 + 2 \cos( \frac{\pi}{4} k) + 2 \cos( \frac{\pi}{8} k) ~~~,~~~\text{ for } k = 0,1,2,...,15}
$$

Which is a real and (circularly) even function of $k$, as the sequence $f[n]$ is (circularly) even-symmetric and real.

METHOD-II

if you prefer a fancier approach, then you can also use this. By inspection, one can see that $f[n]$ is given by a circular left shift of another sequence $g[n]$ (i.e.,  $f[n] = g[n+d]$) where and $g[n]$ is

$$
g[n] = \begin{cases} 1 ~~~,~~~\text{ for } n = 0,1,2,3,4 \\ 0 ~~~,~~~ \text{ otherwise } \end{cases}
$$

and $d=2$. Then the N-point DFT of $f[n]$ and N-point DFT of $g[n]$ are related by (with $N=16$):

$$
F[k] = e^{j \frac{2\pi}{N} d k} G[k] ~~~,~~~\text{ for } k = 0,1,2,...,N-1 \\
F[k] = e^{j \frac{\pi}{4} k} G[k] ~~~,~~~\text{ for } k = 0,1,2,...,15 \\
$$

Ffurthermore, the N-point DFT of $g[n]$ is found as:
$$
G[k] = \sum_{n=0}^{4} 1 \cdot e^{-j \frac{2\pi}{16}nk} = \frac{ 1 - e^{-j \frac{5 \pi}{8} k} }{ 1 - e^{-j \frac{\pi}{8} k} }
$$

Finally the N-point DFT $F[k]$ of $f[n]$ becomes
$$
\boxed{
F[k] =e^{j \frac{\pi}{4} k} \left( \frac{ 1 - e^{-j \frac{5 \pi}{8} k} }{ 1 - e^{-j \frac{\pi}{8} k} }\right) ~~~,~~~\text{ for } k = 0,1,2,...,15
}
$$

Which should yield the same result.


Matlab fast fourier transform (FFT) function and time in hours frequency : hertz or millihertz? I have a function of type

$$C(t)=0.31 \left(e^{-0.267t}-e^{-1.72t}\right)$$

Where $t$ is in hours

The function represents the evolution in a window of 24 hours of a drug blood concentration, after a dose of a drug is given at  $t = 0$ to a patient. I would like to use Fourier transform, to identify th best sampling step to reconstruct the concentration curve. Because it is a human whose blood we are sampling here I would like to reduce the sampling as much as possible, with a reconstruction error of around 20%.

I am using Matlab FFT function to compute the function's Fourier transform and plot the spectrum.

However I am not sure if my frequency are in Hertz or Millihertz. Because I am sampling in hours, I am expecting lower frequencies (millihertz), but I am not sure of the logic.

Here is my Matlab Code

clear all; clc;close all;
%% Time-domain signal
Ts=0.25; % time discretization step
t=0:Ts:24 % in hours
C=0.31* (exp(-0.267*(t))-exp(-1.72*(t)));

figure,plot(t,C);

title('Time domain signal');
xlabel('Time in hours ')

%Frequency domain signal Via Matlab FFT function
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz or Millihertz?')

That function is not periodic. Using the Fourier transform to analyze its behaviour is questionable, to say the least. Maybe you're confusing things with the Laplace transform, which is more useful when analyzing control systems?

Why should the axis in your plot have Hertz as unit?

The discrete Fourier transform¹ does a base transform from a vector in $C^N$ to a vector in $C^N$; its physical significance is that the result is that each bin represents another frequency between $0$ and the the Nyquist rate (which, in the complex sampling case is the sample rate, and in the real sampling case is half the sample rate).

Since you only have a real function, I'd argue frequencies are relative to $\frac{f_{sample}}2$.

$f_{sample}$ is $\frac{4}{\text{hour}}=\frac4{60\cdot60\,\text{s}}= \frac1{900}\,\text{Hz}$ in your case. So each bin represents a frequency $$f\in [-\frac{f_{sample}}2 ; \frac{f_{sample}}2[\,=\,[-\frac{1}{1800}\,\text{Hz}; \frac{1}{1800}\,\text{Hz}[$$.

Now, frequency base vectors are uniformly spread over the nyquist rate, so the $n^\text{th}$ bin represents frequency $\frac{\frac{1}{900}\,\text{Hz}}{n_\text{bins}}$. $n_\text{bins}$ is the number of bins, which is the same as the input vector size; which is $4\cdot24=96$. So each bin is $\frac1{96\cdot900}\,\text{Hz}= 11.5741\,\mu\text{Hz}$. Not very handy.

I think you're somewhat confused what an FFT is. I think a textbook on signal theory will do you very well!

¹ you're using fft, which is a fast implementation of the DFT; not the continuous Fourier transform. That's an important difference.


Differences between OpenCV Canny and MatLab Canny? does anyone know why the MatLab Canny (MLC) is so different compared to the OpenCV Canny (OCC)?
ML-C delivers precise and more connected edges than the OCC, but how is that possible?
The reason why I ask is, that I need to implement my protoype of ML code into C++ and I wanted to use OpenCV.
Exporting code of ML isn't really possible as far as I have tried.

Kind regards,

As suggested above, the Matlab Canny edge detector calculates the gradient using a "derivative of a Gaussian filter" (as stated in the documentation). In other words, Matlab does a Gaussian blur of the image and then finds the gradient of that smoothed image... all using a single fancy filter. [If you want to know the details, just type in edit edge as Andrey suggested, and then scroll down to the smoothGradient() function.]

The blurring operation significantly reduces the amount of noise present in the image, eliminating many spurious edges and leaving behind the good stuff.

Unfortunately, the OpenCV Canny function doesn't let you change the filter kernel it uses via the function parameters. However. You can generate the same results by first blurring the input image, and then passing this blurred image into the Canny function.

This significantly cleans up the resulting edge map. To blur the input image, I personally use OpenCV's GaussianBlur() function with sigmaX=2. This mimics the default sigma in Matlab. The best blurring kernel size can vary from case to case, but in Matlab it is calculated using filterLength = 8*ceil(sigma);, so for a sigma of 2 that would mean a kernel size of (16,16)

Since both the Gaussian Blur and Sobel filters are linear, passing a blurred input image to the OpenCV Canny() function is mathematically equivalent to what Matlab does because of the principle of superposition, as demonstrated in this pseudocode (note: * is the convolution operator):

// The Matlab method: the sobel and blur operations are combined into
// a single filter, and that filter is then convolved with the image
matlabFancyFilter = (sobel * blur);
gradient = matlabFancyFilter * image;

// Equivalent method: image is first convolved with the blur filter, and
// then convolved with the sobel filter.
gradient = sobel * (blur * image); // image is filtered twice

There is an OpenCV Canny tutorial here that demonstrates how to do this using C++. I'm a python guy, so here is what I do:

smoothedInput = cv2.GaussianBlur(image, (7,7), 2);
edges = cv2.Canny(smoothedInput, 25, 50);


DFT/FFT Transfer function I want play and record a sine sweep.
When i have both signals the recorded one and the send one i can create a Transferfunction.
That is what i know so far.

$$
H_0 = \frac{OUT}{IN} = \frac{Y}{X}
$$

Where i'm stuck is that when i read about the Transfer function it is all about the $Laplace \, transform$ (and sometimes the $Z-Transform$)
While i use the $FFT/DFT$. I tried some GNU Octave code to tried an Transfer function and i think that looks what i suspected.
But now i curious

if the $FFT/DFT$ can be used for/ as a transferfunction and
what is the relation between an $FFT/DFT$, $Laplace \, Transform$ and the $Z-Transform$

Your question is fairly broad, let me answer it step by step.

First of all $H(s)$ is indeed called the transfer function, and is the laplace transform of the impulse response. It's useful for finding the poles and zeroes, which is what the fourier transform can't do alone. $H(\omega)$ is the frequency response of the impulse response, not really the "transfer function".

If the frequency response is what you want to attain, then you can use the frequency response (fourier transform). It does not exist for unstable signals though.

The relation between Fourier Transform, discrete Fourier transform and Laplace transform/z-transform:

Fourier transform and DFT are in principle the same. One is for discrete signals, and the other for continuous signals. Because of that there are some different properties like frequency domain aliasing. What both of these transforms do is decompose the signal into complex exponentials (easily converted into sines and cosines (or sinusoids with magnitude and phase), which is the original idea).

Laplace and z-transform are again in principle the same, z-transform is the discrete equivalent to Laplace. It's simply the Fourier transform/DFT with the signal multiplied by a varying exponential $e^{\sigma t}$, or in discrete case just $r^t$. Then you have multiple fourier transforms/DFTs relating to each different value of the exponential ($r$ or $\sigma$). Point of interests are the poles and zeroes, IE. the points where the impulse response in question  multiplied by an exponential has integral that is exactly infinite or exactly zero.

Note that when $r=1$, or $\sigma = 0$, the Laplace/Z-transform is equal to the Fourier transform. In other words, the Laplace/Z-transform contains the Fourier transform in it.

EDIT: I wrote this fairly quickly so please point out any errors.


What should my reference value be when converting FFT bin amplitudes to dB? I want to transform my FFT output values into a dB scale, but I'm struggling to determine the function I should run each bin amplitude through.  My understanding of the decibel scale is that a value representing an intensity must be divided by some reference intensity.  The log of this ratio, times 20, will represent the dB scale value:

$$
\text{dB} = 20 \log_{10}\frac{\text{intensity}}{\text{intensity}_{\text{ref}}}
$$

I should be able to pass each an FFT bin amplitude value into this equation as $\text{intensity}$, and the result will be the dB level of that bin.

Unfortunately, I'm not sure what $\text{intensity}_{\text{ref}}$ should be.  Based on my reading, I believe it should represent the maximum possible intensity.  If the readings were samples from a wave, this would be easy (the highest possible reading for 8-bit audio would be 255; for 16-bit it would be 65535, etc).  However, I'm not sure that FFT outputs have such a predictable range or scale.

One option would be normalizing the amplitude values (dividing each value by the maximum amongst all bins for all FFT windows computed for the signal).  This way, the bin with the highest amplitude would be at 0 dB and all others would be negative values relative to that maximum.  However, this would only work if I pre-computed every FFT window for an entire signal.  What if I was trying to compute dB values for an FFT of a just a single window?

So my question is: what is $\text{intensity}_{\text{ref}}$?  I'm still relatively new to signal processing, so if there's an easy way to guarantee some range of FFT amplitude values or any other obvious solution I'm missing please let me know.

You tagged this dbfs, so, that's basically your answer, because, yes, the output of a DFT is very much predictable!

What you choose, however, as reference amplitude is up to you to define as useful to you.  I can but make a recommendation:

Assume you're transforming a constant vector of $N$ maximum (ADC) values; that'll give you, depending on the normalization that your DFT applies, an all-0 vector with the equivalent energy in the DC-bin. That'd be a useful value for the full-scale maximum; frequency shifting and superposition with other signals can never yield a more powerful bin.


FFT of any sinusoidal signal -- how many point fft to take? suppose I have to measure the DFT of an 60 hz sinusoidal signal and calculate upto 20 harmonics of this.  Value of the 20th harmonic will be 20*60 = 1200 HZ. Then my sampling frequency should be >= 20*60*2 >= 2400 hz

Now how many point FFT i have to calculate to estimate
information about harmonics ?
Is number of points in FFT(N) equal to the sampling frequency ?
How number of points in FFT(N) related to the sampling frequency
?
Is there any advantage with output off FFT if we sample at
higher rate ?

Is bellow three estimates are correct ?

1> If suppose my sampling frequency is fs =  20*60*2 = 2400 then, What i know is that the frequency range [0,fs] is represented by 2400 point.
   I.e each point represent 2400hz/2400 = 1 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

   X[0] is the constant term
   X[1] is 1 hz
   X[2] is 2 hz
   X[1201] is the 1200 HZ

2> If suppose my sampling frequency is fs =  20*60*2 = 2400 then, What i know is that the frequency range [0,fs] is represented by 2400 point. can  i take 4800 sample points ?
   I.e each point represent 2400hz/4800 = 0.5 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

   X[0] is the constant term
   X[1] is 0.5 hz
   X[2] is 1 hz
   X[2400] is the 1199.5 HZ
   X[2401] is the 1200 HZ

Is there any advantage here to take more number of samples to calculate the FFT ? One thing i note is that spectrum is more dense with each point having 0.5 HZ of information, will my calculation of FFT will be more accurate in this case ?

3> If suppose my sampling frequency is fs =  20*60*4 = 4800  =  then, What i know is that the frequency range [0,fs] is represented by 4800 point.
   I.e each point represent 4800hz/2400 = 1 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

X[0] is the constant term
   X[1] is 1 hz
   X[2] is 2 hz
   X[2401] is the 2400 HZ

Now in this case FFT is generating information upto 2400 hz but i need upto 1200 hz. Is there any advantage with output off FFT if we sample at higher rate ?

Please suggest on this.

Many of your questions are answered elsewhere, so I added additional links below but to summarize in one place relative to your questions I offer the explanations below that together with the links should answer your questions:

Required Sampling Frequency:

If the signal is real, then yes the sampling frequency must be greater than twice the highest frequency for a signal bandwidth that extends down to DC. However if the signal is complex then the sampling frequency should be greater than the highest frequency for same case. This is because a complex signal has a unique spectrum over the Nyquist frequency range of $\pm f_s/2$ where $f_s$ is the sampling rate, while for a real signal the negative frequencies are the complex conjugate of the positive frequencies (identical magnitude, opposite phase).

Note that prior to sampling an analog anti-alias filter is required on the signal, otherwise energy if present in higher Nyquist zones (frequencies above the Nyquist frequency range defined above) will alias into the Nyquist frequency range being observed.  (see this post which explains how aliasing occurs: Aliasing after downsampling). It is not feasible to implement a brick-wall filter, so additional margin to the sampling rate is usually added to allow for the characteristics of the anti-alias filter. For this you need to understand the properties of the signal signal of interest prior to sampling (as once sampled there is no way to distinguish an alias from a signal in the first Nyquist zone!). For example in your case, where you are interested in observing the signal out to 20 Harmonics of 60 Hz: The anti-alias filter should pass the 20th harmonics without distortion and then reject the higher harmonics so that you can minimize the sampling rate needed. You would then choose the sampling rate depending on your filter implementation such that no energy folds in to distort the 20th harmonic.

How Many Point FFT is Needed:

This has been answered several times, so see the links given below for further details and caveats, but the main point is that the frequency resolution of each "point" is 1/T Hz where T is the length of your signal in seconds.  Further the bins in the DFT (the FFT is just an algorithm that computes the DFT) as typically returned extend from DC to one bin less that $f_s$. Due to aliasing the DC is equivalent to $n f_s$ where n is any integer, and likewise for all other frequencies in the DFT, so the DFT that extends frmo DC to $f_s$ equivalently extends from DC to one bin less than $f_s/2$ and then from $-f_s/2$ to one bin less than $f_s$ (see fftshift in Matlab). As detailed in the links given, due to the relationship of the sampling rate and the total length in time T of your signal, each bin will have a frequency resolution of $1/f_s$ when no windowing is used.

Note importantly that each bin of the DFT is the result of a Sinc filter, so if your sampling rate is not an integer multiple of 60 Hz, you will also have the effect of "scalloping loss" that will distort the measurement of each harmonic. (see the "Bank of Filters" explanation at this link: Intuition for sidelobes in FFT). By choosing a sampling rate that is an integer multiple of 60 Hz, each 60 Hz harmonic will be properly centered in each bin.

Suggestion for your specific application:

With relationship to choosing the sampling rate and number of bins for your specific application, if the 60 Hz and its harmonics are the dominant signal and you don't have other signal energy of significance then there is a simple solution that will work quite elegantly in accurately estimating each harmonic with minimum measurement distortion. As described above, choose a sampling rate that is an integer multiple of 60 Hz, high enough to allow for the anti-alias filter used, and then choose the number of points such that each bin passes a harmonic of 60 Hz and rejects all the other harmonics. This means the number of samples in the FFT is $f_s/60$ where $f_s$ is an integer multiple of 60 Hz > 2400 Hz (assuming your signal is real). For example 4800 Hz may be a good choice as you suggest in your options (allowing for a fairly relaxed anti-alias filter!), and in this case $N =4800/60= 80$ points.  So in this case $f_s=4800 Hz$ and with 80 points each point is 60 Hz, so 0, 60, 120, 180 etc centered on each harmonic and perfectly nulling the other Harmonics. The actual magnitude measured for each point is the total area under the Sinc function frequency response for each bin, so if there is other energy present at non-harmonic frequencies, this would distort the magnitude measured in each bin. But if the harmonics are dominant, this would be an excellent and simple approach.  This is somewhat similar to the objective in this link in removing harmonic signals (Different way to separate a particular frequency from a signal) as each bin in the DFT is a rotated moving average filter just as done in this referenced link when only one bin is of interest. Note from that link how every other harmonic is nulled!

Please see these other responses which will further help give insight:

Specific frequency resolution

FFT frequency resolution

Number of FFT points required for a specific frequency resolution for an oversampled signal

FFT and number of samples relations

What happens when N increases in N-point DFT


phase offset extraction using matched filter I'm looking for an application of detecting the drift of a received signal carrier frequency from ADC clock of the receiver (in case of no drift, ADC clock should be optimally an integer multiple of signal's carrier frequency), which can be observed as a varying offset component of the instantaneous phase. one proposed solution that I found was to use a matched filter to extract the variable offset. Yet I didn't find any supporting description to understand the mechanism. did anyone come across such an implementation?

This is a common implementation for carrier recovery and tracking.

A matched performs a correlation of the received sequence with the expected sequence as given by:

$$\sum_{n=0}^{N-1}r[n]x^*[n]$$

Which states a correlation over N samples is performed by multiplying and accumulating the time-aligned received sequence $r[n]$ with the (complex conjugate of for complex signals) the expected ideal sequence $x[n]$. Note how any phase shift on $r[n]$ would result in the same phase shift on the result (any phase introduced on the ADC sampling clock would be translated to $r[n]$ directly as the sampled received signal):

$$\sum_{n=0}^{N-1}r[n]e^{j\theta}x^*[n] = e^{j\theta}\sum_{n=0}^{N-1}r[n]x^*[n]$$

Often a frequency offset is involved (which is a phase linearly changing with time), in which case we need to be aware of the resolution bandwidth of the matched filter, which is how the magnitude of the correlation result reduces as a function of frequency offset. For example If the matched filter is on constant magnitude coefficients (such as a moving average filter), the correlation magnitude will go down as a Sinc function with the first null at $1/T$ Hz where $T$ is the time length in seconds of the correlation.

But given the correlation result provides the carrier phase as shown above, if there is a small frequency offset, we can measure that precisely with successive results of the matched filter output given we will be observing a change of phase with a change in time (which is frequency given by $d\phi/dt$).

This is demonstrated below with captured data for live GPS signals where the matched filter (in this case for satellite SV24) is providing a sliding correlation over the 1 ms long GPS course-acquisition code, resulting in the spikes every 1 ms in the plots below indicating when the code aligns with the received signal. However the correlation is going up and down from to the underlying phase rotation due to a frequency offset.

This rotation is clearer on an IQ plot of the same result above, plotting the magnitude and phase on a complex polar plot:

By evaluating the angle over time of the correlation peaks, the frequency offset was determined in this case to be 38 Hz, and after correction the same above plots appear as follows, revealing the GPS data message (we see from the IQ plot that a small residual phase error exists that could also be removed):


How can I calculate what a recording would sound like if it was played through my speakers and mic? I am trying to isolate the sound of someone speaking while music plays back in a noisy environment.

I have multiple speakers and a microphone in my system. I am able to

play back just music and record it, then
play back music while speaking and record it.

After that, I can look at the difference in the signals and get the meaningful information I need.

In order for this to work in real-time, I would need to calculate what that first recording would sound like. If I have the impulse+magnitude response or some other data about the system, is it possible to predict what the microphone will detect?

EDIT:I suspect that some kind of transform could be derived from the impulse response that will give a decent prediction of what the mic will hear. To be clear, this doesn't need to be incredibly accurate. I just want a general idea of what the mic will pick up.

You can use convolution to hear a recording as it would sound through a specific speaker. The linked page describes how to get an impulse response and convolve a signal with it.

https://www.mathworks.com/help/audio/examples/measure-impulse-response-of-an-audio-system.html


Capacity loss of fast fading regarding AWGN The theoric capacity of fast fading channel is:
$$
C_{\rm fading}=E\left\{\log\left(1+\lvert h \rvert^2\textrm{SNR}\right)\right\}, \quad\text{with $E$ is expectance operator.}$$

At high SNR:
$$\begin{align}
C_{\rm fading}&=E\left\{\log\left(\lvert h \rvert^2\textrm{SNR}\right)\right\}\\
&=\log(\textrm{SNR})&&+E\left\{\log\left(\lvert h \rvert^2\right)\right\}\\
&=C_{\rm AWGN} &&+ E\left\{\log\left(\lvert h \rvert^2\right)\right\}
\end{align}$$

The difference regarding AWGN is $\Delta=E\left\{\log\left(\lvert h \rvert^2\right)\right\}$ which is stated as $-0.83\textrm{ bits/s/Hz}$ for standard Rayleigh fading model, or $2.5\textrm{ dB}$ requirement.

Could anyone explain how they come to this result $-0.83\textrm{ bits}, 2.5\textrm{ dB}$ ?
Standard Rayleigh fading is $C\mathcal N\left(0,1\right)$ then $E\{\lvert h \rvert^2\} = 1$. Using Jensen's inequality $$\begin{align}\Delta&=E\left\{\log\left(\lvert h \rvert^2\right)\right\}\\ &\ge \log(E\{\lvert h \rvert^2\}) \\&= \log(1) \\&= 0 \\&> -0.83\end{align}$$ and this $\Delta$ should not be $-0.83$? I think I made some mistakes but I don't know what it is.

The (stupid) mistake I have committed is that log is concave, so the inquality is in wrong side.
A numerical evaluation of the integral proposed by Marcus confirmed the value 0.83 bits.
Thank Marcus.

Standard Rayleigh fading is $C\mathcal N\left(0,1\right)$ then $E\{\lvert h \rvert^2\} = 1$. Using Jensen's inequality $$\begin{align}\Delta&=E\left\{\log\left(\lvert h \rvert^2\right)\right\}\\ &\ge \log(E\{\lvert h \rvert^2\}) \\&= \log(1) \\&= 0 \\&> -0.83\end{align}$$ and this $\Delta$ should not be $-0.83$? I think I made some mistakes but I don't know what it is.

Rayleigh fading is when the amplitudes of the real and imaginary part of $h$ are normal distributed; so $|h|$ follows a Rayleigh distribution, which has a probability density function  of

$$
p(x=|h|) = \frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\text.
$$

Hence, $E\left\{\lvert h\rvert^2\right\}$ is not the variance of the Normal distribution with variance $\sigma^2$.

Instead, the value you're looking for is
$$\begin{align}
E\left\{\log_2\left(\lvert h\rvert^2\right)\right\} &= \int_{-\infty}^\infty \log_2\left(\lvert x\rvert^2\right)\frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\,dx\\
&= \int_0^\infty 2\log_2\left(\lvert x\rvert\right)\frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\,dx\\
&= \frac2{\sigma^2}\int_0^\infty x\log_2\left( x\right) e^{\frac{x^2}{2\sigma^2}}\,dx
\end{align}$$


Why must the Region Of Convergence (ROC) contain infinity and the system function be a right-sided sequence for it to be causal? I need the simple logic on the condition a system becomes causal. We know that causal contains only past values. I can't relate this with the Region Of Convergence (ROC) concept.

A causal impulse response is zero for negative argument:

$$h[n]=0,\qquad n<0\tag{1}$$

Hence its $\mathcal{Z}$-transform is given by

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=0}^{\infty}h[n]z^{-n}=h[0]+h[1]z^{-1}+h[2]z^{-2}+\ldots\tag{2}$$

Note that there are no positive powers of $z$ in Eq. $(2)$. Consequently, $H(z)$ converges for $|z|\to\infty$, which means that infinity is inside the region of convergence.

Specifically you have

$$\lim_{z\to\infty}H(z)=h[0]\tag{3}$$

Eq. $(3)$ is called the initial value theorem of the $\mathcal{Z}$-transform.


How to normalize predesigned IIR filter coefficients to have unity gain? There is a low pass filter with difference equation already provided.

$$
y_n = (x_n - 2x_{n-6} + x_{n-12}) - (-2y_{n-1} + y_{n-2})
$$

$$
a_0 = 1, a_1 = -2, a_2 = 1
$$
$$
b_0 = 1, b_6 = -2, b_{12} = 1
$$

However the filter was designed with some gain as the filtered output is significantly exceeds original signal dynamic range x1 = scipy.signal.lfilter(b, a, x)

Is it possible to normalize filter coefficients to have the same as input signal dynamic range or approximately redesign it from filter frequency response?

In general, divide each of the $b$ coefficients by the maximum gain of your filter. This is equivalent to adding an attenuator before or after the original filter.

In this particular case, the maximum gain seems to occur at frequency 0 (where there is a double zero-pole cancellation).

At frequency 0, the gain of your filter is 36. So divide $b[]$ by 36.

Edit: the math

The frequency response comes from evaluating $H(z)$ at $z = e^{j\theta}$. Substituting and operating, and ignoring the phase, you get:

$|H(e^{j\theta})| = \frac{\sin(\theta/2)^2}{\sin(6\theta/2)^2}$

Taking the limit as $\theta \rightarrow 0$ you get $H(e^{j0}) = 36$.

Note that this is a simpler (and well known) filter applied twice. It corresponds the moving average of 6 samples (with the 1/6 factor missing).

So it is like two filters in cascade, each of which sums 6 consecutive samples. Each is a low-pass filter with gain 6 at 0 frequency.


Why did a digital HPF increase the gain of high-frequency signal? Background: I was trying to remove the DC component of a load cell signal using a high-pass filter. I don't care about that anymore. I'm curious why a digital filter and an analog filter behave differently under the same conditions.

Here (link 1, link 2) is the data I'm working with, and below is my Python code:

import numpy as np
import math as m
from scipy import signal
from scipy import fftpack
from myLib import myFFT
import matplotlib.pyplot as plt

# load data

filename = "raw_data.txt"
fSampling = 1E3
fNyquist = 0.5*fSampling
data = np.loadtxt(filename,skiprows=0)
i1 = 1000
i2 = 5000
t = data[:,0]
x = data[:,3]

# remove DC component

fCutoff = 0.01/fNyquist
bd, ad = signal.butter(4, fCutoff, 'high', analog=False)
xd = signal.filtfilt(bd,ad,x)
ba, aa = signal.butter(4, fCutoff, 'high', analog=True)
xa = signal.filtfilt(ba,aa,x)

# plot results

f, axarr = plt.subplots(3, sharex=True)
axarr[0].plot(t, x)
axarr[0].set_title("raw")
axarr[0].set_ylabel("x [-]")
axarr[1].plot(t, xd)
axarr[1].set_title("digital")
axarr[1].set_ylabel("x [-]")
axarr[2].plot(t, xa)
axarr[2].set_title("analog")
axarr[2].set_xlabel("t [s]")
axarr[2].set_ylabel("x [-]")

plt.show()

Problem The digital filter doubles the values, while the analog filter does not.

Why does the digital filter behave differently? Is something weird about my data, or do digital and analog filters behave differently? I'm new to signals processing (this all originated from me leaving the analog flag with the default "False"), but I thought analog and digital filters were generally different means to the same end.

Edit: I did some research online, and I've learned that I should use a digital filter to filter regularly sampled data. I don't know why yet (digital filters perform better at a known range?), but this is just making even more curious why my digital filter increased the gain of my signal. I analyzed the frequency response of the digital filter, and theoretically I should get a gain of 0:

I want to chalk this up to numerical errors since the digital filter's numerator and denominator are so close, but so are the analog filter's numerator and denominator:

b digital = [ 0.99991791 -3.99967164  5.99950746 -3.99967164  0.99991791]
a digital = [ 1.         -3.99983581  5.99950745 -3.99950746  0.99983583]
b analog = [ 1.  0.  0.  0.  0.]
a analog = [  1.00000000e+00   5.22625186e-05   1.36568542e-09   2.09050074e-14
   1.60000000e-19]

The following coefficients $a$ and $b$ implement a second order digital DC notch filter :

b = 0.9900   -1.9800    0.9900

a = 1.0000   -1.9800    0.9801

Whose frequency response is (y-axis in dB):

Given the raw-input signal $x[n]$ as a vector of samples in $x$ in matlab, you would apply the following line to filter it:

y = filter(b,a,x);

Note that the filter is not linear phase, hence you will have some phase distorions especially about DC frequencies... if you need a strictly linear phase DC blocking filter please indicate.

Given the raw data with the following plot:

The resulting filtered signal plot is the following:

It can be seen that initially there is a transient starting from a value of $0.428$ and takes some samples to stabilize to remove the DC. This transient is a direct consequence of long impulse response associated with the IIR filter.

Shorter impulse responses would provide faster response, but more steady state distortions.

Assuming that a shorter filter would do better, here is another set of DC notch coefficients with shorter impulse response:

b = 0.9100   -1.8200    0.9100
a = 1.0000   -1.8200    0.8281

In which case the output becomes:

Clearly an improvement in dynamic response, but keep in mind that the steadystate response now gives more distortions.

You can even further get a quicker response with the following set:

b = 0.7000   -1.4000    0.7000
a = 1.0000   -1.4000    0.4900

whose result is:

Much faster response, but there could be untolerable distortions on the samples, depending on the application requirements.


edge detection evaluation medical test images MATLAB I programmed an algorithm for edge detection, i want to compare it with other algorithms but with a criteria (Pratt's figure of merit, GMC1) not only a a visual comparaison.

the thing is PFOM needs an image reference (known ground truth) which i don't have, and i tried to find a GMC1 algorithm but i failed.

if any body can help me with reference images (known ground truth) (xray) or implementation would be much appreciated

The closest resource I have come across is the BSD500 database:

http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html

It is a set of natural images that have had the contours drawn by humans. When you download the test images there are also MATLAB files with all the contours.

However, the question is always "what is an edge?" Is it a step in brightness in the local image patch? Is a line an edge? Is it a boundary between semantic objects in a 2D image? The latter is what the above set tests. The best performing approaches for it also consider an abrupt change in texture to be an edge.

I guess the point is the performance of your edge detector is dependent on the application. One could argue that every edge detector is perfect as it perfectly detects edges according to its own definition. But that is just being cheeky.

You could also test some other properties, such as

Is it invariant to rotation? Does the same line at 0 degrees and X degrees give the same result?
Is it invariant to scaling? I.e. zooming of the image
Is it invariant to view point? I.e. warping of the image
Is it illumination invariant? Does it detect true weak edges as well as strong edges?

I also have a new edge detector that I've been working on. One of recent ideas I have had is to use the Mikolajczyk test set. It consists of eight sets of images with various transformations, including the above. Typically the test set is used to evaluate interest point descriptors.

The process would be

Detect edges in each image in the set.
Using the homographies, transform each detection image back to the first image in the set so they can be compared
Compare the edges between images using one of the measures you mention.

Ideally, if the edge detector is invariant to the various translations, the same edges should be found. This approach might have already been done in the literature but I haven't looked.


Decompose an image into better compressible components Our team develops an application that, in some feature, connects to a scanner (via TWAIN or WIA) and scans a document, and makes a PDF file out of it. The application is mostly used to process mail or documents with few images.

We recently noticed that the PDF files generated by our application are four times larger (speaking of file size) than those generated by the scanner software, although similar scanner parameters are used (resolution, color, etc).

So I delved into the generated pdf files and this is what I got:

With MuPDF component mutool.exe I got information about the objects inside the PDFs:

Our PDF contains one image (not surprising, it is made directly from the result of the scan):

[ DCT ]      3302x4668 8bpc DevRGB

The PDF file generated by the scanner software contains two images (aha!, this is more surprising), that differ in size and, that do not have the same color profile:

[ DCT ]      1650x2334 8bpc DevRGB
[ CCITTFax ] 3176x4512 1bpc ImageMask

Still using mutool, I extracted the two images from the PDF generated by the scanner software:

The first image is what we could call the “Background”: here are a thumbnail of it, and a detail at 100% size:

The second image is a pure B&W (1 bit per pixel) layer, with mostly black; here is a detail at 100% size:

So it appears that the software that the scanner comes with uses an algorithm that detects that the document is mostly black and white, and then decomposes the image in two parts that it overlays; the first part consists of small details that can be downsampled (and therefore take less space), and the second part consists of only black or white dots, and with which better compression can be achieved.

My question is: how can we do something similar?

Thanks

Assuming that the image is gray level (not color)

Use Otsu's thresholding method to detect whether the image has bi-modal histogram.
If Otsu tells you with high confidence that the image is bi-modal, use the threshold that it gave you in order to binarize the image. That is your binary mask that can be effectively packed into bits.
Find the gray values of the two histogram modes, create a mask where the "white" color is the gray value of the second mode and "black" color is the value of the first mode.
Subtract previously made mask from the image. These are your small changes that you can compress by using lossy compression.
Save the two modes as two integers together with the document.


Z domain transfer function to difference equation I want to convert this transfer function: $$\ \frac{2\cdot(z-0.5)\cdot(z-0.6)}{z-1}$$ to a difference equation, but I have no idea how.

Can someone help me?

Thanks,

Arjon

If $x[n]$ is the input of your discrete-time system and $y[n]$ is the output, then the transfer fucntion H(z) is written as:
$$H(z)= \frac{Y(z)}{X(z)}$$
where $$X(z)=Z(x[n])\; \;\; \;,\; \;\; \;\; \;Y(z)=Z(y[n])$$
So we get: $$\frac{Y(z)}{X(z)} = \frac{2(z-0.5)(z-0.6)}{z-1}\Rightarrow (z-1)Y(z)=2(z-0.5)(z-0.6)X(z)$$
$$\Rightarrow zY(z)-Y(z)=2z^2X(z)-2,2zX(z)+0,6X(z) $$
At this point we apply Inverse Z Transform and make use of both Linearity and Time Shifting properties.
So we end up with the following difference equation:$$\mathbf{y[n+1]-y[n]=2 \cdot x[n+2]-2.2 \cdot x[n+1]+0.6 \cdot x[n]}$$


Standard deviation of Morlet wavelet's Fourier transform MATLAB has the complex Morlet Wavelet in the following form:

$$\psi(t) = \frac{1}{\sqrt{\pi f_b}}e^{\frac{-t^2}{f_b}}e^{j2\pi f_ct}$$

I arrived at its Fourier transform as shown below (another arrived at the same: Fourier Transform of Morlet wavelet Function?).

$$X(f)=\mathcal F(x)(f)=e^{-\pi^2 f_b(f-f_c)^2}$$

I sought a relationship between $f_b$ and the standard deviation of the Gaussian form of $X(f)$. Is it right? Is there any citable published literature (book or journal article) arriving at this relationship?

Rearranging this into the familiar Gaussian form,
$$X(f) = e^{\frac{-(f-f_c)^2}{2\left({\frac{1}{\pi\sqrt{2f_b} })}\right)^2}}$$

The standard deviation $(\sigma_f)$ of this frequency domain Gaussian is then,
$$\sigma_f = \frac{1}{\pi\sqrt{2f_b}}$$
Rearranging it,
$$f_b = \frac{1}{2{(\pi\sigma_f)}^2}$$

I might have the ingredients to an answer here and here and here.


How does a rectangular window affect the signal? If i am not mistaken a rectangular window of length 10 gives us an array of 10x1:

1 1 1 1 1 1 1 1 1 1

And if we apply this to a signal we are left with exactly what we started with as:

x = x*1

So how does a rectangular window improve or do anything at all to our signal? What is the point of the function?

A rectangular window is typically used when you want to process a smaller slice out of a much larger vector of data, possibly zero-padding within a larger window (for better time locality, because a bigger FFT is too slow, doesn't fit in cache, or some other reason).

In the frequency domain, a rectangular window convolves a longer signal with a Sinc function, producing windowing artifacts that are sometimes called spectral "leakage" if the longer signal is not strictly/exactly integer periodic at the length of the window aperture.

A typical example might be that one wants to use an FFT of length 4096, zero-padded from a rectangular window of data of length 2048, in order to analyze a 43 millisecond chunk (one note? one phoneme?) out of a 48ksps audio recording that is an hour long.


Help filtering a logic level signal without degrading transitions I have a signal that is either noisy "low" or noisy "high", e.g.:

    ######     (high)
####      #### (low)
-------------- (zero)

I want to find transitions from low to high, precisely in time.

However noise is getting in a way:

I can't simply use a low-pass filter (transition accuracy lost)
I can't simply use a derivative (HF noise gets amplified)

I figure there must be a well-known algorithm, I just don't know it's name...

Edit:

Below is an example of "good" noisy step as well as "not so good, but solvable" step.

Edit:
Below is original signal (orange) and median filter (blue).

Because raw signal is periodic at high frequency, median filter produces an artefact of a ramp where I'd like to see a step function (or possibly two).

I've done a lot of work with things like this, hopefully the following will be of some help:

Since you've noted that the computation doesn't have to be performed in real time, this gives you plenty of options.  A natural choice, and I see someone has already suggested it, would be a median filter.  This usually works pretty well, so long as you choose the window length appropriately.  In your example that you showed, it looks like you may need to increase the window size (though it's hard to tell how many samples you have).  The main attraction to a median filter is that they will not smear edges, so long as you choose the window size appropriately; most standard smoothing methods will.

If you're over sampled, a legitimate option would be to just reduce your sampling rate and see where that gets you noise wise.

Other smoothing methods exist, but they'll typically smear your edges, and sort of ruin any edge/clock detection you'll want to do.  Since you have infinite computation time, median filters are probably your best bet in this realm.

Another option to check out would be a Savizky-Golay filter.  This is another data smoothing technique, but it attempts to model the signal using a polynomial.  MATLAB has a built-in that works great:

Savitzky-Golay Filtering in MATLAB

I've used this method on varying types of noisy signals, and if you're not constrained to running in real time, you can certainly tweak the parameters and hopefully find a combination that works for you.

Once you've got your data smoothed, a simple threshold would work I imagine.  You could statistically set this threshold based on the mean of the smoothed signal (this method assumes that the noise is fairly stationary).  Another method would be an edge detection filter, like say the first order derivative of a Gaussian or something.  You would convolve this filter with your smoothed signal, and then utilize some peak picking algorithm to get the start/stops of your transitions; this would be a little bit more robust than just utilizing some sort of differentiation filter due to the weights it implicitly applies.  You can get the filter in MATLAB by entering: diff(gausswin(sz,alpha)), naturally substituting the parameters for size and alpha that best work for you.

Another option that could forego smoothing with the right implementation is some type of hysteresis thresholding with a "hold" function built into it, given that you have some knowledge or expectation of what the values are when the the pulse isn't present.  This would involve having a double threshold, and creating some logic for deciding what values are above and what values are below.  The extra bit on top of a traditional hysteresis threshold would be that you'd need to specify some sort of holding mechanism, i.e. a segment is only valid if it breaks the top threshold for X many samples, and it is only invalid if it breaks the lower threshold for Y many samples.  This option would probably be the most work, as (from experience, I coded this for my own work fairly recently) it would take you a 100 or so lines of MATLAB code to implement.  Setting these thresholds is naturally the challenge here; I've found that a statistical approach is usually robust enough to work.

Last, but certainly not least, you could attempt to do some sort of adaptive filtering if you are able to characterize just the noise component of the signal, and said noise signal is stationary.  If you're able to capture just the noise in the signal, through an auxiliary sensor or some clever high pass filtering, you can use the Wiener-Hopf method to clean up the signal.  Here's an example:

By statistically characterizing the noise variance (not covariance in this case, since we only have one auxilliary signal) and cross-correlation, we're able to come up with an adaptive weight to remove the noise. For some background theory on this, consult Wiener Filter

Now at the signal is all cleaned up, a simple threshold would get you the edges.

If your noise isn't stationary, you can force it to be stationary in small windows through a process called time segmentation.  That's a bit involved though, and hopefully you won't need such a thing.

The MATLAB code that generated that plot is below:

up_s = zeros(1,250)+3;
d_s  = zeros(1,500)+3;
u_s  = ones(1,500)*10;
s = [up_s,repmat([d_s,u_s],[1,50])];
t = linspace(0,1,numel(s));
% set up the noisy signal
nz = sin(2*pi*500*(t));
snz = s+ nz;

% w = R^-1 r, from wiener-hopf method
% zero mean assumed for covariance and correlation;
% k terms cancel and are excluded
w = inv(nz*nz') * nz*snz';
% apply the weight to the signal
w_sig = snz - w'.*nz;

figure;
ax1 = subplot(2,1,1);
plot(snz);
title('Original Signal');
xlabel('Sample'); ylabel('Magnitude');
xlim([1 numel(snz)]);

ax2 = subplot(2,1,2);
plot(w_sig);
title('Signal Cleaned with Wiener-Hopf Filtering');
xlabel('Sample'); ylabel('Magnitude');
xlim([1 numel(snz)]);

linkaxes([ax1 ax2],'x');

tl;dr:
Median filtering is great, you just need to select the parameters appropriately for your sample rate. Savitzky-Golay filtering will also get you good results, given that you get the parameters right, which can be challenging. With any smoothing method, you can expect some smearing of the edges, but median filters really minimize this.

If you can characterize just the noise of the signal, you can use an adaptive filter technique like the Wiener-Hopf method; this method assumes that the noise is stationary, so that's something to keep in mind. With this method, you'll get essentially zero smearing of the clocking edges.


Explanation for the quantization matlab function I found this code on quantization on page 314 in the book Modern Digital and Analog Communication Systems by B.P.Lathi and Zhi Ding. I can understand only half of the code. Please, someone, help in explaining the code.

Code:-

function [q_out, Delta] = uniquan(sig_in, L)
% L = number of uniform quantization levels
% sig_in - input signal vector
% Function Output
%     q_out - quantized output
%     Delta - quantization interval

 sig_pmax = max(sig_in);              % finding the positive peak
 sig_nmax = min(sig_in);              % finding the negative peak
 Delta = (sig_pmax - sig_nmax)/L;     % quantization interval
 q_level = sig_nmax + Delta/2 : Delta : sig_pmax - Delta/2;  %defining q-levels

% Till here I understand the code. I am not able to understand from here.
% Mainly not able to understand what is sigp doing

 sigp = (sig_in - sig_nmax)/Delta + 1/2;   % convert into 1/2 to L+1/2 range
 qindex = round(sigp);                     % round to 1,2,3..... L levels
 qindex = min(qindex,L);                   % Eliminate L+1 as rare possibility
 q_out = q_level(qindex);                  % use index vector to generate output

end

Ebook link - https://www.academia.edu/36135973/Book_Modern_Digital_And_Analog_Communication_Systems_4th_edition_by_Lathi

sigp = (sig_in - sig_nmax)/Delta + 1/2;   % convert into 1/2 to L+1/2 range
qindex = round(sigp);                     % round to 1,2,3..... L levels
qindex = min(qindex,L);                   % Eliminate L+1 as rare possibility
q_out = q_level(qindex);                  % use index vector to generate output
clear
end

sigp is a mid-tread uniform quantizer which is defined using the quantization step size Delta whose details you can find here: https://en.wikipedia.org/wiki/Quantization_(signal_processing)#Mid-riser_and_mid-tread_uniform_quantizers

sigp uses sig_in - sig_nmax instead of sig_in so that the signal can lie between 0 and 2 instead of -1 and 1 since indexing in MATLAB is always greater than 0

1/2 is added to round the first index to 1.

qindex = min(qindex,L);

The highest level can be L so if qindex goes over L, it is truncated to L

q_out = q_level(qindex);

q_out gives the quantization level corresponding to the qindex.


ROC of Z transform Doesnt include a pole on the boundary? I cannot figure out what is going on here. I have an example problem in my book that says the ROC of a certain function is $$ 0.5 < |z| < \infty$$
The function's denominator is $$ 1 - z^{-1} + 0.75z^{-2} -0.25z^{-3} + 0.0625z^{-4} $$
By properties of the ROC, there should be a pole on the boundary of the ROC. Since the example given was a right sided signal, the ROC will be the area outside of a circle boundary. So, the signal should have a pole on the boundary of that circle, which has a radius of 0.5 correct?
But when I plug in 0.5 to the denominator equation, I get 1, not zero, which makes 0.5 not a pole in this case. I have tried graphs, and placing different parenthesis, and calculated it multiple times, but the denominator equation above never reaches zero. According to ROC properties, it has to though!
Where am I thinking wrong here?
Thanks,
-Dom

The ROC $0.5<|z|<\infty$ does not imply that there's a pole at $z=0.5$. What it does say is that there's at least one pole satisfying $|z|=0.5$ (and no other pole with a radius larger than $0.5$), and this is also the case for the given denominator polynomial. The roots of the polynomial are
$$z_{1,2}=0.5\, e^{j\pi/3}\tag{1}$$
and
$$z_{3,4}=0.5\, e^{-j\pi/3}\tag{1}$$
I.e., there are two double poles, and all four poles lie on a circle with radius $0.5$.


Hardware Transfer Function Estimation Suppose you have the ability to inject any arbitrary waveform into a piece of analog rf hardware and collect and digitize the output for analysis.

If you wanted to characterize/estimate the transfer function of the hardware, what waveforms would you use?

You can use adaptive filters for this task. Feed the hardware with uniform white noise and collect the output. Then feed the adaptive filter with exactly the same noise sequence and use the output of the hardware you recorded as the desired response of the adaptive filter. The filter can then be trained with standard algorithm such as LMS. Of course, you will not get the exact transfer function, but the trained adaptive filter's transfer function should be close enough to serve as model of the hardware in question.


Extracting audio track to determine the song beat I'm completely new to this field.

Perhaps I'm calling it by a wrong name, but I want to analyze an audio track and retrieve a timestamp to beat mapping.

Is there an open source 3rd party library that can perform this, or could/should I do it myself?

Sonic Annotator. Good: their beat tracking algorithm is the most sophisticated. Bad: this is not easy to integrate and interface into your own application (you need to load their own flavor of plug-ins called "VAMP"); and the beat tracker is closed source - available as a pre-compiled binary lib.

Aubio. Good: very easy to integrate, plain C with few dependencies. Bad: the quality of the beat-tracking is so-so.

Marsyas. Good: beat-tracking accuracy is OK (but not as good as Sonic Annotator's). Bad: many dependencies, very complex framework.

Depending on what is your target application (is the audio to analyze a recording of a well-known/distributed song, and does your system have internet access), you could also use the Echo Nest APIs.

If you want to roll your own, I suggest you this as a starting point - start porting this rather simple matlab implementation to your language of choice, and you'll be at 80% or 90% of the accuracy of the most sophisticated solution.


Algorithm to convert a GeoTIFF file to PNG image format I have a bunch of GeoTIFF files that I am processing into a TFrecords format for machine learning. Now I can't put tiffs directly into a TFRecords format and so I have to convert the files to either PNG or JPEG. However, I did not want to generate a big bunch of converted PNG files to clutter up my computer and was hoping to do the conversion of the data in-memory.

The GeoTIFF has only 3 bands--it is RGB PanSharpened--so no problem with selecting out the bands I need. When I imported a sample GeoTIFF file into python using rasterio I found that the range for pixel values was between [282, 1560], which seemed rather high. This is satellite imagery, BTW. Now PNG or JPEG files seem to have pixel ranges of between [0, 255] or [0,1], from what I gather. This seems like a simple rescaling of the values to my naive eyes.

So my question is whether anyone could provide the basic algorithm to convert TIFF file pixels into PNG or JPEG pixels. I don't need to worry about the spatial coordinates in this case, so I can treat this as a straight TIFF --> PNG or JPEG conversion.

Could I just divide every TIFF pixel by the max in the range, like 1560? I don't know what the maximum value for a TIFF pixel is, but I could divide by the max otherwise. Or I would have to divide by a factor that maps the range [282, 1560] --> [0, 255], right. So that is just a simple scalar multiplication.

I looked around the web but no one seems to have a basic algorithm for this. I read the file specification for TIFFs and PNGs but could not find the maximum and minimum pixel ranges. All I can find is either people converting the images to new files--using GDAL or gdal_translate or in some cases using PIL. Or other answers fall into using GUI tools like QGIS to again save to new files. But like I said, I would rather not save to a file.

Hopefully someone might know the algorithm to do this conversion manually. Thanks.

I can't put tiffs directly into a TFRecords format...

You can. TFRecords does not know what an "image" is, it simply stores a string of bytes. What this means to higher level applications is up to them.

When I imported a sample GeoTIFF file into python using rasterio I found that the range for pixel values was between [282, 1560], which seemed rather high.

Depending on the kind of data they contain, GeoTIFF images can have higher colour depth per pixel. For example, the SRTM data provide a digital elevation model of the Earth. Each pixel value is associated with the height of a "square" that is 30 by 30 meters (along the equator). The variation of that height over the whole planet is much larger than 256. SRTM's GeoTIFFs have a 16bit colour depth.

Now PNG or JPEG files seem to have pixel ranges of between [0, 255] or [0,1], from what I gather.

PNG and JPEG can save 16bit per pixel colour depth images and you can still work with them with libraries such as osgeo or GIMP.
Depending on the content of the GeoTIFFs though, the data in PNG format (lossless encoding) might not compress as much as you might expect and with JPEG (lossy encoding) you are going to have to tune the compression to avoid losing any necessary features. In the end, it might be better to just leave the data in GeoTIFF or read them in a numpy array and then serialise that straight into TFRecords.

This seems like a simple rescaling of the values to my naive eyes.
So my question is whether anyone could provide the basic algorithm to convert TIFF file pixels into PNG or JPEG pixels. I don't need to worry about the spatial coordinates in this case, so I can treat this as a straight TIFF --> PNG or JPEG conversion.

Indeed, it has to be a simple linear scaling but this is not a good idea for a number of reasons.
If the absolute values are of importance (for example) then the 8bit dynamic range is simply not adequate.
If your subsequent processing is about learning gradients (for another example), then you are going to have to normalise to the global extrema (not the maximum / minimum of each individual file) so that the relative slopes do not change. The danger in doing so is again dynamic range. Although it is possible to squash any numerical range in the interval $0 \dots 1$, when this interval is quantised in 8 bits it is only divided into 256 "steps". If two pixel values differ by less than $\frac{1}{256}$, then they will be assigned to the same "bin". Also, if the dataset contains "spikes", that is the majority of the values are between 20-100 (for example) but there is one value at 10000, then these spikes will "kill" the contrast when you normalise at 8bits. What you will get is a mostly black image with 1 white pixel. Again, it might be possible to remove those spikes but if you have to work with them (e.g. for a multispectral application), then normalising to an 8bit colour depth is not a good idea.
The other thing you would be losing if you went with this option is the transform data. That is, the transform that converts coordinates from image (pixel positions and values) to world coordinates (longitude, latitude, height).
Suppose that your GeoTIFF image ($I$) describes a digital elevation model that has both negative and positive values. The total range of the GeoTIFF would be $max(I)-min(I)$. If $min(I)$ is not zero then your $0$ point in the converted data will be "floating". It might be 20 on one image, 42 on another 100 on yet another and so on.
The transformation between two ranges is straigthforward: $O_{x,y} = O_{min} + \frac{I_{x,y}-min(I)}{max(I)+min(I)} \times range(O)$, where $O$ is the output image, $O_{min}$ is the new desired minimum and $range(O)$ is the desired range.
While the transformation is straightforward, the impact it has on the data and subsequent models can be severe, therefore my recommendation would be to preserve the original numbers as much as possible and only intervene if you are absolutely sure that the transformation will not affect your objectives.
Hope this helps.


Which eigenvector to choose after calculating PCA After calculating the principal component analysis (PCA) of a given data set, we are normally left with a matrix containing the eigenvectors sorted in order of the size of the eigenvalues. Now, in pattern recognition which eigenvector should we choose: the first eigenvector or do we have to do further processing of the eigenvector in order to choose the desired eigenvector?

The eigenvectors are in the order of explained variance. Usually one uses PCA for dimension reduction, that is, selects some N first components. The problem of choosing appropriate N is called dimension estimation. Easiest way would be to take such N that the selected PCs explain 90% or 95% of the variance in the original data. Alternatively, you could use some information-theoric method such as minimum description length (MDL).


Can you get flat response with Linkwitz-Riley cross-over? Many pictures display that the Linkwitz-Riley has $-6\textrm{ dB}$ dips at the cross-over frequencies.

Are these fixed or can these made flat?

Linkwitz Riley filters always have  -6 dB at the cut off frequency. The whole idea of a cross over is to split the signal into multiple bands. The advantage of a Linkwitz Riley is that it adds to be flat (at least in amplitude). So each the highpass and the low pass have -6 dB at the crossover, but the sum of both is flat.


Why do digital modulation when the signal entering the channel is analog? I am learning about digital filters and came across the link Why do digital filters work?

In real world communications, I think the input to the channel is in analog after being modulated by digital modulation techniques such as BPSK, Qam etc. These modulation techniques convert the real valued signal to a symbol space. The following questions may seem trivial and can be found in text books, but it is difficult to extract specific information.

I am unable to understand, is the input to the channel real valued numbers or symbols?
Is the modulated signal passed into the channel? IS the modulated signal digital?
If not, Why should we do digital modulation if the signal entering the channel is analog and in which stage do we need digital filters?

An example with a block diagram would be helpful. Thank you

I think part of your confusion is actually between discrete-time and digital signals. Some textbooks do not clearly distinguish between them.

All signals that exist in the physical world are real. So, the signals that are produced by any transmitter, and the signal that is picked up by the receiver, are real.

In communications, we say we do "digital communications" when the analog transmitted signal conveys "symbols", which are actually numbers drawn from a finite set. This is different than "analog communications", such as AM and FM radio, where the traansmitted information is a (audio) waveform that can take arbitrary values.

In addition to this, these days communications systems are designed as much as possible using digital signal processing techniques. The nature of these techniques is that they deal with discrete-time (sampled) signals. So, a modern software-defined radio receiver will have an analog antenna, a few analog filters and amplifiers, and an analog quadrature downconverter; then, the analog signal is discretized (sampled) and handed off to a digital signal processor for further processing. The digital signal processing can happen on specialized hardware (such as an FPGA or ASIC), or as software running on a DSP or on a general-purpose CPU. This is where digital filters would be used.

Note that a receiver can be implemented using DSP regardless of whether the underlying communications technology is analog (AM, FM, old TV, etc) or digital (BPSK, QPSK, etc).


Different phase response plots using MATLAB's $\tt freqz$ and $\tt fvtool$? I have a butterworth filter in MATLAB:

[B, A] = butter( 2, 3000/16000*2, 'high' );

I plot it in two ways:

using freqz:[h,w] = freqz(B,A,8000);
plot(w);
and using fvtool:h = fvtool(B,A);
h.Analysis = 'phase';
h.Fs = 16000;
h.FrequencyRange='[-Fs/2, Fs/2)';

First plot gives me linear phase, and second one not.

Why do I get different plots?

Citing the help output of freqz:

[H,W] = freqz(B,A,N) returns the N-point complex frequency response
vector H and the N-point frequency vector W in radians/sample of
the filter.

So, your output is not amplitude and phase, but (complex) transfer function and frequency. Hence, to plot the phase using freqz you need to use

plt(w, angle(h));


Maximum likelihood estimator of active time delay and passive time delay A typical time delay estimation problem has the model:
$$
\begin{align}
     x_1(t) &= s(t)+ n_1(t) \\
    x_2(t) &= a s(t-D) + n_2(t)
\end{align}
$$
Where $n_1$ and $n_2$ are considered to be Gaussian noise.
An active time delay estimation requires $n_1$ to be $0$ and $a=1$.

The maximum likelihood estimator, in this case, is the matched filter or the basic cross-correlation method for finding delay. In passive case, $n_1$ is not equal to 0, and the maximum likelihood estimator is HT processor in Generalized Cross Correlation (GCC).

Can someone provide me with an intuition why in active time delay estimation case, the matched filter performs favorably better than GCC, even if I know the spectra of my signal (i.e. the spectral estimates are accurate enough)?

I would expect both GCC and simple CC to give almost equal results. And also, my results vary when estimating time delay by cross correlation in time and frequency domain. Is there any explanation behind this?

Very intuitively, the Generalized Cross-Correlation is a "standard" cross-correlation of the windowed signals (I'll restrict myself to window-GCC, I'm pretty certain there's others, too!). Windowing happens to increase the "peakiness" of the cross-correlation.

It's basically down to the same trade-off between temporal and spectral precision that you consider when applying a window for the DFT. You trade energy in the "temporal edges" (read: start and end of the signal are attenuated by the window) for less leakage in the spectrum.

Now, back to your problem. In the active TDE definition you use,

$$\begin{align}
R_{x_1x_2} =& E\left\{x_1(t)x_2(t+\tau)\right\}\\
=& E\left\{s(t)s(t+\tau)+ s(t)n_2(t+\tau)\right\}\tag1\\
=&\underbrace{E\left\{s(t)s(t+\tau)\right\}}_\text{what contains the peak we're looking for}+ \underbrace{E\left\{s(t)n_2(t)\right\}}_\text{measurement noise}\tag2
\end{align}$$

Now, your GCC with a window $w$ looks something like (how $w$ is determined differs for different GCC approaches)

$$\begin{align}
\text{GCC}_{x_1x_2} =& E\left\{w(t)x_1(t)\cdot w(t+\tau)x_2(t+\tau)\right\}\\
=& E\left\{s(t)s(t+\tau)w(t)w(t+\tau)+\\
 s(t)n_2(t+\tau)w(t)w(t+\tau)\right\}\tag3
\end{align}$$

The step from $(1)$ to $(2)$ is intuitive, since the noise is uncorrelated to the signal we sent and received. (I'm assuming that, but it's a justifiable assumption given the problem statement.)

Now, in $(3)$ I'm tempted to say "ok, lets just pull $w(t)w(t+\tau)$ out of the expectation operator, right?", but I mustn't do that, since usually GCCs choose a $w$ based on $x_1$ and $x_2$; in other words, $w$ depends on $s$ and mustn't be treated like a deterministic function!

Information-theoretically, this means you introduce entropy into your $s(t)s(t+\tau)$ estimate by multiplying it with something that was estimated based on a measurement that includes noise. So, while $(2)$ only incorporates the noise-signal cross-correlation once, and since signal and noise are uncorrelated, with increasing length of observation, the error term in $(2)$ converges to zero.

It's not that easy for $(3)$, especially without a specific noise model, but at the very least we can say that it should (by sheer arrogance and the intuitive application of the Cauchy-Schwarz Inequality), converge slower, if it converges at all towards the actual cross-correlation.

An oft-used practical GCC is GCC-PHAT, where $w$ is a function (usually: a reciprocal) of the Fourier Transform of an PSD estimate of $x_2$, which is pretty close to being the autocorrelation of $x_2$. Now, a reciprocal of an estimate that, should the observed duration of signal and noise be very long and if both the signal and the noise being rather whit, should contain a lot of values close to zero, is obviously rather instable, so this is an especially questionable choice in the case that $x_1$ doesn't contain noise of its own.


How to Find the Probability of Kalman Filter States? I am working on a video object tracking problem. I am using Kalman filter to predict and correct the object position return by an algorithm such as CamShift. I want to adjust the likelihood probability of the states so as to customize the weight of the states before feeding to the Kalman filter. Currently this the code I have:

#Define kalman filter
kalman = cv2.KalmanFilter(4, 2, 2)
# initial position: (x0,y0)
state = np.array([col+w/2,row+h/2,0,0],dtype='float64')
kalman.transitionMatrix = np.array([[1.,0., .1, 0.],
                                [0., 1., 0., .1],
                                [0., 0., 1., 0.],
                                [0., 0., 0., 1.]])
kalman.measurementMatrix = 1. * np.eye(2, 4)
kalman.processNoiseCov = 1e-5 * np.eye(4, 4)
kalman.measurementNoiseCov = 1e-3 * np.eye(2, 2)
kalman.errorCovPost = 1e-1 * np.eye(4, 4)
kalman.statePost = state

Here, 'row' and 'col' are the top-left row and column indices, respectively of the current rectangular box containing the object position, and 'w' and 'h' are the width and height, respectively of the rectangular box containing the object.

Now, suppose (x1,y1) is the position returned by my algorithm (e.g. CamShift,Optical flow etc.).

Now I pass (x1,y1) and the explicitly calculated velocity (v) and acceleration(a) to Kalman filter as:

#Predict using kalman filter
prediction = kalman.predict()

# Update kalman state using optical flow results
pts=np.array([np.float32(x1),np.float32(y1),np.float32(v),np.float32(a)],
dtype=np.float64)

state=pts
measurement = (np.dot(kalman.measurementNoiseCov, np.random.randn(2,
1))).reshape(-1)
measurement = np.dot(kalman.measurementMatrix, state) + measurement
# use optical flow to correct kalman filter
prediction_new=kalman.correct(measurement)

It is all fine till here. But what if I want to give only less than 0.5 probability weightage to the previous state (=[x0,y0,0,0]) compared to the present state (=[x1,y1,v,a]) while computing 'prediction_new'. I have extensively searched the Kalman filter OpenCV library as well as online, but couldn't find an answer. Appreciate help on this. Thanks in advance.

The idea of Kalman is it fuses the data with the optimal weighing to get results.
Optimal requires definition and in the case of Kalman it means that if your model is right this is the optimal.

Since we operate Kalman filter under the assumption the model is right no Kalman Filter will allow you doing what you want.

Now, it seems that in your case what you don't like is the transient phase when there are low number of measurements.
What I can suggest is that for this start phase you will built your own filter based on your weighing and after few iterations send your result to the Kalman Filter as a starting point.


What is the correct block diagram of CDMA communication system? I was reading someone's thesis and now i am more confused about CDMA system. From what i understood is CDMA is multiplexing which allows multiple signals to transmit over single channel. so general structure of CDMA communication system should be like shown below:
Binary Data --> Multiply PN Code --> XOR(spreading) -->BPSK/QPSK Modulation --> RF End
                                                                                |
                                                                                Channel
                                                                                |
Restore Data <-- Despreading <-- Aquistion and lock <-- RF Demodulation <--    RF Reception

And Now I had come across block diagram in the thesis, Pg.27 which basically modulate and they do spreading.

From my understanding it should be spreaded before doing modulation. This is my understanding, please correct me if i am wrong.

But can someone please confirm which one is right ? It would be great if you explain a little bit.

Since spreading and symbol modulation are memoryless linear operations, the order makes absolutely no difference. Both methods of writing the same thing down are equally valid, and the one you choose mainly depends on what is convenient for the rest of your text.


System Identification with a Limited Bandwidth Input Signal and Region of Interest Given a FIR filter $h[n]$. Its action can described as:

$$
\mathbf{y} = \mathbf{H} \mathbf{x} \\
\mathbf{y} = \mathbf{X} \mathbf{h}
$$

where $\mathbf{H}$ and $\mathbf{X}$ is a Toeplitz matrix.
If $h$ is unknown, Least Squares with a white Gaussian input signal $x[n]$ can be used to find the unknown coefficients:

$$
\hat{\mathbf{h}} = (\mathbf{X}^{T}\mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}
$$

Caveat: $x[n]$ must be white; otherwise the regression matrix $\mathbf{X}^T \mathbf{X}$ is badly conditioned.

The frequency domain information is encoded in the coefficients $h$. However, as can be seen above, the LS algorithm estimates the coefficients with zero prior knowledge; the estimation ONLY depends on the input signal $x$. It does not matter if the system to be identified is an allpass filter or has a notch of 200dB attenuation at $\pi/2$.

Now my question: What do I do if I only care about a small frequency range in $h$ and hence my input signal $x[n]$ does not need to be white?

Example: My Nyquist rate is 10kHz. My unknown system is a lowpass with -3dB at 300 Hz. It has some "weird" frequency behavior around 300 Hz which I want to estimate. I do NOT care about anything beyond, say, 500 Hz.
Additionally, my measurement setup prevents me from using a white input signal. I have a bandwidth limitation of 500 Hz. I cannot change the Nyquist rate.

With Least Squares I cannot identify the system because $x$ is not white (persistently exciting). Regularization/SVD does not help me: It provides a biased solution and still gives me $h$ values that try to estimate the entire frequency range. But I really want to say "Give me the $h$ that describes the unknown system best up to 500 Hz with a 500 Hz input signal"

The reason x[n] must be white is because the solution will effectively spectrally weight the channel response based on the amount of energy present in each spectral frequency location. A white noise source provides equal weight to all frequencies. If energy is not present in any particular frequency bin, a proper solution cannot be found for that frequency.

If you only care about a small band of signals, then I would argue that you can still use the least squared approach. The reason for this is to consider a system that provides band pass filtering: for such a system I could excite the input with a white noise source and use the least squares approach (Wiener-Hopf equations) of comparing the input and output signals which would accurately provide a least squares estimate of the channel. So if you have a band limited signal, as long as that signal is white over your band of interest, it will still provide the accurate solution over that band (and you ignore everything else).


How to get main frequency from digital signal I have sampled digital signal with noise. I want to get main frequency of this signal how to do it?
Currently it is smoothed by moving average. Is good way to make Fourier Transform an take the highest peak? Or the first peak?

Thanks.

Generally, yes, you should perform FFT and take the highest peak from the first half of the transformed signal. However, if the frequency your are looking for is constant (or close to constant), i.e. doesn't change (much) throughout the observation time, I would recommend using Welch's method instead of a simple FFT.
If the accuracy is not good enough you can perform interpolation around the peak, or on the entire signal using zero-padding.


What are the limitations of using Magnetic induction rather than acoustic communication in underwater communication In underwater environment, Why acoustic underwater communications are still famous more than magnetic induction communication however in magnetic induction underwater communication, higher data rate is offered, shorter propagation delay and highly environment-independent channel behavior?
In other words, what are they disadvantages of using magnetic induction communications in underwater environment?

Actually, magnetic induction is a current research topic for underwater communications, see http://bwn.ece.gatech.edu/papers/2015/j14.pdf. Also, this paper includes a nice comparison of underwater communication strategies, including electromagnetic waves, acoustical communication, and optical communication.

One major disadvantage of magnetic induction is the short communication range below 100 meters. Acoustical communication has a range that is much larger in the region of kilometers.


Comparing Welch and FFT power spectrum density Based on previous answers from the forum, I implemented a function to compute the Power Spectrum of a 1D time series.

def pow_spect(x, fs):
    nt = len(x)
    power = np.abs(np.fft.rfft(x, nt))
    ff = np.fft.rfftfreq(nt, 1. / fs)
    power = (power ** 2) / (nt ** 2)
    power *= 2.
    power[:, 0] /= 2.
    if nt % 2 == 0:
        power[:, -1] /= 2.
    return ff, power

I compare the output of this implementation to the estimation of the power spectrum density given by scipy.welch :

nt = len(x)
ff_welch, pxx = signal.welch(x, fs, window=signal.get_window('boxcar', nt), scaling='spectrum')

Eventually, I compute ff, p = pow_spect(x, fs) and I compare : p_db_1 = 10. * np.log10(p) and p_db_2 = 10. * np.log10(pxx) for a given time series x. I noticed that all the coefficients of p_db_1 and p_db_2 are very close (up 10^(-8)) except for the first one, p_db_1[0] and p_db_2[0]. Did I make a mistake somewhere ? Thank you for your help.

scipy.signal.welch will, by default, detrend the data by subtracting the mean of each segment. This DC level suppression will only affect bin zero of the resulting spectrum. This is very useful since the DC offset can be so much larger than the noise that it wouldn't fit on a reasonably scaled log plot.


Relationship Between Sampled Continuous and Discrete Time Signals Consider the sketched system below. $x_c(t)$ is an arbitrary, continuous-time signal at the input and $s(t)$ is an impulse train, defined as $s(t)=\sum_{n=-\infty}^{\infty} \delta(t-nT)$, where T is the sampling period. Hence:

$$ x_s(t) = x_c(t)s(t)=\sum_{n=-\infty}^{\infty} x_c(nT)\delta(t-nT) $$

Now I calculate and compare the fourier transform of $x_s(t)$ and $x_c(nT)$.

$$
X_s(j\Omega)=\int_{-\infty}^{\infty} x_s(t)e^{-j\Omega t}dt=\int_{-\infty}^{\infty} \sum_{n=-\infty}^{\infty} x_c(nT)\delta(t-nT)e^{-j\Omega t}dt=\sum_{n=-\infty}^{\infty}x_c(nT)e^{-j\Omega nT}
$$

$$
X(e^{j\omega})=\sum_{n=-\infty}^{\infty} x[n]e^{-j\omega n}
$$

The results look quite similar. In fact, they are equal if $\omega=\Omega T$.

The calculation is no problem, but the issue I have is understanding what this really means. So obviously it is a scaling. If I define the angular frequency as $\Omega_s=\frac{2\pi}{T}$ and plug it into the above formula for $\omega$, I get $\omega_s=2\pi$. So it seems that through the process of discretization the actual frequency information is lost and the fourier transform goes from 0 to $2\pi$. In order to retrieve the information, I need to multiply $\omega$ with the sampling frequency. Is this statement correct? Why does this actually happen? What meaning is behind this? Probably I am missing an obvious thing, I feel like missing the forest for the trees.

You've shown that the (continuous-time) Fourier transform (CTFT) of a sampled continuous-time signal equals the discrete-time Fourier transform (DTFT) of the corresponding discrete-time signal. In both cases, the spectrum is periodic. The independent variables are related by

$$\omega=\Omega T\tag{1}$$

where $\Omega$ is the angular frequency ($\Omega=2\pi f$), $T=1/f_s$ is the sampling period ($f_s$ is the sampling frequency), and $\omega$ is the normalized angular frequency, which is often used as the independent frequency variable when describing spectra of discrete-time signals. Note that you could also use $\Omega$ as the independent variable for discrete-time signals, but it always appears together with the sampling interval $T$. This is an important property of sampled or discrete-time signals, and that's also why the two variables $\Omega$ and $T$ are commonly combined into a single variable $\omega=\Omega T=2\pi f/f_s$.

Note that sampling results in a periodization of the spectrum, and this is why information can be lost. This is true for both spectra $X_s(j\Omega)$ and $X(e^{j\omega})$. The periodicity of the spectra implies that the fundamental frequency interval $f\in [-f_s/2,f_s/2]$, or, equivalently, $\Omega\in [-\pi/T,\pi/T]$ or $\omega\in [-\pi,\pi]$ contains all information of the sampled signal, and all frequencies outside that interval are redundant. This is an inherent property of sampled (or  discrete-time) signals, and it is the cause of aliasing, i.e. the mapping of different continuous-time signals to the same sampled signal.


Concrete use of convolution algebraic properties There are 3 main properties of convolution in every image processing lecture notes I read :

Commutativity: $f\star h = h\star f$
Associativity: $f\star (h_1\star h_2) = (f\star h_1)\star h_2$
Distributivity: $f\star (h_1+h_2) = f\star h_1 + f\star h_2$

I'm looking for examples of when these are used (in image processing preferably).

From my research:

A system with unit sample response $h$ and input $f$ behave the same way as a system with unit sample response $f$ and input $h$.

$\Rightarrow$ We can interchange filter and image when computing convolution.

If two systems with unit sample responses $h_1$ and $h_2$ are connected in cascade an equivalent system is one that has an unit sample response equal to the convolution of $h_1$ and $h_2$

$\Rightarrow$ We can convolve our different filters instead of applying them consecutively. It can save time if some filters are often used together.

If two systems with unit sample responses $h_1$ and $h_2$ are connected in parallel an equivalent system is one that has an unit sample response equal to the sum of $h_1$ and $h_2$

This is where I'm unsure. In which instances we would frequently have multiple systems connected in parallel ?
Certainly in ML while computing some CNN.
Also learned about some sharpening method using blur:
$$
f' = f + c \cdot (f - f * h_{\text{gaussian kernel}})
$$

I can't think of any more common examples so I would be glad if you could enlighten me with some and/or correct my eventual mistakes.
Thank you very much!

For part 3, there is an old trick dating from the time where computing was expensive. If you have a smoothed version $I_\sigma$ of a picture $I$, you can get a crude yet simple approximation of a  Laplacian from $I -I_\sigma$, illustrated here with the Gaussian:

It is indeed generalized in many unsharp filters, with other weights, or in multirate/multiscale (pyramid) operations.
For part 2, you can convolve twice by Gaussians with a single convolution of a single Gaussian. Yet the most classical (yet rarely understood) is that smoothing first then differentiating is the same as differentiating then smoothing, see for instance Laplacian of Gaussian operator. Many people believe that smoothing first reduces noise. However, one can combine derivatives and smoothing, and save a couple of operations, and even design innovative gradients:

For part 1, this is used a lot in proofs. This is something one uses implicitly very often, for instance in cases when you have an image that you need to convolve a lot of times (scale-space filtering for instance) vs a lot of images you convolve when a filter known a priori.
Related question:

How to choose the size of a Laplacian of Gaussian kernel for filtering images fast?


How is the Leslie speaker simulated? What are Leslie speaker / "rotary" style of digital audio processors doing?

I think a simple one would model:

Delay between speaker and a pair (stereo) of microphones.
The tonal change when the speaker rotates (LFO modulated filter and amplitude modulation).

But what (more) is there really?

I highly recommend taking a closer look on great book:

  Zolzer U. - DAFX Digital Audio Effects

As well as his publication on this specific topic:

  Disch S., Zolzer U. - Modulation and Delay Line Based Digital Audio
  Effects

Leslie (Rotary) speaker can be simulated using tremolo and Doppler effect for both change in signal amplitude and pitch.

Basic signal flow diagram that can be found in the article:


How to sketch/plot this discrete time signal? $(n+5)u[n+5]+2nu[n]+(n-10)u[n-10]$ is the question.

The $u[n]$ is messing me up. If I substitute $n$ as $0$, what would be $u[5]$ or $u[-10]$?

Begin with defining the base signal as $x[n] = n ~u[n]$ and recognize the output as a linear compositions (superpositions) of the base signal as:
$$y[n] = x[n+5] + 2 x[n] + x[n-10]$$

Now what you have to answer is the result of translations and scalings applied on the basic signal $x[n]$. This ability is a simple yet an indispensibly necessary skill for a clear understanding and successful applications of all subsequent signal processing theory. Therefore you must train yourself enough for this task.

Incidentally, the signal $x[n+5]$ is an advanced (shift left) version of $x[n]$ whereas $x[n-10]$ is a delayed (shift right) version.

Graphically the result is obtained by summing the three components drawn on the on the same plot. A more detailed approach would use the piecewise definitions of the components to define the sum. It would proceed as follows:

Let $x[n]$be defined as:
$$ x[n] = n ~\text{u}[n] =
\begin{cases}
n , \text{ for } n \geq 0 \\
0 , \text{ for } n < 0 \\
\end{cases}
$$

Then you can define the three superposed signals as:

$$
\begin{align}
x_1[n] = x[n+5] = (n+5) ~u[n+5] &= \begin{cases} n+5, &\text{ for } &n+5\geq 0 \\ 0, &\text{ for } &n+5< 0 \\ \end{cases} \\
x_2[n] = 2x[n] = 2n ~u[n] &= \begin{cases} 2n,~~~~ &\text{ for } &2n \geq 0 \\ 0,~~~ &\text{ for } &2n < 0 \\ \end{cases} \\
x_3[n] = x[n-10] = (n-10) ~u[n-10] &= \begin{cases} n-10, &\text{for} &n-10 \geq 0 \\0, &\text{for} &n-10<0 \\ \end{cases} \\
\end{align}
$$

Modifiy the intervals:
$$
\begin{align}
x_1[n] = x[n+5] = (n+5) ~u[n+5] &= \begin{cases} n+5, &\text{ for } &n \geq -5 \\ 0, &\text{ for } &n<-5 \\ \end{cases} \\
x_2[n] = 2x[n] = 2n ~u[n] &= \begin{cases} 2n,~~~~ &\text{ for } &n \geq 0 \\ 0,~~~ &\text{ for } &n<0 \\ \end{cases} \\
x_3[n] = x[n-10] = (n-10) ~u[n-10] &= \begin{cases} n-10, &\text{for} &n \geq 10 \\0, &\text{for} &n<10 \\ \end{cases} \\
\end{align}
$$

Therefore $y[n] = x_1[n] + x_2[n] + x_3[n]$ is:
$$
y[n] =
\begin{cases}
0 , &\text{ for }     &-\infty<n<-5 \\
n+5 , &\text{ for } &-5 \leq n <0 \\
3n+5 , &\text{ for } &0 \leq n <10 \\
4n-5 , &\text{ for } &10 \leq n <\infty \\
\end{cases}
$$


Wavelet signal analysis for ultrasound echo I am new to the wavelet field, and I am working with ultrasound signal. For my study, I needed to measure material using the pulse echo reflection method.

Basically, when the transducer transmits the signal, the signal will reflected back to the transducer. If material has a thickness, the signal will reflected for two sides which are the front and back. After receiving both signals, I need to differentiate both of them, however the received signals are overlayed.

I think by using a wavelet method, I can differentiate the signal, but the result that I got is unclear.  Can anyone explain about this result?  Is it that both signals can be differentiated?

From what I can see you are not trying to detect defects. For defect detection  (CM,Inspection) wavelets are really popular and they had and still has been used extensively. I am working in the same field for a while but from what I can understand, you are trying to measure the thickness of the material.

In any case, I couldnt quite understand your problem but I strongly advise you to go for spectrogram first. Scalogram is inherently better in terms of resolution (both in time and frequency) but it might be a bit tricky to understand scale-frequency relation as opposed to spectrogram.

Based on your plots, wavelet transform is picking up quite a few of wavelets. This probably is due to the fact that your signals are overlayed. You might want to seperate them. You might also want to go for different wavelets and decomposition levels, I personally know mexican hat and db mother wavelets can provide good results for NDT applications, both for denoising and transforms.


Frequency response and impulse response of black box I am trying to find out the transfer function of a real life continuous-time black box.

First thought is of course to input a delta and get the impulse response as the books have taught us but I think it is not feasible as you can't have a proper delta (infinite at one point).

Q1. How can I simulate a Dirac input in a way that I don't break my system and in a way that it is going to give me a pretty accurate result? Accurate meaning something that will look a lot like the impulse response.

Second thought is to do a frequency response analysis where I will input a few sinusoids at different frequencies and store the output.

Q2. How many periods should I gather for each sinusoidal so that I get an accurate output?

Q3. What effect would a few leading zeros on my input have on my output spectrum? Is it going to be the same as zero-padding?

A3. Yes, because zero-padding the input will lead to zero-padding the output which I will then FFT.

Q4. If the input signal stops abruptly and went to zero after a sufficient signal length time, would that discontinuity affect my output spectrum as well?

If the system is continuous-time then you can't simulate an impulse (Dirac delta).

Generally people take a number of approaches:

Instead of using an impulse, find the step response of the system of interest.  This can be related to the integral of the impulse response as the unit step is the integral of the unit impulse.
Generate a set of spaced sinusoids and measure the output phase and amplitude of each.  The transfer function is found by interpolating these discrete frequency points.
Generate a swept sine wave (chirp) signal over a known frequency range and deconvolve the output to find the transfer function.
Generate a pseudorandom noise sequence, measure the response, and use it to estimate the transfer function.


What principles to use, in order to create field that will be read by OCR? I currently have an invoice, where some numbers are written by hand, and I would need to read them automatically.

My past experience is that people who write numbers will them write in so many different ways (overlapping with field name etc.) that it's really hard to isolate only digits.

Now I would like to redesign invoice, for example, to add field of some rectangular shape, that will be easily detectable (think easy to detect like a QR code is).

I'm not sure:

should I add some ordinary box, or are there any specific designs of fields/boxes/rectangles that would make finding the field and isolating everything in it (from image) simple and 99% correct?

Tech stack that would be used is:

OpenCV
Python
HTML/CSS for creating the design of invoice, that would be exported to pdf for users to fill it with the needed information

There's been forms like that forever: In Germany, the Überweisung (wire transfer) / Verrechnungsscheck (Cheque for accounting) have always been on things like this:

Link

and as far as I can tell, that works pretty well with 1990's tech.

Note how the boxes are part of the "background" and are color-wise probably designed to simply disappear in a first processing step after scanning.

However, while that might be have been at a different level of initial investment, postal services find and read postal codes / ZIP numbers that are handwritten since at least the mid-nineties, so: Yeah, bummer that people aren't good OCR-compatible writers, but maybe your OCR just has to get better: Knowing that you're only looking for numbers will help a lot. Don't take an OpenCV example as a complete OCR solution. There's more to OCR than segmenting individual glyphs and trying to find the closest match to each one from an alphabet in isolation.


SNR After Multiplying Two Noisy Signals If I have two noisy signals where the noise is not correlated and I know the SNR of each signal, how would I find the SNR after multiplying the two signals together?

The multiplication of the two noisy signal gives

$$(x_1+n_1)(x_2+n_2)=x_1x_2+x_1n_2+x_2n_1+n_1n_2=x+n\tag{1}$$

with the desired signal

$$x=x_1x_2\tag{1}$$

and the noise part

$$n=x_1n_2+x_2n_1+n_1n_2\tag{2}$$

Assuming all signals are independent of each other and have zero mean, we get for the signal power

$$\sigma^2_{x}=\sigma_{x_1}^2\sigma_{x_2}^2\tag{3}$$

and for the noise power

$$\sigma_n^2=\sigma_{x_1}^2\sigma^2_{n_2}+\sigma^2_{x_2}\sigma^2_{n_1}+\sigma^2_{n_1}\sigma^2_{n_2}\tag{4}$$

For the total SNR you get

$$\text{SNR}=\frac{\sigma_x^2}{\sigma_n^2}=\frac{\sigma_{x_1}^2\sigma_{x_2}^2}{\sigma_{x_1}^2\sigma^2_{n_2}+\sigma^2_{x_2}\sigma^2_{n_1}+\sigma^2_{n_1}\sigma^2_{n_2}}\tag{5}$$

With $\text{SNR}_1=\sigma_{x_1}^2/\sigma_{n_1}^2$ and $\text{SNR}_2=\sigma_{x_2}^2/\sigma_{n_2}^2$ this can be rewritten as

$$\text{SNR}=\frac{\text{SNR}_1\text{SNR}_2}{\text{SNR}_1+\text{SNR}_2+1}\tag{6}$$


How to delay a signal by D degrees I read in a paper that a sound signal can be delayed by D degrees. How can such a thing be possible? Isn't sound signals expressed only in terms of the time axis?

Thanks in advance

Given that you have limited number of samples, you can perform a Discrete Fourier Transform and represent your sound in frequency domain (frequencies and phases). After that if you want to apply a delay of D degrees you can simply use the known Fourier transform below:

$$ F\left \{ f(t-\theta) \right \} \rightarrow \hat{f}(w) \times e^{jw\theta} $$

Assuming that your sound is a number of sinusoids added together; to apply same phase delay to all frequency components is multiply a complex number $ e^{jw\theta} $ where $\theta$ is your desired phase delay.

On the other hand if you want to apply the same time delay to all your frequency components you should calculate the appropriate phase angle shift to be applied to each frequency component. Assuming $t_{delay}$ is the delay you want to apply, you can obtain $\theta_w$ to apply with
$$ \theta_w = w\times t_{delay} = 2\pi f\times t_{delay} $$

In which, $\theta_w $ should be calculated for each frequency bin of your results of your DFT. And the reasoning behind this can be simply explained as below. Assume a signal with sinusoid components.
$$ f(t) = sin(2\pi f_1t) + sin(2\pi f_2t) $$
If I want to apply the same phase angle shift to both signals (I don't know why anyone would want but anyways):
$$ g(t) = sin(2\pi f_1t + \theta) + sin(2\pi f_2t + \theta) $$
And to apply a time delay:
$$ f(t-t_{delay}) = sin(2\pi f_1(t-t_{delay})) + sin(2\pi f_2(t-t_{delay})) $$
$$ f(t-t_{delay}) = sin(2\pi f_1t-2\pi f_1t_{delay}) + sin(2\pi f_2t-2\pi f_2t_{delay}) $$
As you see the phase angle shift $\theta$ you should apply is different for each frequency component.

edit: If you have the proper MATLAB license you could try checking out the doc and source of delayseq function (it's in Phased Array System Toolbox).


Adding noise to frequency response I have the dynamics of a 2nd order system, mass-spring-damper for example, in the transfer-function format.
For the analysis that I am doing, I am calculating the frequency response of the system by replacing the $s$ with $j\omega$ for a given range of $\omega$ of interest.
Now, I need to examine in the behavior of the system when noise is added. How does one go about adding noise to the frequency response? Is it valid to generate some random noise (constant mean and variance) divided into a real and imaginary part and add it to the frequency response?

I don't think it's a good idea to add random noise to frequency response. In general, noise is added to the input or output signal. You should be clear that at which stage the noise is introduced. According to the different stages of noise introduction, there are two main transfer function estimators, $H_1$ estimator and $H_2$ estimator.

$H_1$ assumes that there is no noise on the input and consequently that all the input measurements are accurate. All noise is assumed to be on the output.

$$
\hat{Y}(\omega) = H(\omega) X(\omega) + N(\omega)
$$
where $\hat{Y}(\omega)$ is the observed output signal, which is polluted by additive noise $N(\omega)$. And you will see that there is a deviation between the estimated frequency response and the real one.
$$
\hat{H}(\omega) = \frac{\hat{Y}(\omega)}{X(\omega)} = H(\omega) + \frac{N(\omega)}{X(\omega)}
$$
This estimator tends to give an underestimate of the frequency response if there is noise on the input.

$H_2$ estimator assumes that there is no noise on the output. Noise is assumed to be only on input.

$$
Y(\omega) = H(\omega) \Big[X(\omega) + N(\omega)\Big]
$$
This estimator tends to give an overestimate of the frequency response if there is noise on the output.
Check here Transfer function estimate