Converting a discrete real time domain to a discrete complex time domain signal To test my understanding of amplitude modulation in the frequency domain I'm trying to replicate the center frequency with 2 side-bands in an Octave plot with the following code

t  = [ 0 : 1 : 4095];
Fs = 2500000;
Fc = 1000000;
Fsig = 100000;
carrier = 20*sin(2*pi*Fc/Fs*t)';
signal = 10*sin(2*pi*Fsig/Fs*t)';
modulated = carrier .* signal;

plot(abs(fft(modulated)))

And I'm getting the following plot

I suppose I'm getting the 2 peaks mirrored instead of a single central peak with 2 symmetrical side peaks because it's an FFT of a real signal and I need a complex signal instead.

How do I convert the real time-domain signal to a complex time-domain signal assuming the center frequency is Fc?

Edit:

To state the question more clearly.
How do I get an FFT which will have one peak at the center frequency and 2 side peaks (3 peaks total) from my AM signal above?

I understand that I'm getting what I'm getting because the AM signal is real.

So I suppose I need a complex time-domain signal to feed to fft instead.

Is it possible to convert my real AM signal to a complex signal somehow for the FFT of the complex signal to look similar to the below (photoshopped)?

I suppose that's what quadrature sampling does in a receiver, but how do I do this in the digital domain?

My updated answer to OP's update

To do this specifically for the case of the FFT of a product of two sinusoids, you must choose a sampling rate or modulation frequency such that the upper sideband of the "positive" frequency component aliases to be on top of the lower sideband of the "negative" frequency component. For instance using a carrier frequency of 1.2MHz and a signal frequency of 50KHz would give you one peak in the center with two sidepeaks. (Although to note with 4095 points this wouldn't match the plot position since the center bin would be 2047 which is closer to where the lower sideband is in the plot).

Otherwise for an arbitrary signal the fft of the following would provide this result:

$$s[n] = Ae^{j\omega_c n}B\cos(\omega_m n)$$

Where

$\omega_c$ is the normalized radian frequency of the carrier ($0$ to $2\pi$)

$\omega_n$ is the normalized radian frequency of the modulation ($0$ to $2\pi$)

$A$, $B$ are arbitrary gain constants for each term (which simply scales the result by $AB$)

The opening part of the question was you are doing this to enhance your understanding of the frequency domain of AM modulation. So I am assuming you really want help understanding what AM modulation is supposed to look like:

When you multiply two sinusoids to get AM modulation, the result is DSB-SC (Double Side-band Suppressed Carrier): The carrier frequency is suppressed and there will only be two sidebands. This is what a cosine wave looks like in the FFT: two sidebands centered on the DC bin (use fftshift command if you want to center DC in the middle of the plot, then showing a positive and negative frequency axis which may help intuition), but since there is no DC in the cosine, that bin is supressed. So a cosine on its own is a DSB-SC signal with the carrier being DC, which is suppressed (we call that the baseband signal).  If you multiply this baseband signal with a exponential carrier $e^{j\omega_c n}$ as I did above, this causes the frequency to shift by $+\omega_c$   , moving what was centered at DC to now be centered at $\omega_c$. If you muiltiply the baseband signal with a sinusoid, you shift in two directions (moving the signal to the positive and negative frequency carrier) since a sinusoid contains two exponential frequencies as further explained below.

What you show in the lower plot is Large Carrier AM, and specifically such a plot could be created from adding a DC offset prior to doing the upconversion:

$$s[n] = (C + A e^{j\omega_c n})B \cos(\omega_m n)$$

Converting the real signal from your product to a complex signal can of course be done, but the result of this would be the expected two sidebands only that you now have in either the positive frequency axis or negative frequency axis. This is the Hilbert Transform that will convert your real signal that has complex conjugate symmetric positive and negative frequencies to a complex signal that only has positive or negative frequencies depending on the sign of Q.

My earlier answer regarding why the upper and lower sidebands are not equal magnitude in the FFT result

Either window your time domain data before taking the FFT or increase the length of your time domain data (and therefore your FFT). You are seeing the effects of “scalloping loss” which occurs when the actual frequency terms are not an bin centers.  Zero padding your time domain data before taking the FFT will also help to make the result look closer to what you would expect since that will interpolate more samples in the frequency domain resulting in balanced sum and difference frequency magnitudes.

Also the command fftshift() will move “DC” to the center of the axis to give you a plot with positive and negative frequencies which may be more intuitive. Keep in mind that a real signal will have positive and negative frequency components in the FFT result since each bin represents an exponential frequency term $e^{j\omega t}$, and as is clear in Euler's identity, a sinusoidal term is composed on two exponential frequency terms. For example:

$$2\cos(\omega t) = e^{j\omega t}+ e^{-j\omega t}$$

Your time domain repesentation would be a sine wave at 900 KHz and 1.1 MHz, so the result if you increase the time domain length by adding more samples or zero-padding will approximate the expected result showing two equally weighted frequency impulses in the positive frequency axis at 900 KHz and 1.1 MHz as well as two frequency impulses in the negative frequency axis at -900 KHz and -1.1 MHz.


Impact of the sample rate of the emitting signal : communication between HackRF One and USRP NI 2921 I am transmitting a 1 kHz cosine wave with a HackRF One. I receive it with a NI2921 USRP at a sample rate of 200kHz. All the process is on GNURadio Companion.

I changed the sample rate at emission to see the impact. I started at 200kHz and lowered of 1kHz each time. The sample rate at reception stays the same : 200kHz. Here are the results :

From 200kHz to 96kHz, nothing appeared. Then, at a sample rate of 95kHz at emission, the cosine spectrum showed two peaks at 13kHz and 87kHz, centered on 50kHz (DC component) giving a cosine of frequency 37kHz.
All the frequencies on the table above are the only ones showing something on the FFT of the received signal.

My questions are :

Why does the sample rate at emission impact the frequency seen at the receiver?
Why are there frequencies that do not make the cosine spectrum appear? (Only the frequencies present in the table work). For example, if I have the two same sample rates at emission and reception (200kHz for our example), it does not work.
Why are the frequencies seen on the FFT plot different from 1kHz, the real cosine frequency emitted with the HackRF One?

The HackRF doesn't support those sampling rates. From reading the firmware over reading libhackrf to reading the gr-osmocom hackrf_sink_c:

The set of supported sample rates is very limited:

osmosdr::meta_range_t hackrf_sink_c::get_sample_rates()
{
  osmosdr::meta_range_t range;

  /* we only add integer rates here because of better phase noise performance.
   * the user is allowed to request arbitrary (fractional) rates within these
   * boundaries. */

  range += osmosdr::range_t( 8e6 );
  range += osmosdr::range_t( 10e6 );
  range += osmosdr::range_t( 12.5e6 );
  range += osmosdr::range_t( 16e6 );
  range += osmosdr::range_t( 20e6 ); /* confirmed to work on fast machines */

  return range;
}

So, use 8 MS/s.

The hackrf is instructed to use a different sampling rate than you think it uses, which obviously "scales" the spectrum.


Shift vector function in MATLAB I am working on my own shift vector function that will be used later to compute the convolution of two signals. The function has to shift the vectors either left or right depending on the magnitude and direction of the shift value. My code is:

%This function shifts a vector either left or right

function sv = shiftvector(v, n) %Where v is the vector and
                                %n is the shift value
vnew = zeros(1,length(v));

for i = 1:length(v)-n

if n>0

vnew(i) = v(i+n);

elseif n<0

    vnew(i+1:end) = v(1:end-n);

end

end
sv = vnew;
end

The issue I have is, whenever I input a negative shift value, MATLAB gives an error. Need help!

Matlab does not accept lesser than $1$ indices. Other languages often accept a $0$-based index. More about Shifting Sequences Using MATLAB is available, including function circshift.m.


Variance of expectation estimator Given a set of $N$ iid samples, $x_k$, of random variable $X$ and the expectation estimator,
$$E[X] \approx \frac{1}{N} \sum_{k=1}^N x_k,$$
how do I find the variance of this estimator (not the variance of $X$) as a function of $N$?
I would only like to know the error in the estimate of $E[X]$ and what I would expect is that more samples (larger $N$) drives the variance down, so that
$$\lim_{N \to \infty} \bigg[\frac{1}{N} \sum_{k=1}^N x_k - E[X]\bigg] = 0.$$
Am I right about this? Perhaps variance is not the correct way to measure the error in the estimate? What I would like to understand is, how many $N$ is sufficient because there is a cost for every extra $x_k$ that is processed.

Without more knowledge on $X$, we cannot say a lot. However, if it is IID with mean $\mu$ and variance $\sigma^2$, the sample mean $\hat{\mu} = \frac{1}{N}\left(\sum_{n=0}^{N-1}x_n\right)$ is unbiased, since:
$$E[\hat{\mu}] = E\left(\frac{1}{N}\sum_{n=0}^{N-1}x_n\right) = \frac{1}{N}\sum_{n=0}^{N-1}E\left(x_n\right) = N\mu/N=\mu\,.$$
And its variance goes to zero when $N$ increases:
$$V[\hat{\mu}] = V\left(\frac{1}{N}\sum_{n=0}^{N-1}x_n\right) = \frac{1}{N^2}\sum_{n=0}^{N-1}V\left(x_n\right) = N\sigma^2/N^2=\sigma^2/N\,.$$
Thus, the expectation converges to the actual mean, and the variance of the estimator tends to zero as the number of samples grows. Under these definitions, the sample mean is a consistent estimator. Note that one could try to use other hypotheses: alternative norms, convergence in law, etc.


Encoding YUV Color Space with Positive Numbers I got 2 hdr rgb textures/signals. When fetched one gets bicubicly interpolated and the other bilinearly. So I was thinking I can store as YUV and store Y from each texture in the bicubic one and the distribution of the rest of the channels are not as important.

Well another restriction is that the hdr output only accept pos float values.

So what I'm wondering are there some other colour space transform similar to standard YUV but which range is positive only?

Thanks!

Why don't you use YUV with a shift?

To make my example simpler, let's use the Color Space YCoCg.
So the channels Co, Cg are in the range [-0.5, 0.5].

When you encode them, make sure you encode Co', Cg' where Co' = Co + 0.5 and Cg' = cG + 0.5.


MFSK vs FSK: Why BER decreases increasing the symbols My question is related to the performance of 2-FSK modulation scheme against M-FSK (M > 2).
From the SNR vs BER curve it seems that as the number of tones (frequencies) increases the BER decrease for the same SNR compared to 2-FSK(pics attached). (I assume maintaining the same symbol duration)

As the number of tones increases the bandwidth increases but i can't figure out why the performance of MFSK also increases.
At first glance i would say that the resilience to the noise should be the same as the 2-FSK.

Which is the benefit of using 2-FSK instead of M-FSK (M>2) given a fixed available bandwidth?

Peter

(I assume maintaining the same symbol duration)

nope! That is usually not the case, because it would be unfair – if you send the same number of symbols per second, you'd get twice as many bits per second using 4-FSK rather than 2-FSK.

Therefore, the comparison is typically done with identical Bit Energy rather than identical Symbol duration (which would be, for constant envelope signals, be equal to identical symbol energy).

Hence, with higher modulation order, you get less symbols per bit, and hence, for constant $E_b$ more $E_{sym}$.

Now, errors happen on a symbol level – and a Symbol Error Rate of $r$ doesn't imply a BER of $r$ – in fact, a wrong symbol has at most its worth in bits in wrong bits, but usually less (see: Gray code), so $\text{BER}\le\text{SER}$ generally holds, and one can in an uncorrelated noise scenario be almost certain that $\text{BER}<\text{SER}$.


Computing of speech signals modulation entropy I got a question regarding to modulation entropy of audio signal. There is some literature I found about this topic, http://www.irit.fr/recherches/SAMOVA/FeaturesExtraction.html#me.

I've read sample .wav file. How do I calculate entropy of this sampled signal of lenght N samples? How should I understand, in this case (in paper equation from above), propability of event?

Ok, I've found what i was looking for.

This paper describes how to measure entropy in signals:

Entropy Based Voice Activity Detection in Very Noisy Conditions, by Philippe Renevey and Andrzej Drygajlo. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.123.6098&rep=rep1&type=pdf

And this question covers the part of coding side of question:
https://stackoverflow.com/questions/18962990/calculate-entropy-of-a-file


Synchronizing Two Audio Tracks I would like to implement an algorithm that synchronizes two audio tracks.

The two tracks are very similar; you can imagine they are two recordings of sound (not necessarily music) from two microphones in the same room. This means they appear as one signal that passed through two different (linear) channels, with different additive noise.

I understand the first idea might be to look for the peak of their cross-correlation, but I am particularly interesting in finding literature on the topic, but all the papers I could find are either synchronizing two musical tracks (using chroma-based audio features), or synchronizing audio with a score. They also assume there might be time-warping, which is an unnecessary assumption in my case.

Chroma-based approaches are used when the tracks to align are two performances of the same material with widely different instrumentations - in which case you need to find the "lowest common denominator" which is harmony, melody, chord progression etc. Chroma are good for that - they totally abstract timbre, recording techniques. But you don't need those in your case.

As you have guessed, cross-correlation of the raw audio would be 1/ not very robust to noise, 2/ not very robust to differences in transducers (the two microphones might have very different responses) 3/ very costly if you have minutes of audio.

I suggest you to:

Extract a sequence of MFCC vectors from both signals. This way you'll have something of lower dimensionality and a bit more robust to noise and differences in transducers.
Optionally normalize the MFCC ("by design" coefficient 0 has more variance than coefficient 1 and so on...)... Ideally you'd want each row of the MFCC matrix to have variance 1.
Compute the cross-correlation of two matrices you get, along the time axis (that is to say, compute the 2D cross-correlation and just keep the data for the time axis).


Are Convolution and Deconvolution Kernels the Same? I need to clarify this and rather confused by this.

Lets say:

  $x = h * g$

  $x$ - measured data

  $h$ - raw data

  $g$ - instrument response function (convolution kernel)

So lets say I have $x$ and $h$. Can I say if I were to deconvolve $x$ with $h$ I will be able to obtain $g$?

Or if I have $x$ and $g$, deconvolving $x$ with $g$, I will be able to obtain $h$. So till here, is my understanding correct?

So can I say $g$ is both convolution and deconvolution kernel? Is this correct?

Question - When Are the Convolution Operator (Kernel) and the Deconvolution Operator (Kernel) the Same?

Discrete Convolution (Cyclic) is described by Circulant Matrix.
Circulant Matrices are diagonalized by the DFT Matrix (Will be denoted as $ F $).
So given a convolution kernel $ g $ with its matrix form given by $ G $ one could have its diagonalization by:

$$ G = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F $$

Where $ \mathcal{F} \left( g \right) $ is the Discrete Fourier Transform (DFT) of the kernel $ g $.

In the general case the deconvolution is given by:

$$ x = G h \Rightarrow {G}^{-1} x = {G}^{-1} G h = h $$

Where $ {G}^{-1} $ is the Deconvolution Operator (Inverse of $ G $ the convolution operator).
In our case we're looking for the case $ G = {G}^{-1} \Rightarrow G G = I $.

Utilizing the diagonalization of Discrete (Cyclic) Convolution Operators by the DFT Matrix (Which is Unitary Matrix):

$$ G G = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F = {F}^{H} \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) F $$

