Exact definition of soft bits in digital receiver At which stage of digital receiver, the bits are considered as soft bits?

I am unable to find the exact definition. I have rough idea how soft bits are different from hard bits but in terms of signal processing, I am confused. The signal after demodulation, phase/frequency/timing offset correction before the hard bits should be considered as soft bits?

I believe adding the samples amplitude of a single bit and averaging it over the number of samples to normalize the amplitude of single bit can be referred as soft bit right before converting it into hard bit?

Soft bits or more generally speaking, soft metrics, are usually taken at the output of the demodulator (e.g. the output of a matched filter sampled at $T_b$). A soft bit can be therefore a real value, but it is usually quantified (which is required if a DSP comes after).

A hard bit is obtained by performing a hard decision on a soft metric. So for instance, if $x[kT_b]$ is the output of the matched filter when a BPSK signal is applied, then a hard-bit can be obtained by simply comparing $x[kTb]>0$, which is $0$ or $1$.

This figure from wikipedia is self-explanatory.

As you see, the output of the matched filter is actually a continuous signal. Some part of your receiver (which is not included in the figure) must perform also timing estimation in order to determine at which instant the output of the matched filter must be sampled. Those samples are the red points in the figure, and might be considered as soft bits. So the higher the magnitude of those samples, the more information they provide about what the original transmitted symbol was. That is the advantage of having soft-bits: they don't only tell you what symbol has been probably sent, they also give an idea about that probability.


Removing an unwanted sinusoidal signal using FIR filter I am learning about the $\mathcal Z$-transform and FIR filters and I do have problem with the following exercise:

  There is given signal $$x[n] = s[n]+\sin\left(2\pi f_n n\right)$$ where $s[n]$ is a useful signal, while sinus is an interference. Let's assume, that we use FIR filter to remove the interference. Where should be placed zeros of a transfer function to completely remove sinusoidal signal?

  Sampling frequency: $f_s=10\textrm{ kHz}$

  Sinusoid's frequency: $f=2.5 \textrm{ kHz}$

Here is my attempt:

$$f_n=\frac{f}{f_s}=\frac{1}{4}$$

\begin{align}
y[n]&=x[n]\star h[n]\\
&=s[n]\star h[n]+\sin\left(2\pi f_n n\right)\star h[n]\\
\implies Y(z)&=S(z) \cdot H(z)+ \underbrace{\mathcal{Z}\left\{ \sin\left(2\pi f_n n\right)\right\} \cdot H(z)}_0
\end{align}

thus

$$0=H(z)\cdot\frac{z\sin\left(2\pi f_n\right)}{z^2-2z\cos\left(2\pi f_n\right)+1}$$

In general FIR filter is given by equation

$$H(z)=b_0+b_1z^{-1}+b_2z^{-2}+\dots$$

hence

$$\left(b_0+b_1z^{-1}+b_2z^{-2}+\dots\right)\cdot\frac{z}{z^2+1}=0$$

And I am stuck here. What to do next?

A very straightforward-answer, without a lot of calculation:

The magnitude of the frequency response of a filter is the magnitude of its z-Transform on the unit circle. Hence, if your filter should filter out a specific frequency, it needs to have a zero at exactly this frequency in the unit circle.

Considering that $f_s=10kHz$, a full rotation on the unit circle relates to $10kHz$. A sine has frequency components at $f_n$ and $-f_n=f_s-f_n$ (considering the periodicity of the DTFT due to the sampling). Hence, there need to be zeros at 2.5kHz and 7.5kHz, corresponding to 90 and 270 degree on the unit circle. Accordingly, the zeros need to be at $z=\pm j$.


Generating sound of particular frequency matlab I want to generate a high frequency (single) tone. I am using matlab (any other platform/tool for the same will also work). I am using the following snippet which I wrote by referring to various sources achieving similar specification.

amp = 10;
Fs = 44100;
duration = 60;
freq = 20000;
t = (0:1/Fs:duration-1/Fs);
%t = linspace(0,duration,Fs*duration);

wave = amp*sin(2*pi*freq*t);
%sound(wave, Fs);

audiowrite('20KHz-pilot.wav',wave,Fs);

draw_spec = 0;
if(draw_spec == 1)
    NFFT = 2048;
    wave = wave(1:NFFT);
    y = fft(wave);
    y = abs(y/NFFT);
    y = y(1:NFFT/2+1);

    % as 0 and Fs don't have any alias while others do
    % we need to double the power of the other frequencies since we are
    % converting from 2 sided power spectrum to single sided spectrum.
    y(2:end-1) = 2*y(2:end-1);

    f = Fs*(0:(NFFT/2))/NFFT;
    plot(f,y);
end

The problem I am facing is that when I load the file in audacity it shows frequencies other than 20KHz in fact 20Khz is not there in spectrogram. Further when I play the sound I can hear it which should not be the case (most adults have audio hearing range upto 18Khz). Also plotting spectrum of the wave in matlab only shows 20KHz.

I have tried to use common conventional variable names so the code should be pretty self explanatory.

I am not able to find where I am making a mistake, is there a mistake in the code or is it that due to sampling rate etc. there are other frequencies also generated (aliasing or similar concepts). I tried playing the sound on different speaker systems to check if it could be a hardware issue but that was not the case.

NOTE: I have browsed through multiple resources with similar question title as mine, unfortunately none of them were able to help me. Thus, I resorted to posting this question.

Try wave = amp*sin(2*pi*freq/Fs*t) instead. Since Fs maps to $2 \cdot \pi$, freq maps to $\frac{2\cdot\pi\cdot {\tt freq}}{\tt Fs}$.

Could be useful for reference on Wikipedia.


Measuring signal phase difference between receivers I am working on a project which aims to locate a transmitter.

I consider doing this using several receivers and measuring the phase differnces between them.

My problem is that arduino measures once every ca $dt=10^{-6}$ seconds, and therefore, to notice a phase difference, the distance between two receivers must be at least $c\cdot dt$ which is ~300 m.

I know that there are methods of finding the phase difference in much smaller systems, but I can't figure out how to solve that problem.

So, how can I measure the phase difference without placing the receivers that far away from one another?

Preface

First of all: Arduino might really not be the tool of choice here; with a 1MS/s max, and some something-in-between-8-and-16-effective-bits ADC, there's hard limits on what you can detect. I will explain the theory of how to do better estimation than just sample-wise, but I won't address these limits here.

Furtherore, I assume you

do real-valued sampling with the built-in ADC (because that is not a simultaneous two-channel/IQ ADC), meaning that you probably have some IF receivers that convert a RF bandpass signal to a directly sample-able IF signal (or one that you can undersample)
keep true to the filtering requirement of Nyquist's sampling theorem,
use the actual received signals, and not just e.g. the time of their amplitude peak, and
have knowledge of DSP basics and some understanding of the DFT.

Confirmation of your Calculation

So, your 300m calculation is correct if you can only compare two sample streams only by shifting one of them integer amounts of samples. And with this, I assume you compare using a cross-correlation function (ie. you find the product of one reception and the other, shifted by all reasonable full-sample delays, and look for a maximum).

But you don't have to do that; you can either

interpolate in time domain and calculate the cross-correllation in time domain, and find the argmax of that
compare the two sequences in frequency domain and estimate the time difference directly.

Interpolation in Time Domain

Since we know that you've adhered to Nyquist sampling theorem (you did, didn't you?), we want to keep it that way. To properly interpolate by a factor of $N$, you just pad every sample you get with $N-1$ zeros, and then low-pass filter the result to the original band (ie. apply an $\frac 1N$th band-filter). That gives you a $N$ times higher resolution.

However, it also means that your correlation  just went in complexity from $L^2$ ($L$ being the length of your sample sequence) to $N^2L^2$, which is really unfortunate on an Arduino. So you'd have to do that on a PC, probably.

Frequency Domain Techniques

You're probably aware of the time-shift property of the Fourier or the Discrete Fourier Transform (DFT). For the latter it states that

$$
\mathcal F\{x[t-\tau]\}(f) = e^{-j\frac{2\pi}{L} \tau f} F\{x[t]\}(f)
$$

So let's name the signal that the "closer" receiversees $x_0$, and let's define the other receiver sees, $x_1$, to be $\tau$ "late"; it directly follows that
$$\begin{align}
\mathcal F\{x_1[t]\}(f) &= \mathcal F\{x_0[t-\tau]\}(f) \\
&= e^{-j\frac{2\pi}{L} \tau f} F\{x_0\}(f) \\
&\implies\\
\frac{\mathcal F\{x_1[t]\}(f)}{\mathcal F\{x_0[t]\}(f)} &=e^{-j\frac{2\pi}{L} \tau f}\\
&=: \varphi(f)
\end{align}$$

I introduced the "auxillary" function $\varphi(x)$ here – so if you look at that,

$$\varphi(x)=e^{-j\frac{2\pi}{L} \tau x}\text,$$

you'll notice that it's a complex sinusoid with frequency $\tau$!

So, by estimating the frequency of that, you get $\tau$,  your time delay. Notice that the variance of a frequency estimate reduces with the number of samples you observe – in other words, the longer your DFT (you'd probably use an FFT), the higher your frequency resolution, and the smaller your delay granularity, and the better your Time-Difference of Arrival estimate.

Remarks

Sampling Rate of $\frac1{1\text{ µs}}= 1\text{ MHz}$: That sounds too high, actually. Where did you get that number? I ask because I have my doubts: what are you going to do with that amount of data inside an Arduino? Spit it out via the serial port? Certainly not, that's much slower. Process it on-Arduino? Don't think so, unless I underestimate the processing power of an ATMega seriously. Store it in RAM and use it later? Well, at 1MS/s, how long will your RAM last?

Correlation in time Domain is seldom done directly, but often by (conjugate) multiplying the DFTs of the sequences, and then reverting back to time domaind – an FFT has a complexity of $\mathcal O(N\log N)$, and thus the whole "transform both sequences – multiply – inverse transform" has complexity of about $\mathcal O(3N\log N+ N) "<" \mathcal O(N^2)$. When you write the fast correlation down, you'll notice similarities to the frequency-domain approach above, if you apply the shift theorem.

Variance of the $\tau$-estimator can't be pushed down indefinitely with real-world signals, since well... real-world signals do have finite length.


Low pass signal filtering using FFT or simple infinite impulse response filter I am currently in need for filtering accelerometer data for an Android application.
First of all, I used a simple low-pass filter (simple infinite impulse response filter) as follows:

for i from 1 to n
   y[i] := y[i-1] + α * (x[i] - y[i-1])

This helped me achieve a smoother result.

Then I decided to play with FFTs. I used a fast-fourier transform to get the signal into frequency domain and then zeroed some of the high frequencies. Then using inverse fourier transform I recreated the signal. This all worked fine and I know that the FFT and IFT implementations are fine. However, the signal wasn't as smooth as the one that I got from before using the simple infinite impulse response filter. I tried zeroing some further frequencies but didn't give me as good of a result as expected.

What is the reason behind this? I though using FFTs and IFT should technically give me a nice smooth graphs. Is this because of the sampling in FFT?

Thanks

There are three potential problems with what you are doing:

Zeroing out the high frequencies causes ringing due to the Gibbs phenomenon.
You are essentially multiplying the Fourier transform with a square function, which is equivalent to circularly convolving the time-domain data with the inverse transform of the square function, which is a sinc function.  Unless you zero-pad the data before transforming it, the circular convolution will cause the data at the end of the FFT block to affect the data at the beginning, and vice versa.
If you were doing multiple FFT's to process a stream of data, unless you use a technique like Overlap-Add or Overlap-Save, there will be discontinuities at the boundaries between FFT blocks.  This problem is related to #2.


Color Tracking a Red Laser Height I'm making a Laser Harp.
The application needs the position of the hand along each of the laser beams to modulate a sound synthesis or modify the properties of a sampled instrument.
I'm wondering if doing a Color Tracking with OpenCV will do.

Make the beam visible with fog machine.
Segmenting the image for each beam.
Measure the height of the beam.
PROFIT.

Other way can be using an infra-red light sensor like Infusion Systems - ReachClose v3.1
What method is more accurate and fail-proof?
Any other ideas?
Diagram3

I can only say something about the computer vision part.

If your doing color tracking you have to take different lighting conditions into consideration, although it will be easier with a laser (compared to other things) since it at least have the same color.

If you are to use this device professionally you should have a delay about 20ms else it will be noticeable. This forces you to use a high speed camera of at least 100fps (10ms per frame) to have time for processing and sending the sound through the amplifiers.

Since you are only tracking straight lines I would recommend to first filter out everything but the laser-color and then hough transform it to find all straight laser lines. And as a final step implement an ad-hoc method to find where the laser ends shouldn't be that hard.


Practical vs ideal lowpass interpolator Consider a signal with a sample rate $f_s = 44.1$ kHz. Let us upsample the signal by a factor of $L = 2$ and interpolate the zeros.

An ideal lowpass interpolator would have a gain of $L$ and a cutoff frequency of:

$$f_c = \frac{f_s}{L}$$

An ideal lowpass filter has an infinitesimally small transition band.

In practice I see real lowpass interpolators have a small transition band centred around $f_c$.

The transition band can be quite large, say, $0.45 f_s$ to $0.55 f_s$.

My question is: why do we centre the transition band of a practical lowpass interpolator around the ideal cutoff frequency? By doing that the practical lowpass stopband is above the ideal cutoff which does not make sense to me as that will allow a small unwanted spectral image from the $0.45 f_s$ to $0.50 f_s$ region to creep into the new signal. The obvious alternative is to make the stopband of the practical lowpass $0.5 f_s$ and put up with a passband starting at $0.4 f_s$ assuming we can't make the transition band steeper. There must be some reason this isn't the way it's done.

When designing a filter, you really care about its behavior in two regions:

Passband: You want little attenuation in this region, and maybe other properties as well, like linear phase, depending upon your application.
Stopband: You want as much attenuation as needed in this region.

Between these two is the transition region. This is treated as somewhat of a "don't-care" band. You don't typically constrain the response too tightly in this area, so you can't really count it being usable. In your example, the passband lies below $0.45 f_s$; after the interpolator, you only plan on using frequency content below this threshold.

This means that you can allow some aliasing in order to simplify your filter design. Your transition region starts at $0.45 f_s$; everything above that frequency in your filter's output will either have a response that is unpredictable given your filter specs (if it lies in the transition band) or one that is highly attenuated (if it lies in the stopband). The takeaway:

  You can't rely upon any frequency content past the passband edge in your filter output anyway. So if it makes the filter design cheaper, why not allow the frequencies above the passband edge in the filter output to contain aliased garbage?

This technique is commonly used in multirate filters as you've noticed, as it allows savings in the required filter order in order to meet a given set of passband/stopband specifications.


Why is it called 'multi-level thresholding', if only the best threshold is selected? I'm referring to Otsu's thresholding, and saw this from his paper:

When he says he's trying to maximize the between class variance of thresholds k* 1 and 2, is it the threshold with maximum value that is chosen? Is that what the "max" means? If it was just about selecting one threshold, why is it called multi-level thresholding?

On the other hand, if multiple thresholds are chosen, then why do we need to find the "max" of it?

ps: I'm trying to implement Otsu's method in a computational intelligence algorithm (Differential Evolution and Particle Swarm Optimization), and not a single research paper has explained clearly, how the thresholds are to be used for multi-thresholding. Not even the original Otsu paper.
Also, after implementing the algorithm, the image with the best fitness threshold is a different (much darker) threshold than what Matlab's greythresh threshold returns.

Using k1 and k2 you can calculate a between class variance value using the Otsu formula. However, there may be many combinations of k1 and k2 for which you have to calculate between class variance values to see which of them gives the highest between class variance value. That's what 'max' means. So whichever k1 and k2 gives you the highest between class variance, you choose those values of k1 and k2 as your thresholds. Same can be done if you want more thresholds.


Dimensional analysis in the Morlet wavelet I'm a complete beginner on signal processing, so please forgive the basicness of this question - but I'm very confused by this. Regarding the Morlet wavelet, Wikipedia says:

  The wavelet is defined as a constant $\kappa_\sigma$ subtracted from a plane wave and then localised by a Gaussian window:
  $$ \Psi_\sigma(t) = c_\sigma\pi^{-\frac{1}{4}}e^{-\frac{1}{2}t^2}(e^{i\sigma t} - \kappa_\sigma) $$
  where $\kappa_\sigma=e^{-\frac{1}{2}\sigma^2}$ is defined by the admissibility criterion, and the normalisation constant $c_\sigma$ is:
  $$ c_\sigma = \left( 1 + e^{-\sigma^2} - 2e^{-\frac{3}{4}\sigma^2} \right)^{-\frac{1}{2}} $$
  ......... The parameter $\sigma$ in the Morlet wavelet allows trade between time and frequency resolutions. Conventionally, the restriction $\sigma>5$ is used to avoid problems with the Morlet wavelet at low $\sigma$ (high temporal resolution).

Presumably, given the comment about the typical restriction of $\sigma>5$, this trade-off parameter $\sigma$ is meant to be a dimensionless value. But then (or, in fact, whatever the dimension of $\sigma$ is), the formula for $\Psi_\sigma(t)$ is not dimensionally invariant under changes of the units in which time is measured! Shouldn't there be at least some extra parameter in front of $t$, so that the formula makes sense in terms of dimensional analysis, if it is meant to be applicable to real-world data (as suggested by comments about "conventional restriction" regarding trade-off)?

(Incidentally, I've seen elsewhere much simpler formulae for the Morlet wavelet, which do make dimensional sense, and where $\sigma$ has the dimension of time - which makes sense, as low $\sigma$ means a narrow time-window of significant relevance.)

Okay, I think I've understood it now: The number $t$ in $\Psi_\sigma(t)$ is a dimensionless number, because in the formula for the wavelet transform it ends up being the ratio between an actual time and a time-scale.


Applying frequency-domain filters on a centered Fourier transform I understand why we shift the Fourier transform such that the 0-frequency is centered for visualization. In the shifted DFT(u,v) of an M*N 2-dimensional image,

the top-left corner of the 4th quadrant is (0,0) frequency or (low u, low v)
the bottom-left corner of the 1st quadrant, (M-1,0) or (high u, low v);
the bottom-right corner of the 2nd quadrant, (M-1,N-1) or (high u, high v); and
the top-right corner of the 3rd quadrant, (0,N-1) or (low u, high v).

Now, when we apply filters (centered again) on the centered DFT, aren't we changing both low and high frequencies closer/further to/from the center?

For example, the centered Gaussian high-pass filter is centered at M/2 and N/2, and supposed to attenuate only low frequencies (is it?). However, applying this filter to the shifted DFT will attenuate not only low frequencies, but also high frequencies in the 1st, 2nd, and 3rd quadrants.

I did a little experiment and can confirm this effect. Applying the centered Gaussian high-pass filter to a centered DFT was not equivalent to applying the non-centered Gaussian high-pass filter to the same non-centered DFT. I had to apply the filter three more times at each corner of the non-centered DFT to get the same result.

I couldn't find any good explanations why this (high-pass filter changing high frequencies around the center or low-pass filter changing low frequencies) is acceptable.

To simplify this question, let's take a one-dimensional image as example. The non-shifted DFT has the 0-frequency on the left edge while the M-1 frequency on the right edge. The shifted DFT has the 0-frequency at the center and the M-1 frequency right next to it on the left side. Attenuating this center in effect changes both the 0 and M-1 frequencies. This is different from attenuating the 0-frequency on the non-shifted DFT.

Thanks in advance!

Note that an FFT of strictly real data is conjugate symmetric.

For a length M FFT of strictly real data, the data in FFT result bin M-1 has the same magnitude as in bin 1 (for a low frequency), but complex conjugated.  So to maintain symmetry after filtering, which is necessary for the IFFT-ed image to remain strictly real, a filter has to modify bin(M-i) by the same ratio as bin(i).


Any Open Source Fast Wavelet transform libraries? I am in need of an open source library for computing Fast wavelet transforms (FWT) and Inverse fast wavelet transforms (IFWT) - this is to be part of a bigger code I am currently writing.

The things I am looking for in the library:

1) Contains a good variety of wavelet families (Daub,Haar, Coif etc.)

2) Ability to run in parallel - VERY IMPORTANT

3) Reasonable documentation, so that I can include new wavelets of my own without having to modify the complete source ( Maybe an OOPS based approach would help?)

I am flexible about the language - C/C++/Fortran90/Fortran77/Python .... Any language would do for me, though I would prefer something which is optimized for speed and in parallel.

So far, I have found PyWavelets - it looks good, but it is in Python ( therefore considerably slow) and it doesn't run in parallel. I am going to use the DWT for processing really huge datasets, so speed is an important concern for me.

Its understandable there may not be anything with all the requirements I mentioned. But I wanted to hear from the community if they have any suggestions.

you can have a look at the LTFAT's wavelet module
http://ltfat.sourceforge.net/doc/wavelets/index.php

it runs in Matlab/Octave with backend written in C. It has fairly large database of wavelet filters and new ones can be added easily.

What exactly do you mean by

  2) Ability to run in parallel - VERY IMPORTANT

Should the computation itself be somewhat parallelized or that it should be possible to run several batches of FWT calculations? I guess that running several instances of Octave should do the trick in the latter case.


Inverse DTFT of left-sided sequence I am pretty new to inverse Fourier transforms and I would like to ask a question. Does anyone know how to bring back to the "sequence domain" this relation below? In other words, get the inverse DTFT of this expression? (Given that C, D and $a$ are constants and you know them):

$$ \frac{C \cdot e^{-j \omega}+D}{(1 - a \cdot e^{ j \omega })^2} $$

There seems to be a mistake in your derivations, the following code shows the correct one for $|\alpha| < 1$ ,  as indicated by @Matt.L $:

N = 128;         % simulation length (choose large enough)
n = -N : 2;      % simulation time span

C = 1.5;          % set some values
D = 2.78;
a = 0.33;         % a is for alpha

% x = -(n-2).*C.*(a.^(-n+1)).*stepfun(-n+2,0) + ...
%    -(n-1).*D.*(a.^(-n)).*stepfun(-n+1,0);

x = zeros(1,130);  % use the explicit code as stepfun() has problem with n.
for n = -N:1
    x(n+N+1) = -C*(n-2)*(a^(-n+1)) - D*(n-1)*(a^(-n));