Namely we need $ \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) \operatorname{diag} \left( \mathcal{F} \left( g \right) \right) = I $ which means $ \forall i, \, {\mathcal{F} \left( g \right)}_{i} \in \left\{ -1, 1 \right\} $, namely the Fourier Transform of $ g $ has the values of -1 or 1.
So there is a set of kernels $ \mathcal{S} $ which is defined by:

$$ \mathcal{S} = \left\{ g \mid \forall i, \, {\mathcal{F} \left( g \right)}_{i} \in \left\{ -1, 1 \right\} \right\} $$

One should notice that the identity Kernel ($ I $) is indeed part of this set of kernels.

The above kernels obey $ \delta \left[ n \right] = \left( g \circ g \right) \left[ n \right] $ where $ \circ $ is the Circular Convolution (Which is multiplication of the DFT).

How could one generate such a kernel?
Well, just generate a vector of numbers which each is composed of $ 1 $ or $ -1 $ and use its ifft().
Pay attention that if you are after real operator the vector must obey the symmetry of the DFT for Real Numbers.

You may find a code to replicate the above in my StackExchange Signal Processing Q19646 GitHub Repository.


Training and testing set for speech recognition I am the moment looking for a dataset similar to the caltech 101, but in this case for different phrases, words and so on... In general just words.

Does someone know if such dataset exist, or I have to create it myself...

There are many open source datasets for speech available for free. You can find most at http://www.openslr.org

Beside that you can check http://voxforge.org


Using pathloss in matlab This is probably a very simple case but I am not sure where I am going wrong. In order to estimate the effect of path loss in Matlab, the function fspl(R,lambda) gives a pathloss value for a distance Rin $dB$.

My question is that, for using this function should one directly multiply the path loss value in $dB$ or convert it to linear scale and then multiply. If I use the $dB$ value, the signal power increases and if I use the linear scale value, there is a lot of difference between the transmitted and received power. Here's a snippet of my code:

SNR =2;
Fs=12e6;
Fc = 2.4e9;
d = 5; lambda = physconst('LightSpeed')/Fc;
data=randi([0 1],3000,1);
hMod = comm.OQPSKModulator('BitInput',true);
modulated_sgl = step(hMod, data);
pathloss_dB = fspl(d,lambda);
tx_sgl = modulated_sgl * pathloss_dB ;

If I try to plot the power spectral density of tx_sgl and modulated_sgl using the following:

hss = dsp.SpectrumAnalyzer('SampleRate', Fs);
step(hss,[tx_sgl modulated_sgl]);
release(hss);

I observe however that the power of tx_sgl is higher than the modulated_sgl. Can someone guide me where I may be going wrong.

The function gives you the path loss in dB. This number is positive, i.e. a path loss of $3$ dB corresponds to a gain of $-3$ dB. So for a given path loss $x$ (in dB), the signal is attenuated by a factor $10^{-x/20}$. Note that multiplication with a dB value hardly ever makes any sense.


Detect and isolate veins in the palm? I am trying to figure out a way to find out the vein pattern in the hand.. Please could someone suggest ways to do the same.? I know first i would need to find the ROI in the image and perform further processing .. Can gabor filter be used as one of the steps ?
Thanks

Ridge filtering seems like a possible processing. I have not tried to clean up the non-veins ridges; perhaps there would not be need for that step in more constrained acquisition conditions.

In Mathematica:

img = ColorConvert[Import["http://i.stack.imgur.com/5og5Z.jpg"], "Grayscale"];
r = RidgeFilter[ColorNegate@img, 3];
Binarize[ImageAdjust@r, .03] // Thinning


What kind of filter I should use to remove the oscillations in this autocorrelation function? I have an autocorrelation function which is shown as following,

I do believe these trailing oscillating wigs are spurious and should be removed by some kind of filter. But I am not familiar with any types of filters. Any help will be appreciated! The original data can be downloaded at

~

Update @ Apr 14 2017:

This result was processed with the filter suggested by Dan.

The following comb filter structure can be used to filter your data prior to computing the autocorrelation which will significantly reduce or eliminate the ringing due to the repetition in your signal. However be aware that it may also be filtering signal content based on the spectral occupation of your signal. The structure of the signal and the frequency response is shown below (for N=101).

Where N is the number of samples to be a 0.5 second delay. Implemented in Matlab using

out = filter([1 ones(1,N-1) -1],1,in)

This filter above as given is a highpass response so depending on the spectral content in the signal of interest, this may not be ideal.

Another option that will equally create nulls at the repetition rate and harmonics of the repetition rate is a moving average filter:

Implemented in Matlab using

out = filter([ones(1,N)],N,in)

Note in this case the filter was normalized, a 2 could be used as the second parameter in the first case to normalize that filter as well.

Note that the first filter removes DC and is a highpass function which may not be desirable, while the second filter has a 1/f frequency roll-off as a sinc function in addition to removing the harmonic nulls which will likely widen the autocorrelation response. If the high frequency content is desired, a third filter option that would be interesting to see can be done based on this post

Transfer function of second order notch filter

This filter also removes DC but can have a much tighter nulling response and flat elsewhere based on a nulling parameter $\alpha$

By inserting zeros in the frequency response we can get the base spectrum to repeat, so in your case you would insert N-1 zeros to perform the following filter function (see referenced post for structure of this filter):

out = filter((1+alpha)/2*[1 zeros(1,N-1) -1], [1 zeros(1,N-1) -alpha], in)

The frequency response for $N = 10$ and $\alpha = 0.9$ is shown below. (Your N of course is much larger). If you vary $\alpha$ less than or more than what I used up to but not equal to 1 (I recommend using $\alpha$ in the range between 0.7 and 0.999); the nulls get sharper as $\alpha$ approaches 1.

If you do implement these on your data, please post the results for all these cases so we can see the difference it has caused.


Quadrature component - condition on equivalent lowpass signal I have a question regarding the equivalent lowpass signal $x_{LP}$ of a bandpass signal $x(t)$ and the quadrature component.

I'd like to find a condition on $x_{LP}(t)$ and $X_{LP}(f)$ such that the quadrature component of $x(t)$ is zero. (for example $x(t)=a(t)cos(2 \pi f_c t)$).