end

M=2000;          % Select a FFT length for demonstration purposes
F=fft(x,M);      % compute the DFT of x[n] so obtained

w = linspace(0,2*pi-2*pi/M,M);              % choose DFT frequency samples
H = (C*exp(-1j*w)+D)./((1-a*exp(1j*w)).^2); % Evaluate the given DTFT

figure,plot(w,abs(H));title('H'); % compare the magnitudes of DTFTs
figure,plot(w,abs(F));title('F'); % Phases will be different as this
                                  % is an acausal sequence.

Derivation of the result for the left-sided sequence  $|\alpha|< 1$ case) is as follows:

Intuition plays a central for getting the answer as clean as possible which involves inverse Fourier transforms. In this example we shall use those three fundamental Fourier transform properties.

1- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=x[-n] \leftrightarrow Y(\omega) = X(-\omega) $

2- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=nx[n] \leftrightarrow Y(\omega) = j \frac{d X(\omega)}{d\omega} $

3- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=x[n-d] \leftrightarrow Y(\omega) = e^{-j \omega d} X(\omega)$

Begin by $$x[n] = a^n u[n]  \longleftrightarrow X(\omega) = \frac{1}{1 - a e^{-j \omega} }$$

Let $y[n]=x[-n]=a^{-n} u[-n]$ then $Y(\omega)= \frac{1}{1-a e^{j \omega}}$

Apply second rule:
$z[n] = n y[n] = n a^{-n} u[-n]$ then $Z(\omega) = \frac{-a e^{j \omega}}{(1-a e^{j \omega})^2}$

Apply linearity:
$  w[n] =  -a^{-1} z[n] = - n a^{-1} a^{-n} u[-n]$ then $ W(\omega) = \frac{e^{j\omega}}{(1-a e^{j \omega})^2}$

Apply shift property on $w[n]$ to get rid of the numerator exponential.
$v[n] = w[n-1] = - (n-1) a^{-1} a^{-(n-1)} u[-(n-1)]$ then $ V(\omega) =\frac{1}{(1-a e^{j \omega})^2}$

Now apply shift property and linearity to $v[n]$ such that the resulting DTFT will be $$H(\omega) = \frac{ C \cdot e^{-j \omega} + D} {(1-a e^{j \omega})^2} = \frac{ C \cdot e^{-j \omega}} {(1-a e^{j \omega})^2} + \frac{ D} {(1-a e^{j \omega})^2}$$

which implies that
$$h[n] =  C v[n-1] + D v[n] $$ hence we get $h[n]$ as:

$$h[n] = C \left( - ((n-1)-1) a^{-1} a^{-((n-1)-1)} u[-((n-1)-1)] \right) - D \left( (n-1) a^{-1} a^{-(n-1)} u[-(n-1)] \right) $$

simplify to get the result:
$$ h[n] = C (-n+2) a^{-n+1} u[-n+2] + D (-n+1) a^{-n} u[-n+1] $$


Analysis of a LTI system using DFT Consider an LTI system $$H(z)=1-\frac{1}{2}z^{-1}+\frac{3}{4}z^{-2}$$
  Let $x[n]=(\frac{1}{3})^n\cdot u[n]$ be the input signal. It is desired to determine the output for $n=0,1...,N_a$. To achieve that, a relevant part of the DFT will be used together with an uniform sampling of $H(z)$ in the unity circle with a separation of $\frac{2\pi}{N_b}$.
  a) Explain the process to achieve this. b) Are there any restrictions for $N_a$ and $N_b$? c) For the same $x[n]$, would the method described in a) work if $H(z)=\frac{1}{1-\frac{1}{2}z^{-1}}$? If not, correct the previous method to find the correct result.

$H(z)$ is sampled in $N_b$ points. $\hat{h}[n]$ (the sequence whose DFT is the sampled version of $H(z)$) has only 3 non-zero values, so... the sampled version of it would have $N_b -3$ zeros (I'm not sure about this). About the size of $X[k]$, I thought that, if I want to know the first $N_a +1$ values of $y[n]$, maybe it is enough to grab the first $N_a +1$ points of $x[n]$ to do the DFT and add the necessary amount of zeros to avoid aliasing? This may be false but I just don't know how to determine the relevant part of the input signal that has to be used. If everything I said above is correct (I highly doubt it), then I thought that the only restriction on $N_a$ and $N_b$ is that $N_a +3 = N_b$.

About c), I have no idea of how to do that.

Any suggestion?
Thanks for your time!

Since the system described by $H(z)$ is causal, it's sufficient to consider the input sequence in the range $0\le n\le N_a$ in order to compute the output sequence in the same range. Each output sample can only depend on the current input sample and on past input samples. Since the length of the truncated input sequence is $N_a+1$, and the length of the system's impulse response is $3$, the length of the convolution of these two sequences is $N_a+3$. If you want to compute this linear convolution using circular convolution (i.e., using the DFT), the DFT length $N_b$ must satisfy $N_b\ge N_a+3$. So you must zero-pad the truncated input sequence as well as the impulse response to that length, then compute the DFTs of both zero-padded sequences, multiply them, and compute the IDFT of the resulting sequence, the first $N_a+1$ elements of which equal the desired output sequence. This answers a) and b).

For part c) it will not suffice to sample $H(z)$ on the unit circle, because the corresponding impulse response is infinitely long. However, for computing the first $N_a+1$ output samples it is sufficient to consider the first $N_a+1$ samples of the impulse response. Just like before, this works because the system's impulse response as well as the input signal are both zero for $n<0$. So you need to consider a system $\hat{H}(z)$ the impulse response of which is a truncated version of the infinite impulse response corresponding to $H(z)$. The length of the convolution of the truncated input signal and the truncated impulse response is $2N_a+1$, so the DFT length $N_b$ needs to satisfy $N_b\ge 2N_a+1$. The rest of the procedure is the same as before.


A Delay Between Two Filtered Chaotic Signals it is a common practice to use a shift of cross-correlation peak to evaluate a group time delay of two signals (chaotic signals are included). Can one synchronously and equally filter these chaotic signals (one is a shifted copy of the other) and then get a cross-correlation between filtered copies of them? What will happen with a group time delay after that filtration? I made a numerical experiment with two simple delta-correlated signals (one is a shifted copy of the other, fft, fft filter and inverse FFT to get a cross-correlation) which indicated that absolutely nothing happens with a delay time after filtration. The same result one can get for correlated signals (I applied a moving average to simple delta-correlated signals and again - nothing happens) The delay time is the same after filtration (with different bandpasses of filters). But when I come to real signals the picture had changed dramatically. The value of the delay time depends on the filter's cut-off limits and the size of the spectral window. Any ideas why? Why do we suppose that we evaluate a delay time correctly when all spectral components are taken into account and can not be sure in the same estimations after filtering?

I assume that your experiment where you saw no effect was with signals that underwent a linear phase vs frequency process in their delay (constant group delay). Under many situations the phase of the channel in the frequency domain is non-linear, causing different frequency components to experience different delays (group delay variation). When you frequency select portions of such signals, you will then see the result of the delay that is dominant for those bands.

This is the same reason it is not best practice to simply use cross-correlation peak to evaluate the group delay. See further details specific to the challenges of doing this, and the better solution for those cases, at this posting here:

How determine the delay in my signal practically

And here:

How to reliably compute the group delay of a comb filter

At the bottom of the first post is a plot of group delay that illustrates my first point. This is an extreme case of delay variation where we see as duplicated below there are frequency selective regions where the delay is as much as 15 samples, while generally elsewhere it is only around 2-3 samples. If we further filtered this signal (with a linear phase filter to not cause further delay variation), we change the weighting between a 15 sample delay and a 2-3 sample delay.


Sinusoidal liftering in implementations of MFCC Some implementations of MFCC apply sinusoidal liftering as the final step in calculations of MFCC. It is claimed that speech recognition can be significantly improved. For instance, if $\text{MFCC}_i$ is a cepstral coefficient, and $w$ is a lifter, then
$$\widehat{\text{MFCC}_i}=w_i\text{MFCC}_i$$

is a liftered cepstral coefficient, where $w_i$ for sinusoidal liftering is defined as:

$$w_i=1+\frac{D}{2}\sin\Big(\frac{\pi i}{D}\Big)$$

When I look at the equation, I understand the sinusoidal function has a shape such that its maximum is in the middle and approaches to zero at edges. Therefore, cepstral vector's first and the last coefficients are is reduced to zero while the middle one is intact.

Why is liftering applied and how does it improve the speech recognition?

This paper$^1$ suggests:

that the reason is because the filtered versions give less weight to the higher coefficients which provide less discrimination than the lower ones.

That paper references this one$^2$ which shows the plot below indicating how the application of window reduces the variability.

References

K.K. Paliwal "DECORRELATED AND LIFTERED FILTER-BANK ENERGIES FOR ROBUST
SPEECH RECOGNITION"
B.H. Juang, L.R. Rabiner and J.G. Wilpon, "On the use of bandpass liftering in speech recognition", IEEE Trans. Acoust., Speech and Signal Processing Vol. ASSP-35, No. 7, pp. 947-954, July 1987.


Can I square signal in the frequency domain to reduce noise? Let's say this is the FFT of my measurements:

You can see that there is a clear peak and the rest is noise. I want to reduce noise, so I square the FFT and normalize it:

The peaks are now more clear, so the (relative) noise has been reduced, at the cost of changing the non-noise signal too. It's now easier to identify the frequency peaks, and of course we can use powers larger than 2 to increase this effect.

Is there a trickery? Is this commonly used? I guess the main problem is that the signal is affected.

Squaring your signal is multiplying it with itself.

Multiplication in frequency domain is convolution in time domain.

So you're correlating the signal in time domain with its time-inverse. I don't think this has beneficial effects aside from a very specific class of signals (namely, the signals which are their own time-inverse).

You forget why you're interested in SNR in the first place. It's not for the sake of SNR itself. It comes from the assumption that your signal is essentially unharmed but overlaid with additive noise. You're breaking that assumption.

Also, the effect is purely optical; you know, squaring positive numbers is a monotonous operation. So if you need to find something "clearly larger" than something else, you can do this with the original graph just as well.

Adapting the Y-axis to be able to spot signals of very different magnitudes is a very common thing to do, however: for PSD-style plots you typically don't look at abs(spectrum) in a linear scale, but convert it to dB, which is a logarithmic scale; in fact, it hence does quite the opposite of your polynomial scale, but it's much more useful for real-world signals.

There's a few cases where you'd square the time signal (e.g. to allow for easy timing recovery of BPSK signals), but from the top of my head, I can't think of anything where you'd work with the squared DFT; might be useful if we know that a signal has a time-inverse part and were looking for that.

Note that Laurent correctly pointed out that "squaring" is a bit of an imprecise term, and in the complex sense you'd normally use it as "multiplying with the complex conjugate", but that doesn't make much of a difference here.


Image segmentation algorithm I wonder if there are algorithms used to segment images efficiently. The ones I tested by now can not segment hair correctly as on this image:

I just want to know if someone can give me some efficient algorithms to segment the hair on this image for example without cutting it and so on (correct segmentation)

Note that I tested famous algorithms such as watershed, grabcut, mean shift ... but no one of them segments this hair correctly.