I know that
$x(t)=Re\{(x(t)+j\hat{x}(t))e^{jw_0t}\}\\
= x(t)cos(w_0t)-\hat{x}(t)sin(w_0t)$

So my idea was that the hilbert transform of x(t) should be zero, in order that the signal has no quadrature component. That would mean that $X_{LP}(f)$ has even symmetry?

I hope that someone can help me with this question!

Your idea with the Hilbert transform doesn't work. The only signal (apart from $x(t)=0$) for which the Hilbert transform is zero, is a constant signal.

A band pass signal $x(t)$ can be written in terms of its complex envelope $x_{LP}(t)$:

$$\begin{align}x(t)&=\textrm{Re}\left\{x_{LP}(t)e^{j\omega_0t}\right\}\\&=\textrm{Re}\{x_{LP}(t)\}\cos(\omega_0t))-\textrm{Im}\{x_{LP}(t)\}\sin(\omega_0t))\tag{1}\end{align}$$

From $(1)$ it is clear that the quadrature component is zero if the imaginary part of the complex envelope $x_{LP}(t)$ is zero, i.e., if $x_{LP}(t)$ is real-valued. In the frequency domain, this is equivalent with $X_{LP}(f)$ being conjugate symmetric: $X_{LP}(f)=X_{LP}^*(-f)$

Note that your equation involving the Hilbert transform is wrong. The complex band pass signal

$$x_a(t)=x_{LP}(t)e^{j\omega_0t}\tag{2}$$

is called the analytic signal. Its real part equals the (real-valued) band pass signal $x(t)$, and its imaginary part is the Hilbert transform of its real part:

$$x_a(t)=x(t)+j\hat{x}(t)\tag{3}$$

where $\hat{}$ denotes the Hilbert transform. So we simply have

$$x(t)=\textrm{Re}\{x_a(t)\}\tag{4}$$

which is the same as Eq. $(1)$.


Is $ H\left\{x_h(t)\right\} = H\left\{x\left(t-\frac{\pi}{2w_c}\right)\right\} $ where $H$ is Hilbert transform Accordind to my understanding of the following:

$$H\left\{X(f)\right\} = -j\;sgn(f)\; X(f)  \qquad(1)$$

my work so far:

since $ x_h(t) = x\left(t-\frac{\pi}{2w_c}\right)$, Their Hilbert transforms are equal in time and frequency domains.Hence,

$$ H\left\{x_h(t)\right\} = H\left\{x\left(t-\frac{\pi}{2w_c}\right)\right\}, \:\text{where $H$ is Hilbert transform.
} $$

Let  $ \: y(t) = x\left(t-\frac{\pi}{2w_c}\right)\ \: $.We have

$$ H\left\{x_h(t)\right\} = H\left\{y(t)\right\} .$$ Therefore,

$$H\left\{X_h(f)\right\} = H\left\{Y(f)\right\} \\ \\ $$

L.H.S

by applying what is stated in the picure, we have:
$$H\left\{X_h(f)\right\} = -j\;sgn(f)\;  X_h(f)   = -j\;sgn(f)\;  \left[-j \; X(f) \; sgn(f)\right] $$

after simplification:
$$\text {L.H.S} = -X(f)$$

R.H.S
$$H\left\{Y(f)\right\} = -j\;sgn(f)\; Y(f) = -j\;sgn(f)\; \left[X(f) \; e^{-j2\pi \left(\frac{f}{4f_c}\right)} \right ] = \text {R.H.S} $$  by simplifying $$ \frac{\pi}{2w_c} = \frac{\pi}{2\left[2\pi f_c\right]}  = \frac{1}{4fc} $$ and by the shifting property:

$$F \left\{x(t-t_0)\right\} = X(f) \; e^{-j 2\pi f t_0}  \\ \\ $$

$$ - X(f) \quad \ne  \quad -j\;sgn(f)\; \left[X(f) \; e^{-j2\pi \left(\frac{f}{4f_c}\right)} \right ] $$

Thank you in advance.

It is quite unclear what you're actually asking, and most of your equations have no meaning. But one thing that you need to understand is that the equation

$$\mathcal{H}\{x(t)\}=x\left(t-t_0\right),\qquad t_0=\frac{\pi}{2\omega_0}=\frac{1}{4f_0}\tag{1}$$

only makes sense for sinusoidal functions $x(t)$ (of frequency $f_0$). For non-sinusoidal signals $x(t)$ the Hilbert transform of $x(t)$ can only be written as an integral:

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\tag{2}$$

where the improper integral is interpreted as the Cauchy principal value. In general, the Hilbert transform is not just a simple time shift.

For a sinusoidal signal $x(t)=\cos(2\pi f_0t+\phi)$ you can indeed show that the Hilbert transform corresponds to a delay of $t_0=1/(4f_0)$, i.e., in the frequency domain the following holds:

$$-j\,\textrm{sgn}(t)\,X(f)=X(f)e^{-j2\pi f t_0}=X(f)e^{-j\frac{\pi f}{2f_0}}\tag{3}$$


What is it called if a system is "in rest" until time $t=t_0$? I am a mathematics student, and am quite new to the whole control and system theory (haven't heard any lectures). I currently work with dynamic systems, but do not know the correct words for every attribute, yet.
I would like to describe a system as follows:

A system is in "rest" until time $t=t_0$ for some $t_0\in\mathbb{R}$ if the system response $y(t<t_0)$ and the state $x(t<t_0)$ are equivalent to $0$ for zero input $u(t<t_0) \equiv 0$.

I am quite sure, that there is a word describing this scenario.
I know the term causal in terms of systems, but I don't think this is the word I need.

EDIT: To maybe make it more clear what i want:
I want to say "A system is in XXXX, when for zero input there is zero output." (~the system does nothing before any input is applied)

There are two things to consider here. First of all, a non-causal system can react to input that hasn't arrived yet. And, secondly, there could be non-zero initial conditions that would contribute to the output signal without any input signal.

So what you probably mean is a "causal system initially at rest" or "a causal system with zero initial conditions".


What is the minimum phase shift for Quadrature Amplitude Modulation (QAM)? I'm a bit of a noob here, looking for a simple formula or just a table which can tell me the absolute minimum phase shift angle possible for different variations of QAM. Particularly QAM 16, 64, 256 and 1024. If the answer is zero (is it?), then what is the minimum phase shift possible in each case other than 0?

Additionally, what is the limiting factor which governs the rate at QAM can switch between different phases? In other words, what is the maximum rate at which QAM can switch between different phases, and how is that determined for a particular device? Is there a relation to frequency of transmission (say a 1 GHz signal), or are frequency of transmission and phase-switching rate completely independent?

The minimum phase shift in 16-QAM, 64-QAM, 256-QAM etc is 0 (same value regardless of whether you are using degrees or radians!) and occurs when the symbol transition is from one constellation point on some radius vector to another constellation point on the same radius vector.  If we take the constellation points as having coordinates $(i,j)$ where $(i,j) \in \{-3, -1, +1, +3\}$ in 16-QAM, then $(1,1)$ and $(3,3)$ lie on the same radius vector (line of slope 1 through the origin) and so the RF phase remains unchanged in the transition from one constellation point to the other. Obviously, the same example will also work for 64-QAM, 256-QAM, and 1024-QAM.

What is the minimum non-zero phase shift?

Well, for 16-QAM, the constellation points $(1,3)$ and $(3,3)$ are distance 2 apart, and together with the origin, form an obtuse triangle with sides of lengths $2$, $\sqrt{10}$, and $\sqrt{18}$. From the standard formula
$$c^2 = a^2 + b^2 - 2ab\cos C$$ relating the sides $a$, $b$, $c$ of a triangle opposite to angles $A$, $B$ and $C$, we get that
$$2^2 = 10 + 18 - 2\sqrt{180} \cos C$$
giving $\arccos \frac{2}{\sqrt{5}}$ as the minimum nonzero phase shift in 16-QAM.
Similar calculations for 64-QAM, 256-QAM and 1024-QAM as left as an exercise for the reader.


Confusions regarding frame rate? We know that when many image are shown in succession they create illusion of video

But what is the minimum rate of frame/images so that still images appear as video?
In this case, images/frames are moving horizontally or vertically?

In this case, images/frames are moving horizontally or vertically?

Neither. They just appear; usually, they stand still for as long as a frame is displayed.

  But what is the minimum rate of frame/images so that still images appear as video?

That depends on a lot of factors like brightness, amount of motion in the material, and mental state of the observer.

Typically, upwards of roughly 20 frames a second, things are perceived as fluid motion, but there seems to be some benefits to going higher in frame rate, especially in high-movement scenes.

Traditionally, cinema and TV used to be 24 or 25 frames a second, but this number has changed in recent years, with digital TV producing whatever frame rate the video codec used does.


How to proceess Image with curved corners and edges to give sharp Corners? I have several Images of shape H, F, I, L.. which have curved corners. How do I process them to give sharp right-angled corners? I tried approximating them with inner rectangles, but this takes a LOT of time.. How do I proceed with it?

This is the sample Image. I want to process it such that I get rid of rounded corneres.

here's how houghLinesP would look like. Not a solution yet, but maybe you can work with that:

    int main()
    {

        cv::Mat image = cv::imread("SharpenCorners.jpg");

        cv::imshow("image", image);

        cv::Mat gray;
        cv::cvtColor(image,gray,CV_BGR2GRAY);

        cv::imshow("gray", gray);

        cv::Mat output; image.copyTo(output);

        // gradient:
        cv::Mat gx, gy;
        cv::Sobel(gray, gx, CV_32F, 1, 0, 3 );
        cv::Sobel(gray, gy, CV_32F, 0, 1, 3 );
        cv::Mat gm;
        cv::magnitude(gx,gy,gm);

        //cv::Mat g_thres = gray < 100;
        cv::Mat g_thres = gm > 100;

        //int lineThres = 5;    // 1st image
        //int lineThres = 10;// 2nd image
        //int lineThres = 25;// 3rd image
        int lineThres = 50;// 4th image

        std::vector<cv::Vec4i> lines;
        cv::HoughLinesP( g_thres, lines, 1, CV_PI/(4*180.0), lineThres, 10, 10 );

        for( size_t i = 0; i < lines.size(); i++ )
        {
            cv::line( output, cv::Point(lines[i][0], lines[i][1]),
                    cv::Point(lines[i][2], lines[i][3]), cv::Scalar(155,255,155), 2, 8 );
        }

        cv::imwrite("houghpSharpen4.png", output);
        cv::imshow("g thres", g_thres);
        cv::namedWindow("output"); cv::imshow("output", output);
        cv::waitKey(-1);

        return 0;
    }

with these results for given thresholds:

1st:

2nd:

3rd:

4th:

Additional idea:

Giving your own comment, that you might try to approximate each contour with rectangles, you could start with houghLinesP results (with a pretty high threshold to get at least all "long" line elements) and try to create the rectangles based on those lines


Random Signals - statistical properties are time dependant? I'm taking a course on DSP and we're being introduced to the random signals, in particular continuous time and discrete time random signals.

We're told that if we repeat a single random experiment at each point in time, the random variable associated with each random experiment when collected together will form the random signal.
We are then told that the random signal/process has statistical properties such as mean that are time-dependant.

Now that makes no sense. Surely, if you're doing the same experiment, the mean of the entire random signal will be the same as the mean of the random variable for one experiment? How will time affect the random signal?

The answer to the question (a counterexample)

Properties of random processes will in general be time-dependent. They are not only when talking about stationary processes. Another related concept (not relevant here but that you might find interesting) is ergodicity. Most of us find it difficult to understand the difference between both concepts, but then you find this brilliant Dilip's answer to clear out all doubts.

Going back to your question: take for example the simple random process

$$X(t)=At$$

where $A$ is some random variable with mean $\mu_a$. The mean of the process would then be:

$$\mathbb{E}[X(t)]=\mathbb{E}[At]=t\mathbb{E}[A]=t\mu_a$$

It's easy to see that the mean depends on $t$, thus it's time-dependent. I hope this easy example shows why, in general, properties of random processes vary with time.

Additional info that might help

I think your confusion arises from thinking that to build a stochastic process you grab many realizations of a single random variable and put them in a sequence to represent a time-signal. That's not how a random process is built in general. If you do this, then the mean will not vary with time because the process will be stationary. But if you don't, it may or may not vary.

I'll use an example here to make myself clear. Suppose you have a random variable $Z$ that represents the result of throwing a fair 6-faced dice. Then $\mu_z = 3.5$. If you throw the dice 3 times and you put them together to form a discrete-time signal of 3 samples, then you will have the following random process:

$$X(n)=Z_1\delta(n) + Z_2\delta(n-1) + Z_3\delta(n-2)$$

where $Z_i$ represents the random variable corresponding to the $i$-th throw of the dice.

Note that each sample of $X(n)$ is a random variable. This is correct for any random process in general. For a given random process $X(n)$, $X(n_0)$ for any $n_0$ will be a random variable.

On the other hand, once you've thrown the dice 3 times, you are not dealing with anything random anymore: you have a signal. That is, a realization of the random process.

If you build a random process such as we have just built $X(n)$, then you are right: the mean of the process will be constant and equal to the mean of the random variable used for each sample (as in this case). In fact, all properties will be constant and equal to those of the variable. But that's just in that specific case, because all samples of the process will be i.i.d., and thus, as @Fat32 has explained in the comments, the random process will be independent of time ($t$). In our dice example:

$$\begin{align}
\mathbb{E}[X(n)]&=\mathbb{E}[Z_1\delta(n) + Z_2\delta(n-1) + Z_3\delta(n-2)]\\
&=\mathbb{E}[Z_1]\delta(n) + \mathbb{E}[Z_2]\delta(n-1) + \mathbb{E}[Z_3]\delta(n-2)] \\
& = \mu_z \qquad \text{for any} \  n\in\{0,1,2\}
\end{align}$$

But this will not always be the case with random processes, such as in the example I wrote at the top of the answer, in which every random variable obtained by evaluating the random process at different values of $t$ is different from all others, therefore preventing the process from being i.i.d. and making its statistical properties time-dependent.


How to automatically detect blurred document images The question I have is related to image quality evaluation. Suppose you are supposed to perform OCR (Optical Character Recognition) on a set of document images. Some images, however, are very blurring, and for these images blind deconvolution is needed. However, there are good images that do not need deconvolution. Then my question is: could I find a criterion to automatically separate these images into two groups: blurred image and clear images and then treat them differently when performing OCR?

Yes, but doing so is difficult. You are basically assuming that when a certain feature is sufficiently present in an image, you can apply a filter and realize an improvement in OCR accuracy. This implies that you can find the right feature, apply the right filter, and evaluate OCR accuracy correctly. Since there are lots of features that might be relevant (edge prominence, high-frequency energy, PSF estimation, unfiltered image OCR accuracy metrics (as @Matt mentions), etc.), lots of filters available (blind deconvolution, among others), and you'll need lots of ground truth images to prove system effectiveness, you are looking at a lot of work and experimentation to create a system that works well.

You also need to make sure you are not working against the OCR software itself, which may do its own image processing to improve its accuracy. In other words, make sure that artifacts introduced by your deblurring do not make OCR accuracy worse, even if the deblurring was done "well" from a deblurring perspective.

All that being said, it should be possible to do this. You might even be able to do this adaptively, by automatically and iteratively adjusting parameters and approaches to optimize the result.


New to DSP and having trouble filtering data I am trying to make program in matlab that counts the number of times the accelerometer was lifted up and down (a rep counter for exercise). Basically my phone is getting accelerometer data on all three axes and then i import this data into matlab. I have tried to filter the data using a simple filter a = [1 -0.9] b=0.1 filter(b,a,data) but that doesn't really get rid of all the high frequency noise.