Thank you very much in advance.

You can find what you need here where a comparaison of lot of matting algorithms is done. You can see the rank of those algorithms' performances as well as online demonstration on a natural and complicated image.

From that study, you can guess that A Global Sampling Method for Alpha Matting is the best algorithm to resolve your problem. The algorithm is conceived by A Global Sampling Method for Alpha Matting
Kaiming He, Christoph Rhemann, Carsten Rother, Xiaoou Tang and Jian Su.

Hope this helps.

Begueradj.


How Does a Separable Filter Work? I am new to signal processing, and was reading about filters. I came upon the sobel filter which has the following forms,

From what I understand to apply this filter to a 2-d signal, place the filter's center on the correct point and multiply the coefficients with the signal at the points of overlap and add them together. So for the following example signal, applying the horizontal sobel filter from above would give:

Signal                 Output (edges are left alone)
|100|100|100|100|      |100| 100| 100|100|
| 10| 10| 10| 10|      | 10| 360| 360| 10|
| 10| 10| 10| 10|      | 10|-360|-360| 10|
|100|100|100|100|      |100| 100| 100|100|

Separable filters are more computationally efficient and the sobel filter is one of them. My problem is I can't figure out how the separable filter is applied. So, what I want to ask is how would the sobel filter be applied separably?

Steve Eddins has a nice page that explains it.

It comes down to being able to do a separate 1-D filters in each direction (requiring $MN(P + Q)$ multiplies and adds), rather than a single 2-D convolution ($MNPQ$ multiplies and adds).  Here the image is assumed to be $M \times N$ and the convolution kernel $P \times Q$.

Of course, this assumes you're doing everything in the time domain, and not using FFTs.


What are the differences in overtones for hard vs soft materials? I've been researching Physical Modeling lately, and so far pretty much everything makes sense from a physical/logical perspective.

The one thing that I haven't yet found a proper explanation on is the difference in overtones between soft and hard materials.  I've seen it mentioned in a few places that this stiffness is the parameter which defines the difference between metallic and nylon string on a guitar, for instance.  I've also caught that one of the two directions (hard vs soft) "spreads the harmonics."  But that's all I've got.

So, what are these differences in overtones for rigid materials?  They seem to be more spread out, I've come to recognize metallic sounds as being rather comb like.  But I'm trying to eventually write code to synthesize these sounds, and need to know some details as to the way they are spread.

Also, here's a great video in which most other aspects of physical modeling are pretty well explained, as well as the mention of more spread out overtones for harder materials (though I don't remember at what time): https://www.youtube.com/watch?v=dUcNzPhZdwk

I am not sure if this is exactly answering your question but it may help you in your modeling efforts.  The metallic strings are going to be fairly simple to model in that the differential equations you use will be only dependent on the shape of string at any given time and the tensile strength of the string.

The nylon strings exhibit non-linear viscoelastic behavior which has a time dependent component to it.  That means that is has a memory and the local velocity of the string will be important.  By memory I mean there is residual strain built in once plucked and the time for recovery will be much greater than the frequency of motion.  I believe this is what gives it a “softer” sound.  Classically viscoelastic materials are modeled with a spring and damper in parallel.

Descretizing the string over a number of elements and using a divided differences approach to approximate the differential equations you should be able to get a theoretical sense of how these differences affect the observed signal.

Hope that is helpful for you in thinking about your problem.


Pointers for digital FM synthesis I am working on a digital synthesis project, in python and c currently. I have made a Sine and Saw wave generator and have now been working on FM between the two wave types. I have successfully made FM work with a Sine wave as both carrier and modulator, but whenever I add a Saw wave into the mix, The carrier frequency increase indefinitely. After a lot of research and trying to understand some more complex DSP concepts (instantaneous phase and frequency), I am feeling a bit lost, and hoping for some pointers to resources where I can learn more at a more basic level of understanding.
Essentially, the goal right now is to modulate the frequency of any wave shape with any other wave shape given only the amplitude of the modulator signal. This would in essence be a VCO emulation or a pitch envelope capable of audio rate modulation.
The code below is currently what I have come up with, which works, but only as intended when the modulator frequency is greater than the carrier frequency (I am assuming because this is phase modulation, not FM).

# self.count == the current sample
# self.fs == 48000
# fmcv == the amplitude value for the given sample number, between -1 and 1, of the modulator

pos = self.count / self.fs
value = math.sin(self.freq * 2 * math.pi * pos + math.pi * 2 * fmcv)
value = value * self.gain + self.offset
self.count = self.count + 1
return value

I have also tried the following which gives similar results.
pos = self.count / self.fs
value = math.sin((self.freq * (1 + (0.9 * fmcv))) * 2 * math.pi * pos)
value = value * self.gain + self.offset
self.count = self.count + 1
return value

Any help and pointers towards good resources where I can understand the concepts better are greatly appreciated.

If the case with two sines is clear (sine wave modulation of a sine wave carrier) then consider expanding any other shape (such as the ramp) into its Fourier components and then superimposing the result of each.
Any waveform with sharp transitions (such as a ramp) requires an infinite number of frequency components, so better to work with smoothed waveforms where the higher frequency components are essentially filtered out.
Another tip when simulating carrier modulation is to not actually simulate the carrier but use a complex baseband signal instead (what you simulate their will match what would occur at any carrier frequency and doesn't require as high of a sampling rate to simulate). This requires understanding the baseband equivalent analytic signal, but would be the recommended next step regarding understanding DSP concepts.


How is multiplexing achieved in spread-spectrum modulations like CSS? I understand why it doesn't interfere with narrowband communication systems but how do different clients communicate in a system like LoRa. Is it time orthogonality like TDM?

disclaimer

Because of the lack of document about LoRa CSS in Internet, the analysis below may be wrong with respect to LoRa system, not to the (linear) CSS principle. Any comment or update are appriciated. Thanks.

LoRa CSS modulation

LoRa CSS uses linear chirp spread spectrum. To transmit $\log_2 M$ bits, LoRa CSS modulation divides bandwidth $[-\frac{B}{2} \dots \frac{B}{2}]$ to $M$ portions.

In the figure, to transmit symbol $-k$ with $-M/2 < k < M/2$, LoRa use the waveform having the instantaneous frequency:
$$
  f(t)=\begin{cases}
               \mu\frac{B}{T}t - k\frac{B}{M} + B \qquad -\frac{T}{2} \leq t \leq -\frac{T}{2} + k\frac{T}{M}\\
               \mu\frac{B}{T}t - k\frac{B}{M} \quad\qquad -\frac{T}{2} + k\frac{T}{M} < t \leq \frac{T}{2}\\
               0 \qquad\qquad \textrm{otherwise}
            \end{cases}
$$

with $\mu = +1$ or $\mu = -1$. The spreading factor is $M = BT$.

By sampling at exactly $B$ Hz, $f(t)$ becomes continuous, we can works directly with $f(t) = \mu\frac{B}{T}t - k\frac{B}{M}$ with $-\frac{T}{2} \leq t \leq \frac{T}{2}$ in linear CSS decoding framework. (In fact it is not simple like that because of the transition phase at the discontinuos point in the figure. However, it does not change the insight of the answer. I will update my answer "later" :). I appreciate the fact that someone interested in math can do it properly)

\begin{eqnarray}
f(t) = \mu\frac{B}{T}t - k\frac{B}{M} = \frac{1}{2\pi}\frac{d\phi(t)}{dt}\\
\phi(t)  = 2\pi\left(\mu \frac{B}{2T}t - k\frac{B}{M}\right)t
\end{eqnarray}

To decode chirp, we use inverse chirp at receiver:
\begin{equation}
c_d(t) = e^{j2\pi\left(-\mu \frac{B}{2T}t^2\right)}
\end{equation}

Decoding by multiplication in $-T/2 \leq t \leq T/2$:
\begin{equation}
z(t) = c(t)c_d(t) = e^{j\phi(t)} \times ^{j2\pi\left(-\mu \frac{B}{2T}t^2\right)} = e^{j2\pi\left(\mu \frac{B}{2T}t - k\frac{B}{M}\right)t} \times e^{j2\pi\left(-\mu \frac{B}{2T}t^2\right)} = e^{j2\pi(-kB/M)t}
\end{equation}

Remember that we have assumed working on the signal sampled at $B$ Hz ? Thus instead of having $z(t)$ what we got is a sequence of sample $z[m] = z(t=m/B)$ with $-M/2 \leq m \leq M/2$. Note that $M=BT$ if you feel lost.

Now DFT of M samples $z[m]$ to have $Z[u]$ with $-M/2 \leq u \leq M/2$, you will have exactly only one peak at the bin $u = -k$ (if no synchronization error, no Doppler, no multipath). Bingo.

Multiplexing

It depends how we define "multiplexing". If we use the term in multiple access level (MAC layer), we can have TDMA and FDMA-like. If we are talking about user communicating simultaneously, the only thing that (I think) we can use is the factor $B/T$, or the spreading factor $M = BT$ if $B$ is fixed.

Let's fix the spreading bandwidth $B$. The chirp duration $T = M/B$ changes if $M$ changes and so do the chirp wavefor that depends on $B/T$ : $c_1(t) = e^{j2\pi\left(\mu \frac{B}{2T_1}t - k_1\frac{B}{M_1}\right)t}$

Thus
\begin{equation}
z(t) = c(t)c_d(t) = e^{j2\pi\left(\mu \frac{B}{2T_1}t - k_1\frac{B}{M_1}\right)t} \times e^{j2\pi\left(-\mu \frac{B}{2T}t^2\right)} = e^{j2\pi(\mu \frac{B}{2}(\frac{1}{T_1} - \frac{1}{T})t - k_1B/M_1)t} = e^{j2\pi(\mu \frac{B^2}{2}(\frac{1}{M_1} - \frac{1}{M})t - k_1B/M_1)t}
\end{equation}

$M_1 \neq M$, $z(t)$ has instateneous frequency varying in time. The DFT result should be something like that, i.e. containing no relevant peak.

LoRa defines 6 spreading factor $2^7$ to $2^{12}$ (7 if we take the special SF $2^6$). Thus, for a fixed bandwidth, LoRa system can have simulatenously 6 (or 7) users.


Sample vs. sample values Is there any difference between the terms sample and sample value in DSP? In statistics, those are very different concepts: a sample is a subset of a population: it is a set of random variables. A sample value is the realization (aka. observed value) of this set.

In DSP, I feel that sample is often used with the same meaning as sample value:

Yes, the confusion comes from sampling theory. Thus, when you sample a signal, you get ... samples. If you want to be clear you can call the other one a statistical sample.


Cut-off frequency in Fourier-Spectrum of a noise message I have a problem involving cut-off frequency.

I have a voice message in an audio file .mat which have an eletromagnetic noise. I am trying to recover the original message by applying digital Butterworth filter with MatLab. To this, besides of design the Butterworth filter, also I have to find the cut-off frequency to eliminate all, or partially, the electromagnetic noise.

The picture below shows the Fourier Spectrum of the voice message with the electromagnetic noise. How could I find the cut-off frequency to eliminate the noise? Thanks in advance!

Sampling frequency = 44.8 kHz.