Another problem i am running into is that when i apply certain filters to the data the data gets shifted so i can no longer integrate it to get velocity (my original plan was to simple filter out high frequencies and then integrate the data to get velocity and then count the number of times it crosses zero and divide by two (which would correspond to the number of times it was at a peak and trough - divided by two for number of cycles).

But the data is so noisy that when i apply a filter not only does it get shifted away from the x axis making it impossible to integrate it, it also still has noise which would give me incorrect zero-crossings.

Is there a good beginner introduction to signal processing and accelerometer data processing? What techniques should i be aware of?

If you want to try to give the data a go here is some sample data, the rep count should be 5, The data was all taken at a 50Hz sample rate. I have to put data in the comments because SE doesnt allow posting multiple links a new poster

You should try with implementing your own linear filter. It is not really difficult, and you will have more control about what's going on behind the scene and it will make it easier for you to tune the filter bandwidth.
Speaking of bandwidth, you can try to plot the frequency analysis of your data using Fourier transform, that will help you in finding a suitable cutoff frequency.
A linear filter implements the equations $\sum_{i=1}^n a_i y[i] = \sum_{j=1}^m b_i x[i]$. As long as you normalize the weights $(a_i), (b_i)$ to sum to 1, you should preserve the range of the data.


Frequency domain vs time domain equalization It is said that a frequency domain equalizer is equivalent to an infinite length time domain equalizer, why is that?

And if this is the case, why hasn't frequency domain equalization been used earlier?

Has time domain equalization any advantages when compared with frequency domain equalization?

Frequency domain equalization requires a certain block structure to be efficient. In OFDM, adding a cyclic prefix to the transmitted block, and using IDFT at the transmitter and DFT at the receiver allows the use of one-tap frequency equalizer. In the time domain, the equalizer length is proportional to the number of channel taps. We can use frequency domain equalization in single carrier systems as well, again with proper block structure, and operations.


Series connection of digital PT1 filters I was trying design $N^{th}$ order filter, I was suggested a series connection of low-pass PT1 filter. Basically, the filter coefficients are like this:

$$a_i = \dbinom{n}{i} T^{i}$$

This works for $N^{th}$ order, while $a_1, a_2, a_3....a_n$ being the filter coefficient.

I don't know how this works, in theory. I tried searching, but could not come up with anything. if anyone has any idea about this. I want to know how this works.

You could have mentioned where you got that formula from. But anyway, I'll try to reconstruct it.

A first order discrete-time low pass filter has the following transfer function:

$$H(z)=\frac{1+c}{1+cz^{-1}},\qquad -1<c<0\tag{1}$$

where the constant $c$ determines the cut-off frequency. If you cascade $N$ of these filters, the total transfer function is

$$H_N(z)=\prod_{k=1}^NH(z)=\frac{(1+c)^N}{(1+cz^{-1})^N}\tag{2}$$

Using the Binomial Theorem this can be written as

$$H_N(z)=\frac{(1+c)^N}{\displaystyle\sum_{k=0}^N\dbinom{N}{k}c^kz^{-k}}=\frac{(1+c)^N}{\displaystyle\sum_{k=0}^Na_kz^{-k}}\tag{3}$$

with denominator coefficients

$$a_k=\dbinom{N}{k}c^k\tag{4}$$

Note that the numerator $(1+c)^N$ is just there to make sure that the scaling is correct, i.e., that the transfer function equals $1$ at $z=1$ (DC).

Following MBaz's suggestion, I add the relationship between the constant $c$ in $(1)$ and the filter's $3\,\text{dB}$ cut-off frequency $\omega_c$. For this we have to solve

$$|H(e^{j\omega_c})|^2=\frac12|H(1)|^2=\frac12\tag{5}$$

From $(1)$ we get

$$\frac{(1+c)^2}{1+2c\cdot \cos(\omega_c)+c^2}=\frac12\tag{6}$$

from which we obtain after some elementary algebra

$$c=\cos(\omega_c)-2+\sqrt{\cos^2(\omega_c)-4\cos(\omega_c)+3}\tag{7}$$


Matching input buffer size to opus frame size I am working on an audio effect that passes it's input though the opus codec. I want the effect to be a VST, so I don't have control of the input buffer size and the sample rate. However, opus requires a sample rate of 48Khz and frames of 120, 240, 480, 960, 1920, or 2880 samples. What is the best way to convert from an arbitrary frame size and sample rate to a known one? A loss of fidelity is OK.
Here is the signal chain as I understand it:

Audio comes from a VST host application such as Ableton Live.
I'm using JUCE as a library for creating the plugin so the audio is then past to the internals of JUCE.
JUCE calls the processBlock method on my plugin with an audio buffer. When testing I can force this audio buffer to be a certain size, but I want to support hosts with arbitrary sized buffers.
I need to somehow convert from the host buffer size, back to a frame size supported by opus. This is where my question resides.
I call opus_encode_float on the buffer to compress the data. This introduces sonic artifacts I want to keep.
I call opus_decode_float on the compressed data. This gives me back a buffer the size of a valid opus frame.
I need to then convert this back to a size the host is expecting (same size as the input buffer).

I ended up solving the problem. Because the input buffer was smaller than 1 frame, I was able to write the input to a circular buffer, wait until I have enough data, then do the processing. The result was the written to a second circular buffer, that would be read every block.


Modelling Unwanted Signal in a LMS Adaptive Filter I'm having some confusion learning about the LMS Adaptive Filter. I understand that the whole model of adaptive filters relies on the fact that we give it a reference signal to which it keeps comparing the input * filter with and the filter coefficients keep changing until the error between input and reference is zero.

Let's say I have a telephone conversation with a time-varying sinusoid added on top that changes frequencies every so often. How exactly do I give it a reference signal? I mean, if I had a reference signal of what I wanted the input to become, I would have just used it? And similarily if I had a reference signal for my unwanted sinusoid, I could have just subtracted it from my original sound.

What am I missing here?

The LMS and many of the variants of Adaptive Filters (In the Linear System context) work in the following settings (Intuitive):

You have access to 2 signals.
One signal is the result of the other one when a Linear System is applied.

This sounds really limiting, yet in practice it is powerful and flexible.

In the settings you mentioned the most known and similar problem is the Echo Cancellation model.

Pay attention that your model can't be formed as a Linear System.
Namely the connection between the clean signal and the corrupted signal can't be described using Linear System.


necessity of patch size in edge detection algorithms I have recently started learning about the image edge detection methods. I studied few of these methods including Sobel and Canny. I observed that edge detection is performed using patch/kernel size as 3x3.

My question is, Is it necessary to fix the patch size as 3x3? OR can we take it as 5x5, 7x7, 9x9 etc.

I have already studied few technical papers that says that patch size is taken as low as possible in order to make computational complexity as low as possible

Let's start to given some information.

We now that the sobel operator perform a aproximmation of the image gradient using convolution ( $ * $ ).

Then,

Being $\mathbf{A}$ an image of $n \times m$ size.

The operator $\mathbf{K}$, commonly, uses two $3 \times 3$ kernels which are convolved with the original image. Thus we have

$\mathbf{G_x} = \mathbf{K_x} * \mathbf{A} $

$\mathbf{G_y} = \mathbf{K_y} * \mathbf{A} $

Being
$$
  \mathbf{K_x} =
  \left[ {\begin{array}{cc}
                -1 & 0 & +1\\
                -2 & 0 & +2\\
                -1 & 0 & +1\\
              \end{array} } \right]
$$

And,

$$
  \mathbf{K_y} =
  \left[ {\begin{array}{cc}
                -1 & -2 & -1\\
                 0 & 0 & 0\\
                +1 & +2 & +1\\
              \end{array} } \right]
$$

So, an extension of this kernel is a $5 \times 5$ matrix

$$
  \mathbf{K_y} =
  \left[ {\begin{array}{cc}
                -1 & -4 & -6  & -4 & -1\\
                -2 & -8 & -12 & -8 & -2\\
                 0 &  0 &  0  &  0 &  0\\
                +2 & +8 & +12 &  8 &  2\\
                +1 & +4 & +6  & +4 & +1\\
              \end{array} } \right]
$$

$$
  \mathbf{K_x} =
  \left[ {\begin{array}{cc}
                 +1 &  +2  & 0 & -2  & -1\\
                 +4 &  +8  & 0 & -8  & -4\\
                 +6 &  +12 & 0 & -12 & -6\\
                 +4 &  +8  & 0 & -8  & -4\\
                 +1 &  +2  & 0 & -2  & -1\\
              \end{array} } \right]
$$

This article demonstrates how to get Sobel gradient operators analytically. For that was using linear approximation of brightness in window 3x3.  Every pixel in window has its own weight. Special weight matrix selection gives Sobel gradient operator.


Why is this recursive DFT algorithm not equivalent to this iterative method? Edit 1/30 - Taking @Fat32's edits into account, it seems like there is still an issue with the scale of the frequency axis. While version 2 correctly identifies the response at 1 HZ, version 1 seems to be covering a frequency range of 1/T.

I'd like version 1 to cover the same range of frequencies as version 2. Can anyone help with this?

(End edit 1/30)

--

I have 2 DFT algorithms that I would expect to return same frequency response. Unfortunately, the more efficient recursive algorithm does not seem to behave as expected.

Version 1 - Efficient Divide and Conquer method, but output is incorrect:

function [ s ] = DFT_ver_1( x )

N = length(x);

if N == 1
    s = x(1);
    return;
else

X1 = DFT_ver_1(x(1:2:N));
X2 = DFT_ver_1(x(2:2:N));
s = zeros(N,1);

for k = 1:(N/2)
        W(k)     = exp(-1j*2*pi*(k + 1)/N);
        s(k)     = X1(k) + W(k) * X2(k);
        s(k+N/2) = X1(k) - W(k) * X2(k);
end

Version 2 - Inefficient, but generates expected results:

function [ s ] = DFT_ver_2( x , f, n)

N = length(x);

b = zeros(N,N);
b_1 = -2*pi*(f/N);

for idx = 1 : N
for jdx = 1 : N
    b(jdx, idx) = b_1(idx) * n(jdx);
end
end

anlyz_fn = cos(b) + 1i * sin(b);
s = single(zeros(N,1));

for k = 1 : N
    val = 2 * sum(anlyz_fn(k,:) * x.');
    s(k) = val;
end

end

I would expect both algorithms to return the same result. However, this does not appear to be the case.

You can see that the bottom plot (version 2) shows the correct frequency response to a 1 HZ signal, but the middle (version 1) does not:

Code to run both:

       T  = single(3.20); % (sec) time window
       dt = single(0.05); % (sec) sample time
        N = T / dt;          % (ct) num samples
        n = 0 : N - 1;       % (ct) bucket index vector
        t = 0 : dt : T - dt; % (sec) time vector
smpl_freq = 1 / dt;          % (Hz) sampling frequency
 freq_res = smpl_freq / N;   % (Hz) frequency resolution
        f = n * freq_res;    % (Hz) freq. bucket vector

x_t = sin( 2*pi*t );

subplot(3,1,1);
plot(t, x_t);
xlabel('Sample Time (sec)');
ylabel('Magnitude');

subplot(3,1,2);

s_v1     = DFT_ver_1( x_t); %, f, n );
s_v1_mag = abs( s_v1(1:N/2) )/N;

stem(f(1:N/2)/T, s_v1_mag);

xlabel('Freq. Response (Hz) v1');
ylabel('Magnitude');

s_v2     = DFT_ver_2( x_t , f, n );
s_v2_mag = abs( s_v2(1:N/2) )/N;

subplot(3,1,3);
stem(f(1:N/2) / T , s_v2_mag);
xlabel('Freq. Response (Hz) v2');
ylabel('Magnitude');

All right, you made two small mistakes and here is the corrected version: (I changed your function's name). Btw, you can get rid of the for loop too... ;-)

function [ s ] = recDFT2( x )
%
% a recursive DFT
%

N = length(x);

if N == 1
    s = x(1);
    return;
else

    X1 = recDFT2(x(1:2:N));
    X2 = recDFT2(x(2:2:N));
    s = zeros(N,1);

    for k = 1:(N/2)
        W(k)     = exp(-1j*2*pi*(k - 1)/N);
        s(k)     = X1(k) + W(k) .* X2(k);
        s(k+N/2) = X1(k) - W(k) .* X2(k);
    end
end


Ideal high pass filter for discrete signal there. I currently get stuck on a question. I was asking to find an inverse discrete-time Fourier transform for the ideal high pass filter.
Here is the question
It is getting more confused after I read the second solution of this post
How to produce a high-pass filter from a low-pass one?
It tends out I can either calculate it by time-shifting property

Or I can calculate it through the definition

But they lead to a different result. My textbook solution said the second one is correct. But I cannot find any problem with the first one.
If we assume $\ 0 \le n < N-1$ for $h[n]$, using N = 61, $\alpha =30$, cutoff analog frequency $=1500Hz$.
(Assume the sampling rate $=44100 Hz$ )
I got these two graphs:
(for first solution)
and

(for the second)
Here are my codes

I really don't know why they lead to different results.

HINT: Your original LPF is $-\omega_c \lt \omega \lt +\omega_c$. When you shift it by $\pi$, your integral limit for second method is $\int_{-\pi}^{-\pi + \omega_c}$ and $\int_{\pi - \omega_c}^{\pi}$.

For the first method, it should be
$$
h[n] = (-1)^{(n-\alpha)}\frac{sin(\omega_c(n-\alpha))}{\pi(n-\alpha)}
$$
because you have to apply shift of $\alpha$ on $n$. This will fix your plot of (1).

For (2), with the changed limits as I hinted at the beginning
$$
h[n] = \frac{1}{2\pi}\int_{-\pi}^{-\pi + \omega_c}e^{-j\alpha\omega}e^{j\omega n}d\omega + \frac{1}{2\pi}\int_{\pi-\omega_c}^{\pi}e^{-j\alpha\omega}e^{j\omega n}d\omega\\
= \frac{1}{2\pi j(n-\alpha)}[e^{-j(n-\alpha)(-\pi+\omega_c)} - e^{j(n-\alpha)(-\pi)}] +\frac{1}{2\pi j(n-\alpha)}[e^{-j(n-\alpha)(\pi)} - e^{-j(n-\alpha)(\pi-\omega_c)}] \\
= \frac{1}{2\pi j(n-\alpha)}e^{-j(n-\alpha)(-\pi)}[e^{-j(n-\alpha)\omega_c} - 1]+\frac{1}{2\pi j(n-\alpha)}e^{-j(n-\alpha)(\pi)}[1 - e^{j(n-\alpha)\omega_c}]\\
= \frac{1}{2\pi j(n-\alpha)}(-1)^{(n-\alpha)}[e^{-j(n-\alpha)\omega_c} - e^{j(n-\alpha)\omega_c}]\\
= (-1)^{(n-\alpha)}\frac{sin(\omega_c(n-\alpha)}{\pi(n-\alpha)}
$$


C implementation of LP filter I have the following code, which I am told is a low-pass Butterworth filter. How do I derive the frequency response and phase for this? I only have experience in analog electronics and so cannot relate this to a real filter response. I'm not sure why everything is scaled up by 2^ 16 - it must be to do with maintaining adequate resolution during the calculations?

#define FILTER_TC(FC, FS) ((411775 * FC) / FS)'

where TC = (2^16) * 2pi * Fc * Fs, with Fc the cut-off frequency in Hz, and Fs the sample frequency in Hz.

#define FILTER_LP(FILTER_OUT, FILTER_IN, FC, FS) \
    (FILTER_OUT += ((signed long)FILTER_TC(FC, FS) * \
    ((signed long)FILTER_IN - (signed long)(FILTER_OUT)) >> 16)'

Thank you for any help with this.

jmv

You can figure out the filter's transfer function by rewriting the body of the preprocessor macro a bit. Let FILTER_TC(FC, FS) (which really just defines the cutoff frequency of the filter) be equal to k in the below expression.

#define FILTER_LP(FILTER_OUT, FILTER_IN, FC, FS) \
    FILTER_OUT = FILTER_OUT + ((signed long) k * ((signed long) FILTER_IN - (signed long) FILTER_OUT)) >> 16)

In the above, I just expanded the += operator to be a bit less terse. By inspection, we can write the discrete-time difference equation:

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

where $y[n]$ and $x[n]$ are equal to FILTER_OUT and FILTER_IN at time index $n$, respectively. This maps to a $z$-domain transfer function of:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{k}{1 - (1-k)z^{-1}}
$$

This is a lowpass filter, although it is not a Butterworth filter. Its frequency response is plotted below for k = 0.25.

In the above frequency response plot, the blue trace shows the magnitude response and the green trace corresponds to the phase response. It shows that the example you gave isn't a particularly good lowpass filter, even for a first-order filter. You could do better with respect to high-frequency attenuation by using a correct Butterworth structure, for example.

The bit manipulations and explicit casts to signed long types are used because the filter is implemented using fixed-point arithmetic. This is a scheme by which you use integers to represent real numbers using a fixed quantization level. Fixed-point implementations of DSP algorithms were almost ubiquitous in the past, but floating-point processing is becoming more and more common. I would expect to find fixed-point implementations primarily in low-power or embedded contexts today.


Noise power and bandwidth $$N=kTB$$
In this formula for calculating noise power, why do we use the total bandwidth of the system?

I mean our signal is at $f_0$ carrier frequency, why is it effected by a large bandwidth and not just only signal’s frequency bandwidth

$N = kTB$ is the expression for the total thermal noise power within bandwidth $B$. $kT$ (where $k$ is Boltzmann's constant and $T$ is the temperature in Kelvin) is the noise density, meaning power/Hz. It is also a white noise process thus this single constant (assuming a fixed temperature) is the noise density at all frequencies. ($kT$ at room temperature is approximately -174 dBm/Hz, a number worthy of committing to memory).

If you want to compute the total power within a bandwidth $B$, you simply multiply by $B$ to get $kTB$.

The second part of your question is not clear; typically we would compute the noise power within the signal's bandwidth in the process of estimating the signal to noise ratio as the well designed receiver would have filtered out the (amplified) noise power beyond the signal's frequency bandwidth. So $B$ in this case should be the signal's frequency bandwidth and not an arbitrarily "large bandwidth".


parametric eq with FIR filters Are there any good published approach for designing parametric EQs as FIR filters.  The standard techniques discussed in DSP textbooks are frequency-domain-sampling or using Parks-McClellan, which this involves passband, stopband, transition band, etc.  This is nice for rejecting a certain portion of the frequency spectrum, or mimicking and existing analog/digital filter.  However what if I want to design a filter directly as an FIR filter, for example a peaking/notching/shelving/low-pass/high-pass filter which takes frequency, gain, and bandwidth parameters?  Is there a standard technique for designing parametric EQs directly as FIR filters?

These filter types are well covered here http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

However these are all single biquad IIR implementations. To turn this into an FIR, you could simply calculate the impulse response of these filters and window to a finite size with the desired amount of accuracy. Another option is to design the filter as IIR, sample the frequency response and then do any of the many FIR fitting techniques.

In general you will find that FIR filters are not well suited for this and the number of taps required will heavily depend on the relationship between the corner frequency and the sample rate and also on the Q of the filter.


CT Slices to Planar XRay I have axial slices of chest CT scans. Now, I want to use all these slices to construct a planar structure like an X-ray with a coronal view (Posterior to Anterior or Anterior to Posterior view).
I have no idea how to proceed with this problem. One way I think is a weighted average of the slices along the coronal plane with more weights to the frontal slices for AP view.
Please share your ideas on how to proceed with this problem. Thanks in advance.

I think that what you are looking for is a maximum intensity projection (MIP).
The general approach has two parts:

Define the view
Obtain the projection

To define the view, you need to either rotate the data in-place or rotate the view plane around the data and then use the coefficients of the plane to read through your data matrix. In either case, you are going to need interpolation.
The most straightforward way to do this is to rotate the data in-place.
If you are interested in obtaining just the three principal views, you don't even need to rotate the data.
Let's assume that you have a matrix like $(y,x,z)$. That is, each axial image's data is in $(y,x)$ and each slice's axial image can be accessed via $z$.
Since these are the axial images, looking at an individual $(y,x)$ image is like looking through our volume from the top. So, we are looking at $(y,x)$ down through $z$.
To obtain an "X-Ray-like" image of the volume as is, all you have to do is iterate through each $(y,x)$ pixel, get all the pixel values that "live" in that column (so, get all the pixel values at location $(y,x)$ from slice $0,1,2,3,..,N$) and select the maximum value from that column. Then move to the next pixel and so on.
This will give you the MIP in the axial direction. To obtain the MIP from another point of view, all you have to do is to rotate your volume to align it with the view and run the same iteration through the image again.
In the special case where you only need the three principal directions, you don't even need to rotate the data, you just change the way you are reading it: Where before we were looking at $(y,x)$ down through $z$ for the axial view, we now need to look at $(y,z)$ down through $x$ for the sagital view...and so on.
Here is a quick example in Python (but generalisable easily to other platforms too):
import numpy
from matplotlib import pyplot as plt

# We need some data, so let's build a phantom.
# This is a sphere that is harded in the middle and contains
# a rectangular block.

v = numpy.zeros((10,10,10)) # This is the size of our phantom

for m in range(0,10):
    for n in range(0,10):
        for k in range(0,10):
            # This is basically the a reversed distance plus a dense rectangular block running through it.
            v[m,n,k]= 20-numpy.sqrt((m-5)**2+(n-5)**2+(k-5)**2) + (10 if m>3 and n>3 and k>2 and m<6 and n<6 and k<8 else 0)

# To get the axial MIP:

z = numpy.zeros((10,10)) # Notice this is a single image, NOT a stack

for m in range(0,10):
    for n in range(0,10):
        z[m,n]=numpy.max(v[m,n,:])

plt.imshow(z);plt.show()

And this will give you something like:

This is a MIP of the axial view, so looking at our volume from the top.
To look at our volume "from the side" (the sagital view), we just change the way we read the volume:
# To get a sagital view MIP:

for m in range(0,10):
    for n in range(0,10):
        z[m,n]=numpy.max(v[m,:,n])

plt.imshow(z);plt.show()

And this will give you something like:

Notice here that the $x$ axis of this image is the $z$ axis of the volume.
To obtain any other view, simply 3D rotate the data in your v to align it to the view plane and then run the z iteration again for the MIP.
You can see here that if the pixel size and slice thickness are not equal, then the sagital view will appear distorted.
If you are looking for a library to achieve this, you can use something like vtk, although, be aware that it requires some level of prior experience to use it and more overhead to get an example up and running.
Hope this helps.


Removing Bounded Noise from an Audio Signal I'm currently trying to remove some noise from a "contaminated" audio file. From what's given, I know that the noise is isolated to a specific frequency band - that is, from ~1400 Hz - ~1800 Hz.

The most straightforward approach is to obviously implement a stop-band filter - which I've done successfully. However, I wanted to explore alternative options being that a band-pass filter also eliminates the desired signal.

Here's what I've done thus far:

In MATLAB, I explored some aspects of the signal, such as it's power spectrum, autocorrelation, etc:

It's clear that the noise is concentrated in the band I mentioned above, and it appears that the noise is Gaussian.

My initial thought was to apply a spectral subtraction method, such as the Ephraim-Malah Algorithm, but one of the assumptions made by this method is that the underlying noise is Gaussian White Noise. In the file I have, the noise is limited to only a particular band, not the entire signal. My concern is that if I were to apply the algorithm, I would unintentionally introduce minor distortions to potentially good portions of the signal.

At the moment, I'm considering splitting the signal into two portions: use a stop-band filter to extract the unaffected frequencies, and use a band-pass filter to capture the noisy region, apply the algorithm to the latter band, and recombine the signals.

Is this an appropriate approach, or is there a simpler and/or better method that I could use?

NOTE:

Being that I have an idea of what the noise is, I could potentially use a Weiner Filter.

If you have some prior information about your noise e.g frequency band, you can simply implement an noise removing algorithm with emphasis on removing the noise on that band (I mean use a somehow weighted noise removing process).

One of nice methods is weighted multi-band spectral subtraction (this is not the only method and in fact these weighted methods have many variations).
The improvement of MBSS is due to the fact that the multi-band speech takes into account the non-linear effect of coloured noise on the spectrum of speech i.e. some frequencies are affected more adversely than others.
Multi-band spectral subtraction (MBSS) method provides a definite improvement over other methods such as the conventional power spectral subtraction method.

Speech spectrum is divided into $N$ non-overlapping bands and spectral subtraction is applied to each band independently. The estimation of the clean speech spectrum in $i$'th band can be obtained from:

$$|x(\omega)|^2=|y(\omega)|^2 – \alpha_i \delta_i |d(\omega)|^2$$

for

$c_i<\omega<c_{i+1}$

Where $ c_i $ and $ c_{i+1}$ are the beginning and ending frequency band (your specific band) respectively. In the $i$'th band, $\alpha_i$ is the over-subtraction factor and $\delta_i$  is factor to control noise removal properties in each band individually. So you can set coefficient to remove the bounded noise from an audio signal.

For more help, you can read this paper :

  Sunil D. Kamath, Philipos C. Loizou, “A multi-band spectral
  subtraction method for enhancing speech corrupted by colored noise,”
   Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., vol. 4, July 2002.


Will the phase of the reconstructed data using a single MRI k-space point show any curvature? Some courses tell that each point in k-space (spatial-frequency domain) represents a certain stripe pattern in spatial domain.

Reference: https://www.sciencedirect.com/science/article/pii/B9780444534859000027
Consequently, if I reconstructed the image using only a single point in k-space data, the phase of the constructed image will change periodically(Exactly speaking, it changes linearly. The periodicity comes from phase wrapping). This means there will be only straight lines in phase map.
Here is a simple MATLAB code snippet for demonstration:
I = imread('cameraman.tif');
ksp = fftshift(fft2(I));

for i=1:256
   k = zeros(size(ksp));
   k(128,i) = ksp(128,i); % choose a single point in ksp

   Ir = ifft2(ifftshift(k));
   phase_map = angle(Ir);
   imshow(phase_map,[])
   title(num2str(i))
end

But some other courses say: the phase map will show curvature if the encoding gradient has a poor linearity. For example, the content of stanford RAD229 https://www.youtube.com/watch?v=uO--lYGjvFM begin from 14:40.
The k-space data itself won't know whether the gradient is good or not, but the meaning of each point wouldn't change. If each point of k-space only represents a certain stripe pattern. The phase of the reconstructed "image" using a single ksp point will change periodically (which means no curvature in phase map) even if the gradient has a poor linearity.
Conflict: Which one is right? Is there any way to generate the phase map showed in RAD 229?
Update:
Here is a small demo in 1D case：
import numpy as np
import matplotlib.pyplot as plt

mx = np.zeros([256,1])
mx[96:128+32] = 1
plt.plot(mx)

ksp = np.fft.fft(mx,axis=0)
plt.plot(np.abs(ksp))

k = np.zeros(ksp.shape);

# voxel #2 and #3 were mis-placed using an imperfect gradient
k[3] = ksp[3]+ksp[2];

rec = np.fft.ifft(k,axis=0);
phase = np.angle(rec)
# it's a straight line after unwrapping
plt.plot(np.angle(rec))

My question is about the variable phase in this code. If the video is right, the plotted phase should be curved. But the result shows that the phase is still  a straight line( after phase unwrapping).

I've acquired some images using a poor gradient unit. The images are distorted obviously but I cannot find any curvature when analyzing these data using the above mentioned matlab code.

Some courses tell that each point in k-space (spatial-frequency domain) represents a certain stripe pattern in spatial domain.

It does. A simpler way to think about what sort of "stripe pattern" would be generated by a point $(m,n)$ in k-space is to imagine a wave that propagates along the direction of $(0,0) (m,n)$. That wave has a frequency that is proportional to the distance of $(m,n)$ from the centre $(0,0)$. Also, you might find this and this useful too.

Consequently, if I reconstructed the image using only a single point in k-space data, the phase of the constructed image will change periodically. This means there will be only straight lines in phase map.

This is not strictly accurate, or maybe it is described in a vague way here. The "phase of the constructed image will" not "change periodically". It will vary linearly.
If the phase varied periodically, the reconstructed image would look like frequency modulation.
Phase is the first derivative of frequency.
Here is an illustrative example in Python (but easily generalisable in other platforms):
import numpy

Fs=100 # Sampling frequency in Hz
T=1 # Duration of signal in seconds

t = numpy.linspace(0,T, Fs*T) # Time vector (in seconds)

p = 2*numpy.pi*t # Phase vector

Here, p is our phase vector. As you can see, it is just a line that increases at a specific rate. If you feed this phase to a sinusoid, you get a specific frequency.
If you want to generate a higher frequency you make the slope that the phase varies steeper. For example y = numpy.sin(10*p). This will give you a 10 Hz sinusoid running for 1 second at a sampling frequency of 100Hz.
By saying that "...the phase [...] will change periodically..." what you imply is that p somehow oscillates. But what will happen if we get an oscillating p? (Frequency Modulation).
Here is a very simple example (it comes as a continuation of the above one):
from matplotlib import pyplot as plt

plt.plot(numpy.sin(20*p+6*numpy.sin(p)));
plt.xlabel("Discrete time")
plt.ylabel("Amplitude")
plt.show()

Which will give you something like:

You can see here that we start at a high frequency, we go to a lower frequency and then back up again to a high frequency. The phase here still "evolves" at a linear rate but on that ramp, we have superimposed a sinusoid to show what happens when "...phase change[s] [...] periodically...".
What you are asking, is actually answered in the Stanford video you cite within the first 5-10 minutes.
To understand the sort of distortion that the gradient brings to the image, you have to understand how does frequency map to space in the MRI.
The key to understanding that is the Larmor frequency. The Larmor frequency is what allows us to map frequency to space (or more accurately space to frequency) by varying the strength of the magnetic field.
In other words, by having a linear magnetic field, we map equal distances in physical space to equal distances in the k-space.
If the magnetic field varies non-linearly (first 5-8 minutes of the Stanford video), then where you excite (or receive from) is not where you think it is.
Therefore, equal distances in physical space do not map anymore to equal distances in the k-space.
This is what distorts the image.
An extreme example in 1 dimension is (again, as a continuation of the above)
plt.plot(numpy.sin(p-numpy.sin(0.25*p)));
plt.xlabel("Discrete time")
plt.ylabel("Amplitude")
plt.show()

This will give you something like:

This was supposed to be a "normal" sinusoid cycle but due to the fact that the phase does not change linearly, what we see here is a distorted version.
Here is the difference between the linear evolution of phase and the profile that causes this distortion:
plt.plot(p);
plt.plot(p-numpy.sin(0.25*p));
plt.legend(["Linear", "Non-Linear"])
plt.xlabel("Discrete phase")
plt.ylabel("Amplitude")
plt.show();

Which gives:

Again, the key to understanding what is going on here is to relate physical space to frequency and field strength through the Larmor frequency. If the gradient is not linear then equal parcels of physical space are not mapped to equal parcels of k-space and that is where distortion is coming from.
Hope this helps.
EDIT:
To the points raised:

...my question is whether the phase is a straight line or curved?

You seem to be referring to the phase spectrum as obtained from the Discrete Fourier Transform (DFT).
This has nothing to do with the reason behind the image distortion and impossible to "probe" from a single image.
The distortion of phase because of gradient non-linearities is something that happens BEFORE the data is captured. By the time the signal is captured it is distorted already.
Therefore, if you do not have a reference image (that is, an image that is distortion free), it is impossible to "see" the phase distortion. You could try infering it too, but that would be valid up to a point and for assumed distortion profiles. The only reason why we are able to tell that an image is distorted is because we know the anatomy.

I just want to reproduce the results in this video.

The results presented in the video come from a simulator. The very first equation presented at the top left is the one that is used to map gradient phase to k-space position. If you don't have $\vec{G}$, it is impossible to know how $k$ gets deformed for a given image.
You can remap the phase yourself with a function (that is, change the balance between the imaginary and real components of an image's complex DFT) and proceed with the reconstruction, or, you can achieve the end result of the distortion via an image geometric transformation. This transformation is literally a remapping of the pixel positions to new positions which is exactly what happens when the gradient contains non-linearities. For more information you might want to see this.


Interleaving over OFDM symbols in same transmission In this question I have asked what are the benefits from interleaving between two codewords over ONE OFDM symbol in packet based transmission made up of more than one OFDM symbol

Purpose of interleaving in OFDM systems

The answer was satisfying. I have another question now and I am looking for your opinion on this.

If I decide to do interleaving outside of a single OFDM symbol, (say for example in packet based transmission) then is there any benefit?

Thanks

Interleaving data that can be contained in one OFDM symbol between $n$ OFDM symbols will significantly reduce your bitrate and spectral efficiency by a factor of $\frac{1}{n}$. The advantage of doing this would be to increase your resilience to short bursts of noise, but only over $n$ symbols, at the cost of reducing your bitrate by $\frac{1}{n}$. I don't think many applications (if any) would benefit from this.

Edit: I think my initial response wasn't right. Suppose you're interleaving $n=10$ symbols with each other over 10 symbols, and only after you receive the whole block of 10 symbols can you read the code. The amount of information transmitted over the channel over time is the same as when symbols aren't interleaved, there is just some additional latency for the receiver to get the message. As long as your data can wait the extra time, it's worth it for your data to be that much more resilient to noise.


Would anyone be able to help out with the following discrete convolution question? I have been trying to attempt this question for quite a while and I am not sure how to approach it. Would anyone be able to help out with the working?

Hint

The simplest way is to use the Z transform property "convolution in time domain is multiplication in z domain".
See Z transform convolution

$$\mathrm{Z}(x[n]*h[n]) = \mathrm{Z}(x[n]) \times \mathrm{Z}(h[n])$$

Then you just need to do inverse Z transform.

Some typical Z transforms (including what is useful for you) can be found at Z transform pairs

Be careful with region of convergence (ROC).

It seems a homework?


What is the normalization factor for the quincunx lattice wavelet transform, and how do you find it? On pages 57-60 (preview was available last I checked, images here in case), there is a quincunx lattice transform described.

Lattice:

o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o

Basically you do these Predict operations on the black points:

x[ m][n ] -= 1/4 * ( LEFT + RIGHT + DOWN + UP )

Where $LEFT = x[ m ][ n-1 ]$, $RIGHT= x[m][n+1]$, $DOWN=x[m+1][n]$, $UP=x[m-1][n]$.

You then do updates on the white points:

x[ m][n] += 1/8 * ( LEFT + RIGHT + DOWN + UP )

Then you will never touch the black values again, so you effectively have:

o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o

You turn your head 45 degrees to see this is just another rectangular lattice, and you label them odd/even again:

o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •

You repeat this again and again, until you have 1 "average" left.

Now in the Haar wavelet transform, there is a power loss in each level that we correct with a normalization factor of √2.

Here, there is a computed power loss factor of about 1.4629 after the first step of the first level (found by running 5,000,000 transforms on random data and finding the ratio of powerBefore/powerAfter and averaging).

I don't know how to show/compute how this power loss is found, and where the 1.46 number comes from.

I don't think there is any single best number for normalisation because it depends on the structure of the values in your lattice.

In the simplest case where all values are equal, the predict operation zeros the black points and the update does not change the white points.  Because each  predict-update pair halves the number of nonzero points, multiplying the lattice by sqrt(2) after each pair of steps would conserve energy.

With all values independent with zero mean and equal variance, the predict step multiplies the variance of the black points by 5/4 and then the update step multiplies the variance of the white points by 281/256 so energy increases at each step.


Given MFCC values for some wav’s, how can I easily tell if all speakers are the same? Context: for a research project I have at regular times a batch of wav files  with recordings of voice. The hypothesis is that all speakers in a batch are the same, but there might be 1 or a few “others” in it. There are 5 to 10 files in a batch.

My goal is to obtain a similarity matrix for all wav files (speakers) in a batch.

What I got so far: I calculate MFCC (16 values) with an existing tool, I have between 2000 and 4000 frames per wav file. Simply taking the average of MFCCx is not meaningful. Next I perform a GMM clustering on each wav and obtain something like this: (MFCC0 discarded)

Wav1
Class   Mean (MFCC1)    Mean (MFCC2) …     Proportion
        1       -44             6               0,09
        2        46           -23               0,37
        3        20           -26               0,10
        4        43           -43               0,06
        5        57           -52               0,06
        6        66           -55               0,18
        7       108           -46               0,15

And here I am stuck, what do I do next to get a similarity matrix, can I calculate some “distance” between the different wavs and how?

Voice is natural talk. Recordings are under various channel conditions (sometimes background noise) and can be with varying emotional status (eg laughter, some excitement,..). Speakers are of random nature (man, woman, old, young, …).

I would like to avoid more difficult algorithms for speaker recognition because that is beyond my current knowledge. Given that speakers are random, it will be rare that the other(s) is very similar to the normal speaker, and some error rate is acceptable, so a basic method should be ok.

Given MFCC values for some wav’s, how can I easily tell if all speakers are the same?

In general, you would deal with this kind of problem using some form of Clustering.

The fundamental premise behind this idea is that the features of similar data items would have similar values.

How you assess similarity, how you assess "a cluster" and what this means for the particular application you are dealing with are very much context specific issues.

In this particular example, the features are the MFCC coefficients. If these are similar between two frames we would tend to think that the frames contain more or less signals produced from the same source.

The simplest clustering to apply would be k-means clustering but this is not going to help too much here, especially if you don't know how many speakers you might get in a recording. So, long story short, you will most likely gravitate to something like hierarchical clustering.

You don't have to extract the means of each feature to cluster them. Doing so would be equivalent to training and using a minimum distance classifier. The only thing that is needed here to complete that is the "ground truth" knowledge that will tell you which centroids belong to which speakers.

Note here that if you have two speakers talking in a mostly silent room, the "feature space" of the clustering is not necessarily expected to be divided into two distinct "point clouds". It could, for example, contain separate regions for background silence, one speaker talking, the other speaker talking, both speakers talking at the same time. This last class might contain further grading depending on how loudly was each speaker talking while interrupting the other.

Hierarchical clustering allows you to focus on details with varying degrees of "attention", provided that the MFCC clusters the feature space with sufficient accuracy.

The other thing you can do, is a Principal Component Analysis on the MFCC matrix which (hopefully) will return one component for each one of these "major events" happening in the file.

Hope this helps.

EDIT:
Forgot to add, if you are looking for the terminology that could return to you some useful results in varying levels of complexity, this is "Speaker Diarisation". See for example here.


Available Online Recordings for Acoustic Echo Identifiablity I am trying to add a few samples for a publication on identification of clean signal from the one which has echo on it. For this purpose I need data set that includes playback of an already available signal in a big room/hall and the recording of the signal with a mic in this hall that includes the echo.

Does anybody know of any available online source for such data samples.

I don't know of any online recordings, but you could try using the Silophone player to create some.  Play your original file, then record what the silophone plays back.
Silophone homepage

The Silophone uses an old grain silo as an echo chamber.  It also plays the echoed audio locally over a speaker outside the silo, so you might want to avoid using anything obnoxious or private or copyrighted.


Increasing the number of bins in DFT with limited points and high samplerate Apologies for the potentially simple question in advance. I have a limited understanding of the subject, but need to apply this to a project on short notice.

My goal is to run a DFT on a signal with N points and sample rate of R. The signal is often much shorter than 1 second, hence R is much larger than N. This causes the resulting DFT bins to skip many frequencies I'm interested in. As I understand, bin intervals are related to R/N, so I need to either increase N or decrease R. I was thinking about trying one of the following:

Stretch the signal by adding intermediate points which are averaged between real points so that N is large enough, then transposing the resulting frequencies back down by the stretch factor.
Appending to the signal with either 0s or another copy of itself to increase N.
Downsample the signal to lower R.

Are any of these ideas reasonable? I can not gather more points of the signal.

Edit: forgot to mention that the frequencies I'm interested in are very low (~20-300hz), so lowering the sample rate (usually 44.1k) shouldn't be an issue.

If all you want is a finer bin spacing, then the most straightforward thing to do is to increase $N$ by zero-padding. Note, however, that this will not increase your frequency resolution. Frequency resolution is limited by the total length of time that you observe the waveform. So, if you want to actually be able to resolve frequencies more finely, you need to observe the signal for a longer period of time.

Your first option (essentially interpolating the signal to a higher sample rate) won't do what you want at all; $N$ will be larger, but so will the sample rate of the interpolated signal. If you double the size of the dataset by interpolating by 2, then you'll have to do a DFT of twice the length in order to get the same bin spacing.

Downsampling the signal to a lower sample rate as you suggested in your third option can make sense if, as you said, you only care about frequencies below 300 Hz and your initial sample rate is 44.1 kHz. However, what I said above applies: this won't improve your frequency resolution at all; you'll need to either live with zero-padding to interpolate the DFT spectrum, or observe the signal for a longer period of time instead.


Why is distorted / quantized white noise still white noise? White noise can be created by setting every sample in the time domain to a random number with a gaussian distribution. Every sample is random, so the frequency response is obviously random. Easy enough.

Next step: Set every sample to a random number between 0 and 1, but the numbers are equally distributed instead of using a gaussian distribution (please ignore the DC offsets, they don't really matter). Samples are random enough, so it still makes sense that it would sound like white noise.

Now, set every sample to either 0 or 1 with a 50/50 chance. This feels like it shouldn't sound like white noise anymore, but it still does.

Finally, set every sample to either 0 or 1 with a 90/10 chance. How does this still sound like white noise when it looks nothing like it?

EDIT
I did some experimenting (with sample rate 48000 Hz).
1/2 chance is basically white noise.
1/4 chance sounds almost like white noise, but if you know what you're looking for it's easy to tell the difference.
1/16 chance sounds like rain smattering against something (you can hear individual impulses), with white noise in the background.

Varying the amplitudes of the impulses did not seem to have any effect on the sound, other than volume.

White noise is defined as having a flat spectrum over the entire frequency range. Using uncorellated random numbers is just one of many methods to generate this, but by far not the only or the best one.
Also, whiteness implies no gaussian or whatever distribution, so your first sentence is wrong twice.


time linearity of LTI system There is this following question.
Consider the transformation H {x}[n] = n x [n]. Does H define an LTI system?
What I understood from the question is:
x[n] -> H -> y[n]
or according to the question,
x[n] -> H -> n x[n]
to prove the time invariance, if I put n-no instead of n, it should be,
x[n-no] -> H -> y[n-no] = n x[n-no]
but it turns out to be,
x[n-no] -> H -> n x[n - no] - no x[n-no]
which means the output depends on extra term no x[n-n0] and is not time invariant.
Incase, the trasformation was any constant * x[n] the it would have conserve the time invariance property, if I am not wrong.
But the proof mentioned in the notes is,

I have difficulty understanding the proof. If anyone could explain it to me or if I understood the question wrong, I'd appreciate it.

Your notes gave a counter example, you just proved it in the general case.
Your note used your x[n] as delta and showed that
y[n-1] unequal to H{x[n-1]} which is exactly what you showed.
In the notes example H{x}=0 for all n, as discrit delta equals zero foll all n but n=0, which means H is zero every where for all input of the kind delta[n-n0]  (becouae of the n factor which zeros y at its only non zero point)
Than you have zero output for delta input. But as you showed time varied output isn't zeros every where.


Inter aural time difference generation using sine waves I have a stereo sound for Inter  aural time difference consisting of sine wave(16 samples as 16 bit PCM Data) for left and right audio. From time to time I would like to advance either the left or right audio so that sound in one ear comes earlier than the other. If I just put another sine wave of different frequency in between my normal sine samples to give an illusion of changing sound in one ear, I hear a click. I am playing a buffer of 1024 samples at a time containing looped sine wave samples. If I drop one sample it generates a click as well. Can someone suggest a way of speeding up the Sine wave without generating a click.
e.g loop the sine wave separately for both ears. Now advance the sine wave in one ear so that it feels like a frequency change and then play the normal sound again in both ears. I have fixed sine samples with me. (I work in java so can't do any heavy computations). I have tried to clarify but the question could still be confusing. Ask for specific details.

Thanks

I understand the point.
If you are outputting signal[i] you cant jump to signal[i+n] instantaneously.
You should just fade it to the new point with some window function. Allow some period of fading time.

For example: If I currently output signal[i] and I want to shift it to signal[i+n] in 1 second. I will use half of Hann window to fade it in and out. So I time scaling Hann from pi to 1 second.
and output $signal[i+n]*Hann[i]+signal[i]*(1-Hann[i])$ for 1 second before I switch to signal[i+n] at full scale. This should not generate a click.

You may want to try other window function or different fading period.


How to identify stable deltas in a noisy signal This is a problem I have had to deal with in several different contexts now.  Given a real time signal which you expect to act as a noisy step function (periods of stability separated by quick changes in magnitude).  How does one identify when the signal has stabilized?  The tension in this problem is between identifying stability as fast as possible and reporting a new stable value with as much accuracy as possible.

My solutions have varied depending on which of those I felt was more important and on how good my transition model was.  I have taken to calling this the stable deltas problem but have been unable to find any references on how to approach it.  Is this a studied problem in signal processing?  Does it have a name?

This sounds like the "change detection" problem of Basseville and Nikiforov, with the parameter after the change unknown (see section 2.4) of the linked PDF.

See, for instance, example 2.4.3:

The usual algorithm for dealing with this, and the one analyzed in the above-linked book is the CUSUM (Cumulative Sum) algorithm.


What is meant by "omega suppression" in a lock in amplifier context? Can anyone explain what is meant by "Omega suppression" in lock-in amplifier discussions or specifications. I'm usually very good at web-searching, but have been unable to find anything that explains rather than mentions this term. I find it mentioned in some manuals but never explained.
Omega and 2 Omega sound like fundamental and 2nd harmonic, but ...?
Some mentions:
In passing only here Setting up a Digitally Tunable Resonator using the HF2LI
Many mentions here MFLI User Manual but meaning uncertain from context.

Looking at this page, the part where they explain "Signal Dynamics and Demodulation Bandwidth" (search for suppress), they give an example of AM demodulation by comparing simple filtering with a locked-in filter. In the text, they say:

For the second set of requirements, frequency components corresponding to the sidebands are rejected by reducing the filter bandwidth to a value smaller than fm [...]

I don't know whether this is what "omega suppression" is about, but it seems to fit the description. What's more, the link you give is for the same datasheet as theirs (newer, to be expected). So my guess is the device has an option for suppressing the sidebands, as described on page 59 on their datasheet (the newer one):

Multi-harmonic  mode  (optional,  simultaneous  measurement  up  to  four  harmonicfrequencies)
Arbitrary frequency mode (optional, simultaneous measurement at 4 arbitrary frequencies

I didn't read everything, just glanced, so based on these, this is my conclusion. Maybe there are others.


Path Center Correction So I have groups of (x,y) coordinates that I'm using to draw lines.  The lines aren't perfectly straight but they are connected.

In a Cartesian plane the middle would be called the origin, but I'm working with a digital image processing program, so the coordinate (0,0) is at the top left of the image.

What I'm trying to do is translate all of the coordinates in each line (group of coordinates) the same distance, so that the first point in each line starts at the middle of the image.

Here is an example of what I want:

Here is what I have:

1) Group of points: $\{ (x_{i},y_{i})\}_{i=1..n}$

2) Centroid of the group of points: $\left(\frac{1}{n} \sum^{n}_{i=1}x_{i},\frac{1}{n}\sum^{n}_{i=1}y_{i} \right)$.

3) Center of the image: $(x_{c},y_{c})$

4) Translation of arbitrary point $(x,y)$ to the center $(x_{c},y_{c})$: $t(x,y)=(x_{c}-x,y_{c}-y)$

Note that $(x,y)+t(x,y) = (x_{c},y_{c})$ ...

Now just substitute 2) in 4) and you have your formula to translate the centroid to the image center. Right?

I am sure that I am missing something here - the translation distance obviously cannot be kept the same for all points. Didn't you mean another type of transform (maybe scaling have to be included as well) ?