I think you consider that pulse in the spectrum as the noise. Then you should remove it by a lowpass filter of cutoff frequency
$$ f_c = 5 ~\text{kHz.}$$

Slightly less than this should be used for guarding against the transition width of the lowpass filter. So you better select something like $4.8$ kHz for the lowpass filter cutoff frequency.

Note that your sampling frequency $F_s$, probably, is $44.1$ kHz and not $44.8$ kHz. Then your digital frequency for cutoff into the Matlab functions will be:
$$f_{cd} = \frac{ 2 f_c }{F_s} $$

with the given values of $F_s$ and $F_c$ your digital cutoff will be:

$$f_{cd} = \frac{ 2 \times 5  \times 10^3 }{44.1 10^3} = 0.2268 $$

or better

$$f_{cd} = \frac{ 2 \times  4.8 \times 10^3 }{44.1 10^3} = 0.2177 $$


Is image padding really necessary while using a filter? I do know that we must pad an image while using a filter. But what I am suggesting is that is it really necessary? Can't we just calculate at run time if the current cell/pixel is an edge cell/pixel and then use the padding value directly over there. It might add an overhead to the function but so does the entire padding function.

Only reason I think the padding function should be used is if the padded version of original image is going to be used in many different function which will overtake the overhead of dynamically calculating the edge cell/pixel values(which might never occur). Am I wrong?

Adding code to your filter function to detect and resolve edges implies code running needlessly throughout most of the image. Image area grows quadratically with image size.

Padding grows linearly with image size (as does the perimeter), so adding padding instead of edge-detection code seems a more efficient option.

Also, when adding padding, you are free to choose what to fill the padding with, according to your particular problem. You may choose to fill in with a constant value, or mirror the image, or get overlap data from adjacent images, etc.


How to make the $\ell_2$ norm of all columns and rows of an $n \times n$ matrix equal to $\sqrt{n}$? I have an $n \times n$ matrix and I would like its columns and rows to have $\ell_2$ norm equal to $\sqrt{n}$. Is this possible?

HINT If we have the diagonal matrix: $$ D = \left[\begin{array}{cccc}
d_1&0&0&0\\
0&d_2&0&0\\
0&0&\ddots&0\\
0&0&0&d_n
\end{array}\right]$$

Multiplying another matrix $$M_r = \left[\begin{array}{c}
r_1\\
r_2\\
\vdots\\
r_n
\end{array}\right]$$to the left with it multiplies each row like this:

$$DM_r = \left[\begin{array}{c}
d_1r_1\\
d_2r_2\\
\vdots\\
d_nr_n
\end{array}\right]$$

If the matrix has the columns $$M_c = \left[\begin{array}{cccc}
c_1&
c_2&
\dots&
c_n
\end{array}\right]$$

Then multiplying from the right with D will give:

$$M_cD = \left[\begin{array}{cccc}d_1c_1&d_2c_2&\dots&d_nc_n
\end{array}\right]$$


Librosa stft + istft - Understanding my output (which always seems too perfect) at varying window lengths I've just started to use Python with Librosa for a DSP project I'll be working on. First thing I've been trying to do is determine my preferred parameters for the FFT window size, and hop-distance.

The domain is music, and my plan is to try a variety of values for the window size and hop distance, and for each of them, do a forward STFT and then an inverse STFT and write the result back out to wav file. I'll then listen to results and choose based on which values I think capture the information in the input the best.

My simple code is as follows:

import librosa.core as lc
import numpy as np
import scipy

_n_fft=80
print(str(_n_fft))
_hop_length=_n_fft/4

data, sampleRate = lc.load("13_Hate_To_See_Your_Heart_Break.wav", sr=44100, duration=20, mono=True)

stftMat = lc.stft(data, n_fft=_n_fft, hop_length=_hop_length, center=True)
iStftMat = lc.istft(stftMat, hop_length=_hop_length)

scipy.io.wavfile.write("testOut.wav", 44100, iStftMat)

powerMat = np.abs(stftMat)
print("powerMat shape = " + str(powerMat.shape))

The behavior I'm experiencing, however, is not what I would have expected.

When I use an incredibly short window length (as in the code above) - I get the correct number of window frames for my FFT length and hop-distance:

powerMat shape = (41, 44101)

44101 window makes sense, and as you can see the frequency resolution is low, with only 41 frequency bands. I would expect the resulting testOut.wav to sound pretty terrible, as the frequency resolution is so low. I can visibly see the effects on a rendered spectrogram as the subtleties in frequency changes are smeared together. Listening back, however, the resulting track sounds great - pretty much like the original input.

Compare this with a much wider window size of 44100 (1 window = 1 second of audio, hop-distance of 1/4*Window size):

powerMat shape = (22051, 81)

Again this output makes sense - in the 20 seconds of audio, with a window length of a second and a hop distance of a quarter second, there would be about 80 window frames. This is pretty poor time resolution, but fairly high frequency resolution with 22051 frequency bins. Again I would expect the resulting testOut.wav to sound poor in the time domain.

Once again the resulting track sounds great - pretty much like the original input. These extreme values, and everything in between, pretty much yield the same output testOut.wav, even though on the real power spectrum I can visibly see the differences when changing the parameters.

Is there a fundamental misunderstanding I'm having with the STFT and it's inverse? Am I simply not understanding the library?

The information you think is missing, due to either poor visible frequency resolution or poor visible time resolution, is actually still contained in the sequence of FFTs, but not in a form that can be easily visualized, especially by just displaying FFT magnitudes or spectrograms.  That's because the information is "hidden" in the FFT phase results.

In the sequence of short FFT windows case, some of the frequency resolution information is represented in the difference in phase between adjacent or overlapped windows.  In the long FFT window case, much of the time locality information is represented in both the absolute phases and in the difference in phases within groups of nearby FFT result bins.  None of that is visible in FFT magnitude plots or spectrograms; and FFT phase and phase delta plots are hard to read.

But complete complex IFFTs do use the phase information contained in the real and imaginary components of the complex FFT results, and thus can reproduce both the time and the frequency information accurately and completely.

This phase information can be used in certain forms of FFT analysis to gain back some of the "missing" resolution.  For instance, the phase vocoder algorithm deduces or uses greater frequency resolution by looking at the phase deltas between adjacent or overlapped FFT windows.  And the time shift property of the FT shows how the phases between DFT result bins rotate relative to each other as an event (an AM or "beat" modulation, etc.) is rotated to different positions within a DFT window.


What does nonnegative zero-phase response mean? I am not exactly sure what nonnegative zero-phase response means. If a filter is zero-phase (i.e. symmetric and non-causal), then what does nonnegative imply? And what are the conditions to satisfy it?
Context: This is an error I get in Matlab when trying to design min-phase filters:
Error using firminphase>checknonnegative (line 79)
The filter has to have a nonnegative zero-phase response.
Error in firminphase (line 38)
checknonnegative(b);
Error in firgr>genfilter (line 323)
hRet = firminphase(hRet, s.zero, 'angles');

What they mean here is that the real-valued amplitude function of the linear phase FIR filter that you provide to the function must be non-negative, because it is interpreted as the desired squared magnitude of the minimum-phase filter. This amplitude function is referred to as zero-phase response, which is the frequency response of the filter when the coefficients are shifted such that they are centered around the origin, i.e., $b[n]=b[-n]$ is satisfied. Note that the filter length must be odd.
If $N$ is the (odd) filter length, the zero-phase response is given by
$$B_{zp}(e^{j\omega})=B(e^{j\omega})e^{j\omega(N-1)/2}$$
where $B(e^{j\omega})$ is the frequency response of the causal filter.


The Phase Information of the DFT As known, we can get a 'phase spectrum' from the DFT of a input signal. Assume, we do a DFT on a given equal interval sampled 50Hz AC signal, the signal is a complete whole cycle, but the start sample point may at any point. And I get the amplitude and phase from the DFT, but what's the reference point of the phase?

If you assume that for a given DFT length $N$ your signal is given by

$$x[n]=A\cos(2\pi n/N+\phi),\quad n=0,1,\ldots,N-1$$

i.e. its period equals $N$, then the first DFT coefficient (starting to count from $0$) is given by

$$X[1]=A\frac{N}{2}e^{j\phi}$$


High Pass Filtering and Noise I have the following code:

sigma_th = 0.08;
ph_noise = sigma_th*randn(1,2^22);
[b,a] = butter(1,4e6/1e9,'high');
ph_noise_out = filter(b, a, ph_noise);
figure(1);
[psd1,w] = pwelch(ph_noise,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd1/2));grid on;
hold on;
[psd2,w] = pwelch(ph_noise_out,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd2/2));grid on;
[psd3,w] = pwelch(ph_noise_out - ph_noise,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd3/2));grid on;

Here, I generate a noise signal ph_noise and high pass filter the noise using a butterworth filter to get the output ph_noise_out.
If I check the power spectral density of the ph_noise and ph_noise_out, I see that the output spectrum is filtered at low frequencies but is exactly the same at high frequencies. But when I plot the difference between the two noises, I see that the difference falls down as 20dB/dec. Why is the difference falling so slowly? How can I make it fall faster?
PS. The blue curve is the ph_noise, orange is the filtered noise ph_noise_out and yellow is the power of the difference signal.

The reason is the decorrelation due to the varying phase shift between the input and output. See the plot below showing the phase of the OP's high pass filter, and how closely it follows the cancellation when plotted on a log log scale (the higher frequency components of the two signals are better aligned in time and therefore have higher correlation leading to better cancellation).

The correlation between signals vs phase is proportional to the sine of the phase (signals in quadrature are uncorrelated consistent $\sin(\pi/2)=0$. For small angles the sine of the angle is the angle in radians hence we see the close match between the response of the phase and the cancellation of your signal as in your plots.

Due to the small angle, the OP's cancellation signal, assuming perfect amplitude matching, will have a magnitude that is directly proportional to the phase itself in radians between the two signals. Consider the cancellation of two signals with identical magnitude but a difference in phase:

$$y_1= Ae^{j\phi_1}$$
$$y_2= Ae^{j\phi_2}$$

$$y_1-y_2 = Ae^{j\phi_1}- Ae^{j\phi_2} = Ae^{j\phi_1}(1-Ae^{j(\phi_2-\phi_1)})$$
$$= Ae^{j\phi_1}e^{j(\phi_2-\phi_1)/2}(e^{-j(\phi_2-\phi_1)/2}-e^{j(\phi_2-\phi_1)/2})$$
$$ = -Ae^{j\phi_1}e^{j(\phi_2-\phi_1)/2}2j\sin\big((\phi_2-\phi_1)/2\big)$$
$$ = -2Aje^{j(\phi_1+\phi_2)/2}\sin\big((\phi_2-\phi_1)/2\big)$$
$$ = 2Ae^{j(\phi_1+\phi_2-\pi)/2}\sin\big((\phi_2-\phi_1)/2\big)$$

For a small difference angle $\phi_2-\phi_1$, $\sin\big((\phi_2-\phi_1)/2\big) \approx (\phi_2-\phi_1)/2$, so the above cancellation becomes:

$$\approx A(\phi_2-\phi_1)e^{j(\phi_1+\phi_2-\pi)/2}$$

And thus for small angles the magnitude of the cancellation is directly proportional to the difference in phase. And the cancellation ratio as the magnitude of the cancelled signal compared to the magnitude of the input signal would be the difference in phase directly:

$$|y_1-y_2|/|y_1| =| \phi_2-\phi_1|$$

To achieve what the OP is trying to do, the phase response of the filter must be compensated for before doing the subtraction. This would be much easier if a linear phase filter was used instead of the Butterworth filter as the phase compensation of the input signal could then be accomplished with a simple single tap delay. Of course, to completely cancel the signal the amplitudes would need to be exactly matched in the region of cancellation as well.

In practical application even with close matching a finite delay will exist between the two identical signals prior to subtraction, resulting in a high pass response with a cutoff of $1/T$ given a finite delay $T$. With a fixed delay $T$ the lower frequency components much less than $1/T$ in frequency would have a higher cancellation, and the decorrelation will go proportional to $f$ as the frequency approaches $1/T$. Signals much higher than $1/T$ would be decorrelated so would result in a $3$ dB increase as expected for the the sum of two uncorrelated signals.


Compare two words from different speakers I have a set of audio files which contain a word each one, and I'd want to be able to compare a new word (coming from a smartphone mic) to check if it resembles to a specific file and how much it resembles to that audio file. I don't need to know which word it is, or doing a STT system, actually, just checking how similar two words are. The words may come from different speakers, male, female, old people, kids, even non-speakers for a specific language.

Now, the question is, which should be the processes to do something quite simple? I've been reading about MFCC, DTW and I can ensure I understand little to nothing about how I should proceed. I'd need to create an algorithm to match those words, but I have not a single clue on how to proceed after the "Conception steps".

I have, however, tried an Android library called Musicg, which should compare sounds based on:

A frequency range for the whole sound
An intensity range for the whole sound
A standard deviation range (yes, for the whole sound)
A High pass and low pass limits
Range of times the sound crosses the "zero" line

But I am pretty sure this might not be enough to make sure two sounds are similar or different. The examples they handle are only onomatopoeii, (for those who code, here's an example to detect Claps).

Now, the question. Is there a magic algorithm that does all the heavy work of comparing two audios from different speakers that is able to give a similarity score? If so, can you put me back to the tracks with a simple (well, as simple as it can be, of course) example, please?

Bonus track: Is any of the algorithms able to do the processing with precomputed values from the set of audio files?

There is no (known) magic solution, as the sound of the same word from very different speakers are often not actually similar to one another in term of any simple characteristics such as you list (frequency range, intensity range, pitch, etc.)

To often, humans can only guess at what word is said from an unknown speaker based on context.  And you don't have that with single word recordings.


How can I find this all pass filter's phase? How can I find this all pass filter's phase?

$$H(e^{jw}) =\frac{−a* + e^{-jw}}{1 - a(e^{-jw})}$$
where $a = ρe^{jϴ}$ is the pole.

Then I want to find group delay, I know that the group delay is the derivative of the phase with respect to w but i don't know how to find the phase

phase of a filter is the angle (or "argument") of the complex number that is $H(e^{j\omega})$ at each value of $\omega$ .  you may need to "unwrap" the phase by adding or subtracting integer multiples of $2 \pi$ .

BTW, to post math to these pages, surround the $\LaTeX$ expression with single dollar signs for inline math and double dollar signs for math expressions on their own line.  if you need a quick guide for LaTeX, my suggestion is the one that lives at Wikipedia: https://en.wikipedia.org/wiki/Wikipedia:Math .  i'm still learning how to do other markup here.


How to get rid of noise and oscillations using Kalman filtering I have not studied signal processing at all, so please forgive any ignorance in the following discussion.

I have some noisy position measurements that I've been trying to smooth. I've been attempting Kalman filtering but getting either hf noise or lf oscillation depending on the constants I put in. What should I do?

The full story:

I've used a program which is basically like this (I read briefly about how the Kalman filter works and made my own version):

for z_measurement in samples:
    expected.r = old_estimate.r + old_estimate.r * delta_t
    expected.v = old_estimate.v + old_estimate.a * delta_t
    error = z_measurement - predicted.r
    new_estimate.r = old_estimate.r + gain.r * error
    new_estimate.v = old_estimate.v + gain.v * error

    old_estimate = new_estimate

This worked fine for constant acceleration, but when I have a slowly changing (uncontrolled) acceleration (drag which appears slowly decrease) as another state variable, I had a trade off. As I lowered the position gain, I got rid of high frequency noise, but I got even worse low frequency noise which overwhelmed the signal in the acceleration. The position gain seemed to act like shock absorbers, but let through too much high frequency noise.

I've been trying to learn real Kalman filtering to see if that could help. I do not have a control system, so I'm ignoring B and u. It seems to turn into a computer program like this:

for z_measurement in samples:
    x_prediction = A_update * x_old_estimate
    P_int_covariance = A_update * P_old_covariance * A_update + Q_process_error
    K_gain = P_int_covariance * H_measure * (H_measure * P_int_covariance * H_measure + R_measure_error)**-1
    x_new_estimate = x_prediction + K_gain * (z_measurement - H_measure * x_prediction)
    P_new_covariance = (Identity - K_gain)* P_int_covariance

    x_old_estimate = x_new_estimate
    P_old_covariance = P_new_covariance

This works fine for calculating averages, but looks non-trivial to define P, Q, and R right for a real system.

From what I now understand of Kalman filtering:

A is the prediction model of the system. It corresponds exactly to my previous program:

expected.r = old_estimate.r + old_estimate.r * delta_t
expected.v = old_estimate.v + old_estimate.a * delta_t

H is a meaningless function used to access the measured parameter.

expected.r

Q and R are error matrices that define the expected noise/error at each stage, and P defines how error of one state variable increases the error of other variables.

However Q and R are constants, and P and K_gain quickly converges to constants when Q exists. In fact K_gain is the only useful matrix (the others are only used to calculate the eventually constant K_gain).

It seems to me that rather than using dark arts to guess Q and R, you can do it directly by playing with K_gain until it works.

Which a few days later leads me back to my original program.

...

Will persisting and correctly defining Q and R, (or K directly) actually give me gain constants that get rid of my noise and oscillations, or would I be better off trying non-linear gain functions. Alternatively should I put up with the high-frequency noise and use a separate low-pass filter afterwards, but then would it be better to do simple low pass filtering and totally forget the Kalman filtering?

Are there any other filters I should try? I'm afraid that linear filtering will loose too much signal. I tried some quadratic filtering and ended up with oscillations with periods related to the length of the filtering window.

First, I think your motion equations are wrong. Assuming a constant acceleration model, the transition matrix, $\mathbf{A}$, should look like that:
$$
\mathbf{A} = \left[\begin{array}& 1 & \Delta_t & 0.5\Delta^2_t \\ 0 & 1 & \Delta_t\\ 0 & 0 & 1\end{array}\right]
$$

Next, for your question, if the measurement noise (defined by $\mathbf{R}$) and the model noise (defined by $\mathbf{Q}$) are Gaussian, then the Kalman filter is the linear MMSE. That means that it is the optimal linear estimator in the MMSE sense.
Playing directly with the filter's gain yields a simple exponential averaging of the data.

If you have access, I would recommend Tracking and Data Association by Bar-Shalom and Fortmann. This is an excellent book that covers Kalman filtering.


Sampling signal after root raised cosine I am generating series of 1 and 0. then pulse shaped them with root raised cosine and after than deciding to do match filtering and recovering the bits however it turns out that I am not getting the correct bits.
What I do is as follow:

generate 0 and 1.
pulse shape them with RRC.
match filter using the same transmited filter
Remove samples up to peak using [val,pos]=max(mysignal)
sample using mysingal(pos:overample:end)

But what I get is not correct. I simulated this using 10 bits but I get only 2. I know there is something wrong but I do not know why. My method works when I only have one bit.

I am wondering if I have more than one bit should I do this procedure for every  bit.How can I approach this problem.

This is also my code assume RX2 is the pulse shaped version of my signal and os is the oversampling factor and 0.25 is the roll of factor. I am wondering how to detect all the bits like all 10 or 100?

RX1=filter(srrc(os,0.25),1,RX2);
[val,pos]=max(RX1);
RX=RX1(pos:os:end);
Bits=RX;
One=Bits>=Th;
BER=sum(One)/length(Bits);

The first one is the signal after passing through matched filter using conv command and no built in command like filter. the original bits are 0 0 1 1 1
Can anyone tell me how can I sample signal to get these bits back.
what command I should use. simply downsampling will not work.

The second one is for the same bits but using filter command for 0 0 1 1 1 can anyone please tell me how to obtain correct bits back from any of these figures.

I can visibly see from your figure that you have done the two cascaded RRC's properly. I circled the samples I would expect to see for your 00111 data sequence after RC filtering. Clearly you are using 8 samples per symbol, and a longer data sequence would continue in similar fashion after the initial delay that you have as shown.

To demodulate this, you typically would have a timing recovery loop to establish the ideal sampling locations. There are several approaches to establish timing, but one recently discussed on this forum is the use of the Gardner TED (Timing Error Detector) which would create an error signal that you could use to control a timing sample location (increment up or down a position pointer based on the accumulated error, forming a timing recovery loop. This would be done using a (discardable) training sequence with lots of transitions, typically as part of the header for the data packet that would be significantly longer than your short data set.

Another synchronizer I like is the M&M synchronizer (Mueller & Meuller) as well as excellent approaches I learned from fred harris using polyphase filters for timing synchronization (thank you fred!).

Below are links to recent discussions on the Gardner TED, but in particular I want to copy here the figure for the eye diagram of the transmitted symbols after they have gone through a root raised cosine (RRC) filter ("Before RRC Filter") the received signal after it has been passed through a second RRC filter as you are doing. Conveniently this was done with same alpha= 0.25 roll-off factor you have used, and what it shows is the continuous time representation of every possible trajectory between symbols based on previous symbol patterns. In your case, you have eight samples per symbol, landing on one of these trajectories.

Links to recent discussions on the Gardner TED:

Location of Matched Filter

Gardner Timing Recovery for Repeated Sybmols


How do I properly do gamut mapping for an image? I asked this question earlier. To summarize, dithering an image with a specific palette was causing the image to look smeared and very discolored.

Thanks to Olli Niemitalo, I found out my issue was because the colors in the original image (particularly the near-blacks and near-whites) were outside the gamut of the palette I had. As a result, Floyd-Steinberg dithering was piling up tremendous error levels which were being carried over to the next pixels. I confirmed this to be true by replacing every out-of-gamut color with the nearest color from the palette (which are guaranteed to be in the gamut, by definition). Trying to dither after doing this yielded better results.

What I'm having trouble with is finding the nearest color inside the gamut of my palette, given another color. Here's my current method:

Convert palette colors from sRGB to XYZ.
Interpret them as a collection of 2D points on a grid (discarding Z)
Draw a convex polygon around the points to determine the gamut.
For each pixel in the image

Convert the pixel from sRGB to XYZ
If the pixel's XY values are not in the polygon, find the nearest point that is inside it.
Convert this result back to sRGB (using the new X and Y values, and the original Z value)
Replace that pixel in the image with the new color

After all that, dither with Floyd-Steinberg.

This is the palette I'm working with (in (R,G,B) format):

(162, 143, 66)
(128, 148, 168)
(120, 99, 100)
(31, 39, 97)
(126, 116, 103)
(203, 35, 9)
(57, 81, 43)
(104, 101, 98)

Original image:

Gamut mapping (attempt):

Dithering:

That certainly looks better than the original problem I was having, but it's still troublesome. I'm sure I'm getting the gamut mapping procedure wrong somehow, and I suspect it has to do with discarding Z, but I have no idea how to proceed. Googling this yielded only Photoshop instructions, which isn't what I'm looking for.

Wikipedia page Gamut puts it nicely:

  Generally, the color gamut is specified in the hue–saturation plane,
  as a system can usually produce colors over a wide intensity range
  within its color gamut; [...]

But you don't have the general case. Without black for example, you can't produce darker shades of all colors even if they are within the gamut projected onto the XY plane. So you will need to determine the 3-d gamut polyhedron. The gamut is a polyhedron in an intensity-linear (gamma=1.0, linear gamma) RGB or XYZ color space, because intensities of color components mix additively in additive color mixing.

You have additive color mixing in the kind of dithering that you do. If someone looks at the picture from afar, they only see the average intensities of the color components of the pixels. This is not the same as for example average squared intensities, so gamma matters. In dithering, the error calculation and the residual error diffusion arithmetic should be done in an intensity-linear color space.

3-d gamut mapping is not going to be trivial. First you need to find the convex hull of the set of color coordinate points in your palette, in an intensity-linear color space. Then you need to check if a wanted color is inside or outside the convex hull. If it's outside, you need to map it to the surface of the convex hull, which is (the union of) a set of triangles in 3-d space. Here the easiest mapping might be to find the nearest point on the convex hull.

Better results might be had by using a perceptual color difference metric like CIEDE2000. However, this makes finding the nearest point a non-trivial optimization problem, possibly non-convex. Possibly the perceptual color difference metric would also give better results in dithering, just at the step where you choose from the palette the color to use.


Filter for DC component in FFT2 of image I am trying to get the main frequency of the Fourier Transform of the image below, which is the intereference pattern detected by a CMOS sensor.

When I apply fft2, I got the DC peak at 0 frequency, that I know it can be suppressed if I do FT = fft2(image-mean2(image)). However, I only obtain a 0 in the (1,1) position, but the "neighbor" frequencies in the first row and column have still a higher magnitude that the peak I want to detect, which is in the (6,11) position.

Is there a formal way to filter the image, meaning without imposing FT(1,:) = 0, FT(:,1) = 0, in order to be able to extract directly that peak?

When processing noisy data, it is usually beneficial to filter high frequencies, this can be done with a small kernel as shown in answer by Dan Boschen.
Also if your data has an envelope it will give you some high intensity low frequencies, the low frequencies require a long kernel to be filtered, and it is better to do in the frequency domain.
The topic of filter design has been extensively explored in the past but for 1D filtering (time domain). I use the 2D information to construct a direction independent filter.
Also, since the idea of this filtering is just to emphasize certain frequencies and we don't need to implement the transfer function in any recurrence equation, I just multiply the magnitude of the filter leaving the phase unchanged (I know your problem).
img = img - mean(img(:));
IMG = fft2(img);
[M, N] = size(IMG);
figure, imagesc(abs(IMG(1:64, 1:64)));

The frequencies you are interesteed in are problably those in the second spot close to 10 in the vertical axis. But we have a some frequencies that are detected higher represented at the top left corner of the image. Those need to be attenuated so that you can capture the frequencies of interest.

I will create a transfer function for a wide bandpass filter by cascading a
a high pass filter with cutoff in 10/v, and a lowpass filter at 200/v
The low frequencies are due to the envelope
f1 = 10/(N+M);;
f2 = 200/(N+M);
H = @(f) abs(f1 ./ (f1 + 1j*f) .* (1j.*f) ./ (1j*f + f2));
% see the filter response
x = linspace(0, 0.5, 1000);
figure, loglog(x, H(x));

Assuming each pixel corresponds to a square in the space (a rectanble with width = height) we can compute a frequency measured in cycles per pixel (where pixel corresponds the physical distance that separates two pixels)
% Spatial frequencies measured in cycles per pixel in each direction
phi_y = angle(exp(2j*pi*(0:M-1)/M))/(2*pi);
phi_x = angle(exp(2j*pi*(0:N-1)/N))/(2*pi);

% Construct a matrix of freq
F = sqrt(phi_y'.^2 + phi_x.^2);

If you have the actual length of each pixel you can simply divide phi_x and phi_y by the corresponding lengths.
As I said we have two options, one is to apply the analytic filter, but that changes the phase, or simply apply a real filter, by multiplying the magnitude.
Let's double check if this filtering does not modify too much the images by plotting them.
figure,
subplot(221);
imagesc(real(ifft2(IMG .* H(F))));
colormap('gray'), colorbar();
title('Analytic filter');
xticks([]);
yticks([]);
ylabel('real');

subplot(223)
imagesc(imag(ifft2(IMG .* H(F))));
colormap('gray'), colorbar();
xticks([]);
yticks([]);
ylabel('imaginary');

subplot(222);
imagesc(real(ifft2(IMG .* abs(H(F))))), colormap('gray'), colorbar();
xticks([]);
yticks([]);
title('Magnitude filter');

subplot(224);
imagesc(imag(ifft2(IMG .* abs(H(F))))), colormap('gray'), colorbar();
xticks([]);
yticks([]);

In my opinion the real filter is good, it does not leak the signal to the imaginary part and keeps the appearance of the image the same, the information you want is still there. If you want want you can increase f1 so that it filters the low frequencies even more, and you could for instance multipy by the square of the filter magnitude.
Now the last thing is to check if the peak we expect is the higher than the ones at the top left corner.
figure, imagesc(abs(IMG .* H(F))(1:64, 1:64)), colorbar();


Convert RGB Image to Grayscale and Display It (Python + Matplotlib) I just started learning image processing and I was trying to read a RGB image then convert it to grayscale. I was hoping for something like this:

However, what I get was:

I tried using both scipy and PIL but they yield the same results. Am I lacking of understanding about grayscale image here?

Using scipy:

from scipy import misc
car = misc.imread('image.jpg', mode="L")
plt.imshow(car)

Using PIL:

from PIL import Image
pil_im = Image.open('image.jpg')
pil_imgray = pil_im.convert('LA')

img = np.array(list(pil_imgray.getdata(band=0)), float)
img.shape = (pil_imgray.size[1], pil_imgray.size[0])
plt.imshow(img)

You need to set the cmap properties of the imshow() function:

https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.imshow.html

Use:

plt.imshow(img, cmap = "gray")


Signal and Space/Time I have been reading a couple of definitions on signal such as:

  Any quantity measurable through time over space or any higher
  dimension can be taken as a signal.

And this,

  In the physical world, any quantity exhibiting variation in time or
  variation in space (such as an image) is potentially a signal that
  might provide information on the status of a physical system

Taking human voice as example,

A sound signal is converted to an electrical signal by a microphone, generating a voltage signal as an analog of the sound signal, making the sound signal available for further signal processing. Sound signals can be sampled at a discrete set of time points.

I can understand where times comes in but where's the "space" come here according to the definition of signal?

My understanding of your question is that you have a small misunderstanding here. In the definition it is not said time and space. A signal can vary in time or space. Some signals vary with time, as your example human voice that varies over time in air pressure (or equivalently voltage), some vary with space, like image and some vary with both time and space such as video. I feel,  your misunderstanding in human voice example is that you believe there should be an space dimension too in the game.


Restricted Isometry Property (RIP) in Compressive Sensing What is the meaning of Restricted Isometry Property (RIP) condition in Compressive Sensing for Sparse Signal Analysis?
How can we define Restricted Isometry Constant (RIC) for the RIP condition?

Thanks in Advance!

The restricted isometry property states that:
$$\begin{equation}
  (1-\delta_S)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_S)||x||_2^2
\end{equation}$$
for any $S$-sparse vector $x$. The restricted isometry constant is $\delta_S$, $0 < \delta_S < 1$.

This means that the matrix $A$ is guaranteed to only change the length of any vector $x$ "very little" as long as the vector $x$ is at least $S$-sparse (has at most $S$ non-zero coefficients).

Suppose we have arbitrary $\frac{S}{2}$-sparse vectors $x$. In order to be able to reconstruct such vectors in general, from measurements taken as $y = A x$, we need to be sure that it is possible to distinguish between measurements $y_1 = A x_1$ and $y_2 = A x_2$ of any two such vectors. If $y_1 = y_2$ for any two such vectors $x_1$ and $x_2$, we would not be able to distinguish them and reconstruct them unambigously. Therefore we need to make sure that measurements of any two $\frac{S}{2}$-sparse vectors are "sufficiently different".

If we calculate the difference between any two $\frac{S}{2}$-sparse vectors, their difference can be at most $S$-sparse. So for reconstructing any $\frac{S}{2}$-sparse vector correctly from measurements taken with $A$, the restricted isometry property quantifies how well $A$ lets us do that (the smaller $\delta_S$, the better).

For an early introduction to compressed sensing and the restricted isometry property (and other concepts), see Candès & Wakin, 2008.


Generating a noisy sine wave in Python, efficiently (Questions after code block.)

For now, I am generating the sinusoid and the noise separately and then normalising the amplitude so I can write the signal to a .wav file (16-bit, 16kHz).

def genSine(f0,fs,dur):
    sinusoid = []
    for t in range(dur):
        sinusoid.append(math.sin(2*math.pi*t(*freq/fs)))
    sinusoid = normalise(sinusoid)
    return sinusoid

def genNoise(dur):
    noise = np.random.normal(0,1,dur))
    noise = normalise(noise)
    return noise

def normalise(x,MAX_INT16):
    maxamp = max(x)
    amp = floor(MAX_INT16/maxamp)
    norm = np.zeros(len(x))
    for i in range(len(x)):
        norm[i] = amp*x[i]
    return norm

def writeWav(w):
    # audiophilic magic

if __name__ == '__main__':
    f0 = 440
    fs = 16000
    dur = 1*fs                      #seconds
    MAX_INT16 = 32767
    sinusoid = genSine(f0,fs,dur)
    noise = genNoise(dur)
    sum = x + y
    sum = normalise(sum,MAX_INT16)
    writeWav(sum)

Is there a more efficient way to sum the two signals (sine + noise), perhaps bypassing/incorporating the normalisation step (it is currently called three times, in genSine, genNoise and main)?
How can I ensure set the amplitude ratio between the sine and noise signals?

I'm new to Python and stackexchange so any help is appreciated!

The are some ways to improve your solution.

First, use numpy and vectorize your code:

def genSine(f0,fs,dur):
    sinusoid = []
    for t in range(dur):
        sinusoid.append(math.sin(2*math.pi*t(*freq/fs)))
    sinusoid = normalise(sinusoid)
    return sinusoid

becomes

import numpy as np
def genSine(f0, fs, dur):
    t = np.arange(dur)
    sinusoid = np.sin(2*np.pi*t*(f0/fs))
    sinusoid = normalise(sinusoid)
    return sinusoid

You will see the biggest improvment in efficiency when vectorizing all your code like this.

Second, concerning the normalization, you should consider the mathematics of your problem instead of doing it 3 times. You have 3 pieces of signal, the sinusoid, the noise, the mix. The first is already normalized (between -1 and 1), the second also (hence the name normal distribution for describing the noise). So you should only normalize the third.

Hope this helps!


Frequency vector and fft I have got a question concerning the definition of the frequency vector for an fft operation.

Generally, I work with a frequency vector, f, with power of 2 elements (2048, 4096, 8192, ...).

Given a certain simulation analysis time, time (e.g. 600s), I should define f as follows:

% Frequency defition
t = 0:dt:(time-dt);
df = 1/(time);
fn = Nfft/time;

$$ f = -f_{n}/2:df:f_{n}/2-1; $$

where $ f_{n} $ represent the Nyquist cut-off frequency.

Actually, for:

computational reasons
symmetry of the power spectra along f axis
not throwing away real or imag part of the fft

I aim to define only half of the frequency range, as for example

$$ f = 0:df:f_{n}/2-1; $$

After calling the fft of my input signal, I would get the desired time series as

ouput = [real(fft) imag(fft)];

But, this way, I count the 0 frequency term twice and the -fn/2 is completely discarded.

How would it be possible to emcompasses the whole standard frequency range starting from only half of it?

The proper way to define your frequency vector after a DFT is as follows. Let $N$ be your DFT length, and $f_s$ be your sampling rate in Hz. Furthermore, define an $N$-length frequency vector $\bf{f}$, where each element $f_i = i$, for $i = 0, 1, 2, ... N-1$.

Now your frequency vector in hertz is simply going to be $\bf{f}$$\frac{f_s}{N}$

Now, assuming you are DFT'ing a real sequence, simply pick all elements with frequency values less than $\frac{f_s}{2}$, and you are in business!

You can also see my answer here for actual code.


What Is the Best Color Space to Work with for Hand Detection: RGB, HSV, or YCrCb? I'm working on hand gesture recognition using C#. I'm using OPENCV for skin detection but its performance is not optimal. So I decided to do some filtering based on color space, but I do not know what color space is the most suitable for this purpose?

Intuitively, HSV is the place to easily define Skin Color Hues.

Yet there is a broad work on that and even articles about the optimal Color Space for Skin Detection.

Yet, you should have a look at OPTIMUM COLOR SPACES FOR SKIN DETECTION (Alternative at IEEE - Optimum color spaces for skin detection).

According to them there is no difference in the performance as long as you create the Optimal Detector.
This actually makes sense, since all Color Spaces in the article are One to One functions it means the real thing here is to develop the correct classifier.


Why DFT exact for periodic signals? I have in my course notes:

$$
\sum_{k=0}^{N-1}y(k)e^{-j\omega_nk}\approx\sum_{k=-\infty}^{\infty}y(k)e^{-j\omega_nk}
$$

Where the $\approx$ for some reason becomes $=$ when $y(k)$ is periodic. Could you please explain why this is the case?

Thank you so much!

It's not correct to say that both expressions become equal when the signal is periodic. What is correct is that - obviously - everything one can know about an $N$-periodic signal is represented by the expression on the left of your equation, and consequently, both expressions are equivalent in the sense that one can be obtained from the other.

Note that that the infinite sum on the right-hand side does not converge in the conventional sense if the signal is periodic. The Discrete-Time Fourier Transform (DTFT) - that is what the right-hand side represents - of an $N$-periodic signal can be expressed by using Dirac delta impulses:

$$X(\omega)=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X[k]\delta\left(\omega-\frac{2\pi k}{N}\right)\tag{1}$$

The coefficients $X[k]$ can be obtained from the expression on the left-hand side of your equation be setting $\omega_k=2\pi k/N$. In that sense, both expressions convey the same information.


Show That a 2D Linear Transform $ T \left( \cdot \right) $ Is Homogeneous By my understanding, a transform T is homogeneous if T[0] = 0.

Then to prove that a linear transformation is homogeneous we say that:

T[ax(n1, n2) + bx(n1, n2)] = aT[x(n1, n2)] + bT[x(n1, n2)]

What I want to know is the difference between the following two proofs and which one is "more correct" or is there perhaps a even better answer?

let x(n1, n2) = 0

then we know that,

ax(n1,n2)- ax(n1,n2) = 0,
T[ax(n1, n2) + ax(n1, n2)] = aT[x(n1, n2)] + aT[x(n1, n2)] = 0

VS

Simply showing that both sides are equal to 0 if x(n1, n2) = 0 through algebra.

On a side note, can someone also comment on why we would even care if a transformation is homogeneous vs heterogeneous? In reference to a transformation applied onto a 2D signal (such as an image transformation).

Maybe the better question to ask would be to prove if a transformation is heterogeneous, is it linear? (if so can someone help me define what it means for a transformation to be heterogeneous? Like how a homogeneous transform means T[0] = 0.

The idea here is that any Linear Operator is Homogeneous though not every Homogeneous Operator is Linear.

The classic proof indeed is to build the zero term as a sum (Difference) of 2 other elements in the domain. Then use Linearity to show it must be zero:

$$ T \left( 0, 0 \right) = T \left( {m}_{1} - {m}_{1}, {n}_{1} - {n}_{1} \right) = T \left( {m}_{1}, {n}_{1} \right) - T \left( {m}_{1}, {n}_{1} \right) = 0 $$

  Maybe the better question to ask would be to prove if a transformation is heterogeneous, is it linear?

As I wrote above, if an operator is heterogeneous is must not be linear.

Think of the simplest Linear Operator in 1D: $ T \left( x \right) = a x $. Clearly it is a linear function and indeed homogeneous.
Let's look on Affine Function: $ T \left( x \right) = a x + b $. Usually this is also called Linear Function while it is not unless $ b = 0 $. Why? It is not heterogeneous.

It is easy to see that $ T \left( {x}_{1} - {x}_{1} \right) \neq T \left( {x}_{1} \right) - T \left( {x}_{1} \right) $.


Delay signal according to normal distribution I am trying to delay a signal by a certain number of samples. This certain number, I would like to generate using a random normal distribution in Matlab. I am using the following for generating the number of samples r = randn(1,1)/sampleTime. This however generates negative numbers too. In order to avoid the negative numbers, I am taking the absolute value of r but I know that this is not the right way, since in this case I would be ignoring the left part of the distribution. Is there a way I can use just positive samples without shifting the mean too much?

Normal distribution regardless of whatever mean and variance you provide, can still produce large positive and negative numbers(albeit with a small probability).

Do you want to specifically use normal distribution, or all you need is a random delay generator? If latter is the case, you can use rand() from MATLAB, also you could use randi if you need integer delays.

Also you could use something like this to generate uniformly distributed random numbers in a specific range.


Can playing a wav file theoretically damage my speakers? I've been fooling around with some audio transformations and it's useful for me to listen to them when debugging.  Is it possible to create/play a wav file that will damage my speakers?

I'm working on a macbook pro and using the python module scikits.audiolab to play the sounds.

Generally yes, but only if the amplifier is oversized for the speaker or there are other design faults. Active speakers where the amplifier needs to be generously sized to deliver the expected momentary performance are usually protected by peak limiters and filter circuits.

For laptop speakers the probability of damage is extremely low.


extract the tree structure from a noisy data set I am currently doing a project that attempts to reconstruct the 3D structures from slices, i.e. the tomographic approach.  After many steps of image processing, the last task is to connect the “points” (the centroid of the sections) to create the skeletons. The skeletons have been confirmed to be a forest of binary trees.

The problem here is that the data are very noisy. Would you please suggest some techniques to perform data fitting to extract the trees?

I guess a very good resource about that is here (Implementation might be tough, but even reading it will give you a lot of insight about the extraction of 3D centerlines)

Also, please take a look at 1 and 2.


Extracting sound pressure from wav File I have recorded an ultrasonic wave from a microphone using MATLAB command. After plotting that recorded wav file, the graph shows the amplitude between $1$ and $-1$.

Can anyone tell that what is does mean?
How can we convert that data into pressure?

I read from the detail of audiowrite() here that the data should be normalized between $1$ and $-1$. But my question is in what form it reads and stores in the wav file. Either it is voltage or something else.

Why it is normalized?

I'd say it's extremely difficult to obtain the sound pressure from a WAV file.

Consider what happens to the sound signal between the microphone and the WAV file. Pressure is applied to the mike. It is converted to an electrical signal, with power that depends on the properties of the mike transducer. Then it is attenuated on its way to the sound card, where it is amplified, filtered, sampled, quantized, encoded, and further massaged by the sound card driver and the OS sound system.

(Furthermore, if the signal is ultrasonic, the frequency response of your sound-card front-end may very well come into play too).

Then, in software, the samples are further processed to fit a particular WAV format (the details of which are amply documented in Wikipedia and other sources).

One possible way to proceed is to obtain experimental readings with your particular set-up: generate a set of known, controlled sound pressures on your microphone, and correlate them with the numbers you get in the WAV file.


low pass filter after fft in matlab I'm new to signal processing and I'm learning on my own. Therefore I might ask questions not worthwhile your time.

At the moment I'm trying to figure out the concept of a lowpass filter.
Using matlab I created a complex signal x:

x = rand(1,1000)+rand(1,1000) .* 1i

And i created a filter Num with the fdatool, the specifications of Num are :

Fs = 1000;  % Sampling Frequency

Fpass = 460;             % Passband Frequency
Fstop = 500;             % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor

From what I understood I had to do now. I tried to do the Convolution of Num and x in the frequency domain by multiplication of the vectors.  Step by Step I did this:

First I padded Num to the same length that x has.
Then I multiplied x and Num, to get y= Num * x;

Now y should have nearly no frequencies above 500, right?
If I check the signal there doesn't seem have happened a lowpass filtering.

It appears you're getting a bit lost between the time and frequency domains. If "Num" is the numerator of an FIR filter (as it appears), these are the coefficients (or impulse response) of your filter. The impulse response as a signal exists in the time domain.

If you want to perform a filtering operation using this impulse response (or filter kernel as some call it), you should perform a convolution (conv) between it and the signal you wish to filter in the time domain.

--

The impulse response you've created has an equivalent frequency domain response (which you will have seen in fdatool).  The time domain convolution described above is equivalent to multiplying this frequency-domain response by the frequency-domain representation of your input signal.

Since your signal x is complex I assume you mean it to a frequency-domain representation of some 'noise' signal. This representation in itself is strange, and it doesn't represent a 'white' noise signal, but a purely random arrangement of frequency and phase that I find it difficult to conceive.

If you want to observe low pass filtering I would suggest the following:

Create a pure real (no 'i' component) noise signal in the time domain using rand (remember that rand produces values in the range 0->1, you will need to scale and shift it to create a signal in the range -1 -> 1).  This signal will in theory be purely white, and hence have a flat magnitude response in the frequency domain (i.e. it would be a straight line).
Convolve this signal with your filter impulse response.

You could consider this operation to be equivalent to multiplying the response of your filter (that you saw in FDATool..) by the frequency domain of the input signal (which is effectively constant). The result should be a scaled version of your familiar low-pass filter response (although some other effects such as the trade-offs in your filter design come into play).

The output of this procedure is still in the time domain, so if you perform an FFT or some other transform to convert to the frequency domain, you can then plot it, and observe the low-pass filtering effect.

Good luck for your learning :)