Measuring THD using consumer sound card I am trying to measure (estimate) THD of an audio path using a consumer sound card. The idea is that I do two measurements, one without the system being measured in the path and one with it:

MEASUREMENT 1:

Sound card DAC ---> Sound card ADC

MEASUREMENT 2:

Sound card DAC ---> Measured system ---> Sound card ADC

There is the same sine wave generated in both measurements. In each measurement, the ADC data is transformed using FFT (with Hann windowing function) and the THD is calculated by the standard THD(F/R) equation (https://en.wikipedia.org/wiki/Total_harmonic_distortion) from the amplitudes from the peaks of harmonics on the FFT.

The idea is that by doing two measurements, one with and one without the measured system, is that I could eliminate the TDH introduced by the non-ideal DAC, ADC and sine wave generated and (only) measure TDH of the measured system.

Since the THD is not a linear phenomenon, I assume this is not a correct way to estimate THD of the measured system:

THD_measured_system = THD_measurement_2 - THD_measurement_1

So my question is what would be the correct approach in eliminating the THD of the measurement system and if this is even a possible approach?

EDIT:

Would checking the difference between amplitudes of the harmonics between measurements 1 and 2 and calculating the THD from this difference work?

The main prerequisite for measuring anything meaningful is that the target parameter (THD) is much better in your measurement system than in the system under test. The higher the difference, the more accurate your result will be. If your measurement system is orders of magnitude better, that the system under test, the measurement error can be neglected. If it is not, your result will have a considerable unsecurity. "Eliminating" the error introduced by your measurement system is not possible, because, like you said, it's inherently non-linear.

For example: if the THD of your sound card, is just 10dB better than the THD of the system under test, the uncertanity will be at approximately 30%. With 20dB it will be 10% and with 30dB it will be 3% and so on.

Regarding your edit: No, this would not work well. As we said, THD is non-linear, not only in the sense, that new frequencies are generated, but also regarding the phase of those new components. E.g., if both systems have the same THD individually, it would not just add up (i.e. be 6dB higher), but could be anywhere from +6dB down to -inf. All information, that can be garnered from comparing the overall THD with the measurement THD comes with uncertainity regarding the individual values. That is a direct consequence of the non-linearity.


Normalizing white noise to match pink noise's amplitude at a given frequency? I have a pinking filter which I am running white noise through to get pink noise. I would like to then mix the white noise with this pink noise, where the amplitude of the white noise will match the pink noise at a given frequency.

My pinking filter is +/-0.5 dB above 10 Hz.

So let's say I want the white noise to match the amplitude of the pink noise at 100 Hz. Or at least be standardized to close to this frequency's amplitude.

I believe in terms of power which is amplitude squared for pink noise we have: power = 1/freq.

So would it be: amplitude = sqrt(1/freq)?

And thus amplitude of pink noise at 100 hz would be 0.316? I know that's not correct but I don't know what is.

Basically I'm looking for a multiplier for the white noise so that at that given frequency it will match the amplitude of the pink noise.

Thanks.

Okay, when discussing white noise or pink noise (or red noise or brown noise or flicker) or some other random process, there is this property called the power spectrum, in which white noise has a constant value for all frequencies.  But we integrate the power spectrum over all frequencies (from $-\infty$ to $+\infty$) to get power.  Integrating a constant, non-zero function over an infinite frequency range gives you infinite power.

So in reality, white noise does not exist.  (It puts out more power than the Sun.)

But if you bandlimit white noise to a finite bandwidth, you have a finite power.

Sampled signals (in discrete time) must have a finite bandwidth in order to be sampled properly.

If you were to synthesize bandlimited white noise by use of a random number generator that produces good random numbers with a uniform distribution between the values of $0$ and $\Delta$, you would have a mean of $\frac{\Delta}{2}$ and a variance of $\frac{\Delta^2}{12}$.  If you subtract from each uniform random number the mean of $\frac{\Delta}{2}$, then you would have the same variance, but a mean of zero.

Now if you were to add $N$ of those independent zero-mean uniform random numbers together, the variance would increase to $\frac{N \Delta^2}{12} \triangleq \sigma^2 $ and the distribution of the random numbers would closely approach that of a normal distribution (sometimes called Gaussian).

Now variance of a zero mean (or DC free) random process is the power, which is the area (or integral) of the power spectrum.  So from $-\frac{f_\mathrm{s}}{2}$ to $+\frac{f_\mathrm{s}}{2}$, you would have a flat power spectrum with an amplitude of $\frac{\sigma^2}{f_\mathrm{s}}$ so that the area or integral comes out to be $\sigma^2$.  (The amplitude in the power spectrum is zero for all $|f|>\frac{f_\mathrm{s}}{2}$.  This is what makes it bandlimited.)

So now you have "white" noise that is finite energy.  When you run that through a pinking filter (is that what you are doing to make pink noise?), then you have to multiply that constant amplitude of $\frac{\sigma^2}{f_\mathrm{s}}$ times the magnitude-squared of the pinking filter to get the frequency-dependent amplitude of the pink noise.  Then you have to integrate that from from $-\frac{f_\mathrm{s}}{2}$ to $+\frac{f_\mathrm{s}}{2}$ to get the power of the pink noise.  Whatever is the ratio of that power to the power of the bandlimited white noise going in, you compute the square-root of the reciprocal, and that is the scaling constant to multiply times your pink noise to make it have the same power as the bandlimited white noise going in.


Difference between "scaled white noise " vs. "unscaled white noise" I have been searching Matlab user-guide for Wavelets. In their graphical user interface for denoising, there is an option for "unscaled white noise" and "scaled white noise" for the noise structure. What is the main difference between the two? Does the scaled white noise imply that the noise is amplitude changing with time? I have searched Google books and Google Scholar, nobody defines it. Thanks.

Wavelets, as bases or frames, are linear like-decompositions. As such, signals and processes can be decomposed into a linear or weighed sum of coefficients multipling (wavelet) vectors:

$$s[n] = \sum a_n e[n]$$

The indexing above is not wavelet-specific. Traditionally, there was some natural order: low $n$ denoting low-scales, or high frequencies, and vice-versa. So traditional linear filtering often consisted in keeping indices above or below some cut-off indices (with weighting), assuming that data above and below could be noise (typically low-pass filtering).

This linear (order-related) scheme remains practically useful. However, and some cases, a non-linear selection can prove more efficient. Especially when indices, across time and scale, don't exhibit clear ordering anymore.

$$s[n] = \sum a_{j,k} \psi\left[\frac{n -k2^j}{2^j}\right]$$

The main idea is to select "cleaned signal" coefficient not within indices, but with respect to magnitude: the highest the magnitude, the most likely it belongs to data instead of noise.

Selecting the highest coefficients amounts to setting a threshold below which you discard coefficients. This threshold should be selected as cleverly of possible. nd here is the catch. In the early works on wavelet thresholding, for proofs and practice, it was often considering that the noise was Gaussian with (known) unit variance. This would yield nice theories and practical algorithms.

Yet, this requires a knowledge on the noise, which is often unknown (in shape and amplitude). If we forget the shape (Gaussian or not), if the noise level is $\sigma$ (instead on one), this should be taken into account, esp. by changing the threshold accordingly. Caveat, estimating this $\sigma$  remains tricky.


Can an audio recording tell whether it is made in a large or small room? The conversation of one or more people will be analyzed. I think need to look for something similar to reverberation or echo in an audio signal. How can this be implemented?

Update: small room is about 10-15m2, large room 50-100m2

There could be a number of approaches for this task, depending on the number of channels (microphones) in your input and structure of available data.

Given a labeled data set from a previously known set of rooms, you can train a neural network (NN) to classify in which one of them a new sample was recorded.

A second approach would be to estimate the parameters of the room impulse response (RIR). Certain methods could be under the limitation of a single microphone. A microphone array would supply more information like in this paper. Depending on the shape of the array, different methods can be utilized.

Also, do you have any other prior knowledge on the room? For example, can you get a clean (without speech) audio segment? Do you know where is the microphone and sound sources are located on the room? These priors can point you to methods that utilize this information and improve the estimation.


Minimum number of Poles and zero of transfer function H(z)? Suppose $G(z)=H(z)(1-\frac{1}{2}z^{-1})$ now in question its saying ROC of G(Z) is entire Z plane except Z=0,so here we need not to add anything because G(Z) already a right sided signal with ROC entire z plane except at z=0.

now  $H(z)=\frac{G(z)}{(1-\frac{1}{2}z^{-1})}$, for the input $x(n)=(1+j)^n$,
  H(z) should become zero.

Now here i am not getting what should be the value of G(z) so that its output becomes zero for the input $(1+j)^n$.

Restricting the discussion to LTI systems with rational transfer functions, there are various conditions on the placement of poles and zeros  of $H(z)$ implied by several properties of the system under concern. Using some of them helps in this question.

Prop 1: Total number of poles equals total number of zeros (including those at $z=0$ and $|z|=\infty$) ,i.e., There are equal number of poles and zeros for any $H(z)$.

Rule 2: For real systems (which has real impulse response $h[n]$) poles (or zeros) exist as either purely real or as complex conjugate pairs.

Now applying these two properties to your system which is real and causal (how?) we reach the following conclusion:

1- If ROC of $(1-\frac{1}{2} z^{-1}) H(z)$ is the entire $\mathcal{Z}$-plane (exc $z=0$), then $H(z)$ has a pole at $z=0.5$

2- If $x[n]=(1+j)^n$ produces zero output, then $z=1 \pm j$ are two zeros of $H(z)$

Are there any more poles of $H(z)$ other than those at $z=0$ ? No. Because if there were so, then part one would indicate it in the ROC of $(1-0.5 z^{-1})H(z)$ which is given as entire $\mathcal{Z}$ plane indicating that there are no more poles of finite, nonzeros values.

Are there any more zeros of $H(z)$ ? The given information is not sufficient to conclude that. There could be more zeros. If they were real then they could be even or odd numbered, if they were complex then they must be in complex conjugate pairs. We only know one complex pair of zeros at $z=1 \pm j$. Note that for each added zero to $H(z)$ there would be one more added pole at $z=0$ (in addition to the existing pole at $z=0$)

Hence we conclude: $H(z)$ has two poles at $z=0$ and $z=0.5$ and at least two zeros at $z=1+j$ and $z=1-j$. If $H(z)$ had any more zeros at $z_k=\alpha_k$ then there would be equal number of poles at $z=0$. And if the zeros were complex they would be in complex-conjugate pairs (and the added poles would be doubled)


transient response of simple digitalized RC low-pass filter In http://en.wikipedia.org/wiki/Bilinear_transform#Example, digital version of simple RC low-pass filter is presented:

$$\frac{1 + z^{-1}}{(1 + 2RC / T) + (1 - 2RC / T) z^{-1}}$$

where $T$ is sampling interval. Let me reorganize by $A = 2RC/T$, and get the following:

$$y(n) = \frac{x(n) + x(n-1)+(A-1)y(n-1)}{A+1}$$

But I cannot see how I would get transient and steady response back from this difference equation. Can anyone help here?

If you assume a zero initial condition (i.e., if the system is switched on at $n=0$, then $y[-1]=0$), then the system is fully described by its impulse response, i.e. its response to the unit impulse $x[n]=\delta[n]$ (which is zero everywhere, except at $n=0$, where it is $1$).

Write the system equation as

$$y[n]=b_0x[n]+b_1x[n-1]-a_1y[n-1]\tag{1}$$

In your example $b_0=b_1$. If you apply $x[n]=\delta[n]$, you get the following response:

$$\begin{align}y[0]&=b_0\\y[1]&=b_1-a_1b_0\\
y[2]&=-a_1(b_1-a_1b_0)\\
y[3]&=a_1^2(b_1-a_1b_0)\\\vdots\end{align}$$

So we get for the impulse response

$$h[n]=\begin{cases}0,&n<0\\b_0,&n=0\\
(-a_1)^{n-1}(b_1-a_1b_0),&n>0\end{cases}\tag{2}$$

The general input-output relation is described by the convolution of the input signal and the impulse response:

$$y[n]=\sum_{k=0}^{\infty}h[k]x[n-k]\tag{3}$$

E.g., the step response can be easily computed from the impulse response:

$$a[n]=\sum_{k=0}^{n}h[k]\tag{4}$$


How do I convert this simple Laplace equation to Z-domain? A basic model of coupled strings (eg. piano) is provided here as:

The principle is that it has two identical string simulations each formed by a delay line and LPF. The outputs of these are summed at the "bridge" in the center. The sum is multiplied by $-H_b(s)$ to get the bridge output, which is then subtracted from each of the string loops to simulate damping from energy lost to the bridge.

The bridge equations are:

$$ V_b(s) = H_b(s) \big( V_1^+(s) + V_2^+(s) \big) $$

$$ H_b(s) \triangleq \frac{2}{2+ R_b(s)/R} $$

Where $R_b(s)$ is the impedance of the bridge and $R$ is the impedance of the string.

So in order to implement this, I need to convert $H_b(s)$ (and/or $V_b(s)$) from a Laplace equation to a Z-domain equation.

I am starting to understand some of the principles of the Laplace and Z-domain but I have never converted anything before and I am not sure how to do this. How I get $H_b(z)$ for this usage?

Once I have $H_b(z)$, do I just multiply it by the sums of the two string outputs and then subtract that result from the string loops?

Also, why is the impedance of the string simply $R$ while the impedance of the bridge is $R_b(s)$? Is impedance of the bridge a Laplace function to imply that its impedance may vary with frequency? (Whereas the string impedance would be frequency independent.) If so, how might it simplify things to treat the bridge impedance as frequency independent too?

Thanks.

I believe, mike, that the answer to your question is likely the Bilinear Transform.  Make sure you identify the significant frequency (likely the resonant frequency of the LPF1 and LPF2 or $H_b(s)$ and apply prewarping to that frequency, so that the digital filter hits at the same place that the analog filter.  Remember the transfer functions of the two delay lines are $z^{-N_1}$ and $z^{-N_2}$.  Or, in the s-domain $e^{-N_1sT}$ and $e^{-N_2sT}$ , where $T$ is the sampling period or the reciprocal of the sample rate.


Analysis and algorithmic complexity of the Constant Q transform Here I am again (Implementation of the constant Q transform + property questions) :)

I was reading this paper on a constant-q transform algorithm, where they managed to turn the following equation for the constant-q transform:

into this:

(where X[k] is the fourier transform of x[n], and K*[k, kcq] is the complex conjugate of the fourier transform of the w[n, kcq])
- claiming it to be an efficient algorithm.

While I'm still new to algorithmic complexity, I tried to derive the big-O complexity. I'm assuming the transforms of the windows (ie. K*[], the kernels) are pre-calculated. Let N be the size of the transform, and K the number of components (kernels / filters / results). Would this be correct?

So thats the FFT (N * log2(N)) plus K loops of N multiplies (the convolution) + additional stuff done at each loop (done K times).
Also, since the kernels are precalculated, that incurs both a certain amount of memory storage needed + a lot of memory access during the calculation. With their implementation, they store the kernels in full size (that is, N), so the required storage memory for the calculation would be K kernels of N size.

I'm defining M() to be the amount of memory used during an calculation:

That is, we need to load the whole signal (N), also K number of N-sized complex (2) kernels. That is quite different from algorithms like the standard fft and direct evaluation of the constant-q transform which requires no working storage, and as of such are defined to be:

-Since they only need to load the signal. If we define the complexity of the FFT and the direct evaluation of the constant-q transform to the following:

(This assumes direct evaluation of the CQ-transform which incurs at least NK operations, together with NK operations calculating the discrete window for each component.)

we can create a small comparison table for an example. If the transform is of size 2048, and the wanted number of components is 1048 (N = 2048, K = 1048), we get the following results:

It would seem the cqfft algorithm (theirs) incurs a rather extreme memory cost (that in practical appliances could dominate computational speed) for a mere ~2x increase in computation speed. Is this correctly derived?

Here's a reference-pseudocode of their implementation i wrote, as i understand it. If i have misunderstood their algorithm, it will probably show here:

std::vector<complex<float>> fft(float[], size_t);

class kTransform
{
    vector<complex<float>> cqfft(vector<float> signal)
    {
        vector<complex<float>> fftdata = fft(signal, size);
        for (int kernel = 0; kernel < kernels.size(); ++kernel)
        {
            complex<float> accumulator;
            for (int k = 0; k < size; ++k)
            {
                accumulator += fftdata[k] * conj(kernels[kernel][k]);
            }
            fftdata[kernel] = accumulator / size;
        }
        return fftdata;
    }

    void compileKernels(vector<float> freqs, size_t dftSize)
    {
        auto numKernels = freqs.size();
        kernels.resize(numKernels);
        size = dftSize;
        for (int kernel = 0; kernel < numKernels; ++kernel)
        {
            kernels[kernel] = fft(Window(freqs[kernel]));
        }
    }
    size_t size;
    vector<vector<complex<float>>> kernels;
};

Thank you for any insight. As a bonus question, does anyone know what the complexity of a fast implementation of the Gabor transform is?

I didn't verify each exact step of your analysis, but it seems that you're overestimating the memory part by a huge margin. A lot of the kernels are identical except for the exact tiem moemnt at which they're applied. In fact, you only have one kernel per level. So for a typical music application you have only 12 levels per octave, and thus 12 kernels. Over a 10 octave range (20 Hz-20 kHz) that's still only 240 levels. The paper you reference uses quarter-tone resolution (24/octave) and 5 octaves. That's just a design choice.

Next, you state that "With their implementation, they store the kernels in full size". I doubt it. The paper explicitly mentions the effect of trimming the kernels (see Figure 3) and states that the kernels are trimmed to 6 values max (and sometimes just 1 !). This saves both memory and time. The exact amount does depend on the trimming applied.


Can you give example applications using single bin sliding DFT? I need to detect a single tone signal in real-time. The amplitude and the phase of the single tone signal are changing in time. I want to detect both the amplitude and the phase of the signal. I do not want to do FFT since I only need to see a single frequency bin and also FFT is not suited well for real-time applications. I see that there are sliding Goertzel and sliding single bin DFT methods. These methods seem easy to apply and both need very fewer computations compared to DFT. However, they seem like have stability issues. Some damping factors are used to solve the stability problems. However, I still have doubts on the sliding DFT method. Is there any point that I need to pay attention while applying single bin sliding DFT? Can you give examples of real world applications that are using single bin sliding DFT?

I do not want to do FFT since I only need to see a single frequency bin and also FFT is not suited well for real-time applications. I see that there are sliding Goertzel and sliding single bin DFT methods

FFT is just an implementation of the DFT. So, for the same reasons you don't want the FFT you wouldn't want any other DFT.

Goertzel is basically a way of getting the same result any DFT would have given you, but only for a single bin, exactly. Goertzel's latency, when implemented in parallel hardware, is the same as the overall DFT's, just with less ressource usage. In software, you'd use a bunch of operations.

To detect a single tone, my recommendation here would simply be a correlator (which happens to be the same as a matched filter) with the tone of interest. Basically, take the tone you're interested in, and use it as filter taps (apply windowing to taste). You can easily shift in single samples into such a filter. If you don't window (==rectangular window) and convolve a complex sinusoid with a whole number of samples as period, you get what you'd get as a single bin DFT: it's the same formula.

The whole truth, though, is that you don't really get a useful estimate for the presence of a periodic signal with every new sample – no matter how you twist and turn this, it's always a rate-reducing operation ($n$ samples go in, 1 goes out).

So, you might want to consider your definition of "real-time detection". How could one sample change the fact that there's a let's say 200-samples-long periodicity in the signal?

Depending on what you implement this on (an FPGA? ASIC? DSP? MCU? PC-style computer?), you'd want to pick a different implementation, and you'd be getting samples in buffers, anyway – so, I'll go out on a limb here and claim that if you are writing software, you can probably get away with a simple narrowband FIR (and you'd tune that bandwidth and transition width to make your detector robust to slight deviations in frequency), and it's often most efficient to implement that using fast convolution based on the FFT.


Direction of sound using microphone array I want to estimate the direction of arrival of sound (angle in a 2d plane) using a microphone array. I'm hoping to follow a time-delay based method using a 4 mic array. The method requires me to find the time delay of the signal upon arriving at the different microphones. For example, in the 4 mic array, I'm required to find the time delay of the signal between mics 1 and 2, 1 and 3, and 1 and 4 and so on for all other pairs. It would be very convenient if I could have a 4 mic array that could simply be plugged into my laptop and I guess an NI Labview algorithm would take care of the rest. And some one somewhere suggested that the microphone array can be created by using 4 usb (mini) mic and connecting them to a USB hub. First of all, can this be done? Secondly, since its a USB mic, the signal will already be digitized right? Any suggestions regarding how to extract information about time from this? I do not have any prior experience with signal/sound processing.

If you connect four discrete USB microphones a challenge will be to synchronize them. Even when you start them simultaneously (from a programming point of view) you cannot be sure that they all start the audio acquisition at exactly the same time so there will be an inherent unknown time delay between the microphones. This can be solved by having a reference signal coming from a known direction. Then, by comparing the expected delays with the actual delays, synchronization can be achieved. With your setup I would expect the synchronization delays to be slightly different every time the system is started so that calibration needs to be done every time.


Noise addition for aesthetic 'film-grain' modification of RGB channels in GLSL forgive me for the lack of knowledge on this, but I'm hoping someone can point me in the right direction.

I'd like to know the best technique for adding noise to RGB channels individually so that magnitude can be varied individually.

What algorithm of noise is the most aesthetically pleasing, or organic looking?

Many thanks.

Film grain is visible in images and is due to the "analog pixels" that make up the film emulsion. These are crystals of silver halides at various sizes, the volume of which determines the film's "speed".

The larger the crystals covering the film's surface the higher the chance they will capture photons, the lower the spatial resolution they allow to be recorded and vice versa. So, a 400 ASA film is more likely to be used in low light conditions (e.g. natural light in night scenes) than a 50ASA film. However, during development, the 400ASA film will have visibly more "grain" and consequently its spatial resolution will be low.

It is impossible to accurately control the size and distribution of crystals on the film emulsion. That is, "accurately" with respect to the fine control we have with CCD sensors. But they are very small and at huge numbers over a film surface and this is why film resolution is unparalleled.

And that is where the "aesthetic pleasure" of the "noise" is derived from: The volume and spatial distribution of tiny little crystals.

To cut a long story short, if you want to add realistic grain to a digital image, shoot a grey card at varying exposures, develop the film (without adjusting for exposure at the developing machine), scan the negative and then, literally, add it over the top of an image, possibly with a bit of "multiply" to modulate your original image.

Now, obviously, the grey card will produce a grey image. An ideal grey image has identical values in the R,G and B channels. But depending on what you are simulating, that might not be the case. This type of film (B&W) has a single photographic emulsion and responds to total light flux, but this type of film has multiple coatings and therefore it is more likely to "add" some coloration. In other words this second one might have a different "character". In fact, different processes used by different manufacturers lead to films with different characteristics anyway. You will know an ektachrome picture when you see it for example. And this includes the characteristics of the noise.

OK, so the other question is "How do I simulate this?".

Well, if you have to, then the key idea is to take some noise distribution, apply a low pass filter over it and then take this as your "noise overlay". Varying the parameters of the distribution and low pass filter you will achieve different qualities of noise (but this is not "film grain").

So, take a uniform random distribution image, pass it through a low pass convolution matrix (or blurring filter) and you have some noise overlay.

A more aesthetically pleasing result might be achieved by using a pink noise distribution. You can generate this very easily via filtering your uniform noise image (from above) with a spatial filter that has a response of $\frac{1}{f}$, where $f$ is the spatial frequency. This is close to the physical process, from the point of view of the distribution of the volume of the crystals on the film. In other words, fast film, larger crystals, film more sensitive but lower spatial resolution, i.e. higher spatial frequencies (small objects) less likley to be recorded by the film (that is, the top end of $\frac{1}{f}$). But this is not film grain.

Other distributions may be used similarly with different results. If you want to simulate some additional coloration (as noted above), of course, you will need to multiply the R,G,B components of the grayscale noise overlay with different factors depending on what you are trying to achieve (e.g. to simulate tinting).

Hope this helps.


MSE in adaptative filter I'm trying to filter some motion noise from an ECG signal. To do that I'll try to implement an adaptive filter using the LMS algorithm.

I think I have to calculate the MSE of this:

E[e^2 ] = E[(s + n )^2 ]+ 2E[(s + n)X ]W^T + WE[X^T X ]W^T

in which s+n is the noised ECG,

X (n) = [ 1, Accx (n), Accy (n), Accz (n)] are the values from the accelerometer,

W = [w0, w1, w2, w3] the coefficients.

Wt +1 = Wt + 2αeX,

I'm trying to understand how to calculate alfa and minimise the MSE.

http://www.mathworks.com/help/dsp/ref/adaptfilt.lms.html

can anyone give me a simpler example? I'm kind of lost here.

This scilab script implements a simple LMS adaptive filter.

M = 50;
LMS = zeros(M,N);
LMS(:,1) = zeros(M,1);
ERR = zeros(1,N);
y = zeros(1,N);
mu = 0.0005;

for t=M:N,
    Uwindowed = u(t - [0:M-1]');
    y(t) = LMS(:,t)'*Uwindowed.';
    ERR(t) = d(t) - y(t);
    LMS(:,t+1) = LMS(:,t) + mu*Uwindowed.'*ERR(t);
end;

The inputs to it are u(t), the noisy signal and d(t) the desired signal.  For the standard LMS algorithm to work, you need to have an idea of what the noiseless signal looks like.

From your equations:

$\alpha$ appears to be similar to mu (or mu/2).
u is your s+n.
I am not sure what your d is.
Your equation for the error is, er, erroneous. :-)


Discrete Time Fourier Transform Pair Discrepancy When taking the DTFT of a discrete-time signal $x[n]=a \times u[n]$, the book "Signals and Systems" by Chi-Tsong Chen claims that the result is $\frac 1 {1-a \cdot e^{-j \cdot \omega \cdot T}}$; however, many other sources claim the DTFT is $\frac 1 {1-a \cdot e^{-j \cdot \omega}}$ without the T term.

The book I have states that $x[n]=x(nT)$ by definition and the DTFT by definition is

$$X(\omega)= \sum_{n=-\infty}^{+\infty}x(nT)e^{-j \omega n T}$$.

Can someone explain the reason for the differences?

Thank you for your time.

This is just normalized frequency versus un-normalized frequency.

Check the capitalization of omega. The capital omega $\Omega$, usually refers to the normalized frequency (with $-\pi$ and $\pi$ corresponding to $-f_s/2$ and $f_s/2$) and the normal omega $\omega$, to the regular angular frequency.

Your two sources just use different convention. There is no contradiction.


Circular buffer wrap around with chorus effect I am developing an embedded DSP audio processor and am trying to implement a chorus effect. To my understanding, the chorus effect is multiple delayed versions of the original signal, where the delay itself is random. In my program I have set the delays to be randomly defined upon initialization and then remain static during the program execution.

The challenge I'm faced with is that in using a circular buffer to store my past signal, I cannot simply wrap delayed signals back around the circular buffer as there is multiple delays present.

With a single delay I can simply check if the delay required exceeds the current position in the circular buffer, if so, I can jump to the end of the circular buffer subtract the delay and add my current position in the circular buffer.

if (Delay > Current_Position_in_Buffer)
{
    Circular_Array[Length_of_Array + Current_Position_in_Buffer - Delay];
}
else
{
    Circular_Array[Current_Position_in_Buffer - Delay];
}

Is there an efficient way to check and adjust the position in my circular array of each delayed signal, without separately checking each delay with an if statement?

If I understood correctly what you are asking, then this is probably more of a programming question... But the solution would be to use modular arithmetic to compute the position of the delay. With a length equal to a power of 2, this could be accomplished with an AND instruction.

const int Mask = (1 << 16) - 1;
Circular_Array[(Current_Position_in_Buffer - Delay) & Mask];

In this example, the length of your buffer is $2^{16}$.

More generally you can use the modulo operator, but it will be slower.


TTS - Diphone database I am working on text to speech conversion. I have done speech synthesis by

concatenation of phonemes.The transition between phonemes is not

smooth.I want to shift to diphone synthesis.  Where can I find diphone

database for the English language. Any other alternative might help.

Have you had a look at the diphone synthesizer resources in the Festival Speech Synthesis framework?

The state of the art for Text2Speech (though still not practical for production due to GPU compute requirements) is Deepmind's WaveNet.

Their results are definitely smoother than concatenative or parametric methods.  Good luck!


Why does SNR have an impact on the bit rate? I have some confusion regarding the impact of SNR on the bit rate. For example, if I have a QPSK-OFDM system, then my bit rate will be 2bits mutiplied by the number of sub-carriers of the OFDM modulation over the symbol time. How can I account for the SNR in this case?

Yes, try to understand what is the concept of SNR  and it depends on what. Then, it will be easier to understand why SNR and BER increase and decrease in all cases.

In your case, OFDM, the relationship between the SNR and BER is given:

$\textrm{BER} = \frac{1}{2} \left(1 - \sqrt{ \frac{L/N \times \textrm{SNR}}{2 + L/N \times \textrm{SNR}}} \right)$

Where $N$ is the number of subcarrier and $L$ is the channel diversity. $\textrm{SNR} = p/\sigma^2$ .. where $p$ is the transmitted signal power, and $\sigma$ is the variance.

so now, you check the relationship between the type of used modulation with the above equation, (QPSK and N) , it's easier to explore the behavior of BER. remember also that probability of error when using QPSK is different about when using for example QAM and so on

good luck


Why does Simulink generate this code for a PID controller? For the Simulink PID Controller model

The Simulink generated code (rewrite for better understanding) is:
#define PERIOD 0.005
double PID_step(double err,double P,double I,double D,double N){
  static double Integrator,Filter;
  double POut,IOut,DOut,Out;

  POut= err * P;
  DOut = ((err * D) - Filter) * N;//LineA
  Out = (POut + Integrator) + DOut;
  Integrator += (err * I) * PERIOD;
  Filter += PERIOD * DOut;//LineB
  return Out;
}

Line A and B are the code for the differential part, err is the input and DOut is the output. The Z-domain transfer function is $$\frac{DN}{1+\frac{NT_s}{z-1}}.$$
I want to derive the code from the transfer function.
Let $$\frac{Y(z)}{X(z)}=\frac{DN}{1+\frac{NT_s}{z-1}},$$  then $$zY(z)+(NT_s-1)Y(z)=DN(zX(z)-X(z))$$
Use Z-transform formula $$z^{-n}Y(z) \leftrightarrow  y[k-n],$$ we have $$y[k+1]=DN(x[k+1]-x[k])-(NT_s-1)y[k]$$
Here x is the input and y is the output.
So my code for the PID differential part is
double y0=0.0;
double Ts=0.005;
double PID_DifferentialPart(double x,double D,double N){
  static double x_last,y_last;
  static int first=1;
  double y;
  if(first){
    first=0;
    x_last=x;
    y=y0;
  } else {
    y=D*N*(x-x_last)-(N*Ts-1.0)*y_last;
  }
  y_last=y;
  return y;
}

Using the same input (e.g. input sequence {0 2 3 4 2.5}, with D=2.0,N=100.0,Ts=0.005), the outputs are different. The output of my code is {0 400 800 1200 1100}. Output of Simulink code is {0 400 400 400 -100}. Where did I go wrong?

I'm not sur exactly where you went wrong because your haven't fully explained your approach but for the derivate component of the PID, you should convert the C++ lines to this :
$$ D_{out}[n] = N(err[n]D  -F[n-1]) \\
F[n] = F[n-1] + T_s*D_{out}[n] \\\\$$
In the z-domain : $$
F(1-z^{-1}) = T_s*D_{out} \\
F = \frac{T_s*D_{out}}{1-z^{-1}}\\
D_{out} = N(err * D-Fz^{-1}) \\
D_{out} = N(err*D-\frac{T_sD_{out}}{1-z^{-1}}z^{-1}) \\
D_{out}(1+\frac{NT_sz^{-1}}{1-z^{-1}}) = N(err*D) \\
\frac{D_{out}}{err} = \frac{N*D}{1+\frac{NT_sz^{-1}}{1-z^{-1}}} \\
\frac{D_{out}}{err} = \frac{N*D}{1+\frac{NT_s}{z-1}}
$$
Edit :
You posted C++ code. Why is x_last only updated once? It needs to be updated every period. Try my solution below
double y0=0.0;
double Ts=0.005;
double PID_DifferentialPart(double x,double D,double N){
  static double x_last,y_last;
  static int first=1;
  double y;
  if(first){
    first=0;
    x_last = 0;
    y=y0;
  } else {
    y=D*N*(x-x_last)-(N*Ts-1.0)*y_last;
  }
  y_last=y;
  x_last=x;
  return y;
}


How to distinguish between the different frequency domains? Sometimes the terms 'Fourier domain', 'complex frequency domain', 'Frequency domain' and 's domain' are used interchangeably.

Take those answers here for example: http://answers.yahoo.com/question/index?qid=20090924230321AAxvnJg

Can you really use them interchangeably all the time, without being technically wrong? So, could you describe what would be wrong if I would replace 's domain' by 'fourier domain' for example? Or replacing 'complex frequency domain' by 'frequency domain'?

This answer is primarily opinion. It seems to me that the most sensible thing to do is to make the transform of interest clear in your text. Then you can use "frequency domain" exclusively and unambiguously.

For Laplace transforms, "s domain" will work as well, and I do believe they are interchangeable. I do not believe "Fourier domain" would be correct usage in this case. "s domain" may be preferred if you wish to emphasize that real values of $s$ have a clear aperiodic interpretation (exponentials in the time domain).

If the transform is not yet clear, you could say "Fourier domain", I suppose. In that case the assumed spectrum would be the complex (possibly conjugate symmetric) form.

Still, if it is within your control, be explicit about the chosen transform.


What is the state of the art for vocal elimination/reduction? I have no background in signal processing but I would like to know what is the actual state of the art for vocal elimination/reduction. Obviously, I am not talking about the infamous "stereo elimination", but other, more sophisticated techniques that not everybody is aware of.

I understand that for many projects there might be no software or online application available, but I'd find very useful to hear a demonstration.

Some recent ideas (2-3 years old, I haven't followed recently):

Design an elaborate generative model of the spectrogram of music signals, combining a non-negative mixture of a small set of background frames (models the background music) + a harmonic comb filtered by a non-negative mixture of smooth basis functions (the source filter model of the lead instrument), with constraints on the pitch transitions. Estimate model parameters from the audio clip. Use Wiener filtering to resynthesize the background music/vocals from the estimated spectra. See A musically motivated mid-level representation for pitch estimation and musical audio source separation by Durrieu et al. Such Bayesian approaches can be made more complex ad nauseam, by integrating prior knowledge about note transitions/durations - at the cost of increasingly hairy inference procedures. Demo time!
Techniques inspired by background subtraction in video analysis. The idea is that the accompaniment (especially in "bland" pop music) is invariant by a translation of a few bars (the duration of the "loop" or "riff") - while the vocals are not. Once the structure of the piece of music is known and can be described in terms of repeated groups of bars (a well-studied topic in itself), one can take the median of translated versions of the spectrogram to suppress the vocals. See REpeating Pattern Extraction Technique (REPET): A Simple Method for Music/Voice Separation by Rafil & Pardo. Demo time!

In the commercial space, you can download a trial for ADX Trax.


Are there algorithms to invert a magnitude spectrum? I have a speech signal and have calculated the magnitude spectrum for the signal. Is it possible to get a good estimate of the original signal from the magnitude spectrum, seeing that phase information has been discarded after the FFT?

I have looked at the Wiener-Khintchine theorem, but I am unable to find any code examples of it in use, and I have not managed to get it to work.

Apart from the stochastic literature, there are a set of deterministic algorithms that claim exact(?) recovery of a signal form its partial Fourier description, collectively referred to as:

  signal reconstruction from DTFT phase or magnitude alone

They are generally iterative in nature and exact convergence is hard to achieve in practice. Monson Hayes [et al] published relevant IEEE papers to propose several algorithms that reconstruct a real signal from its DTFT phase function alone, which is based on the phase sufficiency claim. The magnitude is a little more diffcult as I know. Have a look at those papers for your benefit.


Artefacts from playing ultrasound through PA system I encode data in audio by energising 5 frequencies between 17.5kHz and 19.5kHz:

I play the WAV on a sports stadium (20k seating, so BIG) PA rig, placing a microphone somewhere in the spectator seating area.

This is the data (I did it at vol=0.6 and vol=1):

As can be seen, there is some kind of 'ringing' occurring. Now the frequency channels are using are all integer multiples of a fundamental frequency (~100Hz), and the ghost frequencies appear at exact harmonics. If I crank up the volume I also get energy midway between harmonics.

There is one further effect: Ghosting at 1/2 f. This is most visible if I analyse a sweep (vol=0.8:

As can be seen, much of the energy at frequency f has been redirected to 1/2 f.

Can anyone explain these unwanted artefacts?

PS Often I observe a perfectly valid SE question attracts several instant down/close-votes before making it back to the surface. It always makes me scratch my head a little. It seems rather feeble to downvote without supplying a comment to summarise the decision logic.  e.g. Is it OT / malformed / inarticulate / lazy / etc, and how so?

PPS Someone has now suggested _intermodulation_ in a comment. Now this is valuable information. The question together with this keyword, once catalogued by Google, adds to the effectiveness of the Internet as a resource for Information. Yet the question is teetering on 4 (un--publicly-justified) close-votes. I do not think this is right. :/

First of all, like Cal-linux mention in their comment, I would strongly suggest intermodulation distortion for the frequencies that are not multiple integers of the fundamental.
Now, regarding the sweep tone artefacts. Since you are talking about high frequencies we are talking about a high-frequency driver. Nowadays, most high-frequency drivers are domes with some kind of compression driver attached to them. It has been proposed by Hubbard that compression drivers generate subharmonics. By visual inspection (as you also state) part of the energy of the sweep seems to be present at the first subharmonic $\frac{1}{2} f$, which could potentially be somehow related to the results Hubbard presented.
Further investigation is needed to be able to say for sure but this is indeed a very interesting observation. I work as a P.A. engineer, so I will try to duplicate the effect and work more on it. I will post back if I manage to reach any conclusions.


Why we sometime use some values of fixed point and floating point both in dsp? I have been writing kernels for a communication dsp device, which uses cmac/cmad operation on fixed and floating point numbers. I am not very clear why both types is used in the same operation.

The transition from text-book math to product code can be a lot simpler with floating point than with fixed point. The latter requires some specialized knowledge that, while not rocket science, may be more common for electronic engineers or embedded engineers than computer science or math people. Simple stuff tends to be a lot simpler in floating point than in fixed point.
If you really need to know the numerical behaviour of your code (as in certain recursive algorithms), it tends to be the other way around - the «cruise control» of floating point also makes it more tedious to figure out what is actually going on, while the explicit requirements of fixed point makes the details more plain sight.


Question About Kailath's Paper - An Innovations Approach to Least Squares Estimation Part I: Linear Filtering in Additive White Noise I'm reading the paper at the link below and I was following it for about 2 pages until I hit a road block on the bottom of page 648 where the author says:

  putting together 9-11, we obtain

and  gives  equation 12. I understand equations 9, 10 and 11  but don't see how 12 comes about from them.

Also, if anyone knows a book or a paper that provides an exact mathematical definition of white noise, it's appreciated. Definitely, I'm seeing how important it is to understand white noise, particularly in the context of this paper. Gaussanity ( which is what I'm used to seeing  in the many derivations of the Kalman filter ) is clearly quite different ( Dilip, I always believed you but now I'm seeing the importance )  and doesn't cut it here. This paper is actually the first one of a four part-er but, if I can follow this one by the end of 2018, I'd be quite satisfied.

An Innovations Approach to Least Squares Estimation Part I: Linear Filtering in Additive White Noise.

They used to have different convention for writing this stuff back then.
But actually what you saw is really simple.

It's all based on the Orthogonal Principle of MMSE.

They say in (9) that the Additive Noise is uncorrelated (By defining its Auto Correlation by Delat Function).
In (10) they say the optimal estimation of $ x \left( t \right) $ is given by $ \hat{x} \left( t \mid t \right) $ is a linear combination of $ y \left( t  \right) $ and $ z \left( t \right) $.
In (11) they exactly use the Orthogonal Principle which the estimator must obey.

In (12) they just derive the Correlation Matrix between the processes.
Since the linear estimator represent the correlation it is not surprising to see that it is the linear combination defined by the optimal filter $ g \left( \cdot \right) $.

The trick here is the fact (Using article notations) $ \overline{ \left( x \left( t \right) - \hat{x} \left( t \mid t \right) \right) v' \left( s \right) } = 0 $  which means (By Linearity of the Expectation) $ \overline{ x \left( t \right) v' \left( s \right) } = \overline{ \hat{x} \left( t \mid t \right) v' \left( s \right) } $.
Expand the right hand term by (10) by plugging $ \hat{x} \left( t \mid t \right) $ and you get (12).

By the way, Kaliath has much better written, in my opinion, paper on the subject called A View of Three Decades of Linear Filtering Theory.


Simplest All Digital GMSK Modulator What is the simplest way to implement an all digital GMSK modulator with no multipliers? The implementation approach should be applicable to GMSK modulation with BT=0.5 or less (where intentional ISI is introduced to obtain higher spectral efficiency).

By "all digital" I am referring to the baseband modulation. The result will be IQ samples that fully represent the GMSK modulation at baseband.

The acceptable answer will be the simplest implementation with no multipliers needed that will generate a GMSK modulated signal that has a correlation coefficient with a "perfect" GMSK signal > 0.999 and a spectrum that will meet the spectral mask for GSM as shown in the plot below (which has a symbol period of 48/13 μs (i.e., 270.8333 kHz), that I got from this website: http://etutorials.org/Mobile+devices/gprs+mobile+internet/Chapter+1+Introduction+to+the+GSM+System/Radio+Interface/:

Please preface your answer with spoiler notation by typing the following two characters first ">!"

I don't actually have an answer, just ideas, but since this is a riddle, I might of contribute these as (hidden as spoiler text) hints for others, which might have ideas that I miss:

So, MSK is often introduced as

    represented by a half-symbol time-offset BPSK on the I and Q components. That comes in very handy here (avoiding multiplicators) – BPSK can be implemened as $\pm 1$ on the respective axis. The idea is that at the boundary between symbols periods, the phase is continuous, which w.l.o.g cannot be realized in any other way than by having zeros exactly at the point where the other component carries the BPSK constellation point. Thus, the minimal viable MSK modulation has

   $$\begin{align} I&= [- {s_1} && 0 && -{s_2} && 0 && -{s_3}&\dots ]\\Q&= [ 0 && s_1 && 0 && s_2 && 0 & \dots ] \end{align}$$

   as samples, with $s_n$ being the bits to be transmitted represented as $\pm 1$ differentially. That doesn't look very much like CPFSK at all, but one has to realize that the advance between two consecutive samples takes two different values:

   For samples that have an even index $n$, the following sample always is -90° further, no matter which value $s_{n/2}$ has

   For samples that have an odd index $n$, the phase advance is -90° if $s_{\lfloor n/2\rfloor}=s_{(n+1)/2}$, and +90° in case the bits differ.

   Considering the phase difference of consecutive samples as the instantaneous frequency, we thus see that we get an alternation of a constant $-\frac \pi2$ between samples and either $0$ or $-\frac \pi2$, carrying the information. I'm not even convinced any more this fulfills the criteria for being MSK.

   Since we'd probably want spectral symmetry, we'd employ the usual tricks to shift the signal in frequency by half the sampling rate (i.e. "multiplication" with $[-1\,+1\,-1\,+1\,\dots]$).

Oooh and I'd have an approach for the Gaussian pulse shaper, but that simply applies

   the fact that repeated convolution of something with itself converges against a Gaussian – same reason the CLT works.

 That means I can just cascade a set of moving average filters – all of which work without multipliers to get an approximated Gauss filter.

Sadly, that doesn't inherently solve the issue

   that my MSK approach above basically only allows 1 sample per symbol – but we might at least recreate the inter-symbol dependence characteristic of GMSK trellises.


Gonzales Digital Image Processing and the Laplacian Operator In Gonzalez' book Digital Image Processing, third edition, he writes on page 182 that

  It can be shown that the simplest isotropic derivative operator is the laplacian operator

where the Laplacian is given by
$\nabla^2f = {\frac {\partial^2 f }{\partial x^2}} + {\frac {\partial^2 f }{\partial y^2}}$, and the particular paper he links in reference is not available for free.

It is :
Rosenfeld A, Kak A.C. Digital Picture Processing, vols. 1 and 2, 2nd ed., Academic Press, New York.

I was looking for insights to why:

a)The Laplacian is isotropic

b) What he means by "simplest". I am assuming this means in terms of computational complexity of things like how many differentiation operations, additions ops, multiplication ops, etc.

Any help appreciated.

Isotropic operator is an operator that has no bias to any particular direction in the image.
If you see the below Laplacian mask, Edge detection is taking place in every direction.

If you see the below image, three filter kernels are applied. The first kernel detects edges in x direction, the second in y direction and the third is Laplacian kernel detecting edges in all directions, hence it is isotropic.


What distance metric can I use for comparing image features like elongation and solidity of each image? What distance metric can I use for comparing image features like elongation and solidity of a contour of each image?
Except Least Square and without using a support vector machine because i do not know at which class images belong.

If I understood you correctly, each contour is described by a 2-element feature vector $f = [e, s]$, where $e$ is elongation and $s$ is solidity.

In that case, you might want to try the Mahalanobis distance, which is defined as follows:
$$d(f_1, f_2) = \sqrt{(f_2 - f_1)C^{-1}(f_2 - f_1)}$$
where $f_1$ and $f_2$ are the feature vectors that you are comparing, and $C$ is the covariance matrix of your data set.


I/Q Sample demodulation to input data I am trying to convert some raw unscaled I/Q sample data back to the original input data. I've been reading a lot of different posts and articles in I/Q demodulation and the focus is almost solely on the process of combining the two separate I and Q pieces on the transmitter side into one and then separating them again at the receiver side. I haven't found any helpful documentation on how to do the final piece of this.

A lot of programs exist that output data directly in an I/Q raw format. In my case, I have 16-bit signed integers in the format of $i[0], q[0], i[1], q[1], ...$ . Presumably, the signal originated from a binary output on the transmitter end. After the first demodulation part, the separation of the I and Q values on the receiver end, what final steps are needed to complete the cycle back to that input binary value?

I realize this is most likely extremely simple, especially after all of the math that had to be done to get the two waveforms back from the single RF one. Nevertheless, I have yet to find a decent explanation of this final process.

In addition to hotpaw2's answer:

  Presumably, the signal originated from a binary output on the transmitter end.

I'd say that is most probably wrong.

The typical digital transceiver based on complex baseband looks like this:

Binary data -(0/1)->
 Symbol mapper =(complex-valued symbols)=>
  pulse shaper =(complex-valued, "smooth" digital baseband signal)=>
   synchronous dual-DAC =(I & Q signals)=>
    filters =>
     I: mix with cos(2 pi f_RF) --\
                                   subtract--> filters -(RF signal)-> Antenna
     Q: mix with sin(2 pi f_RF) --/

On the transmitter end, you'll see some RF signal – and that is always analog, never binary.

  I haven't found any helpful documentation on how to do the final piece of this.

You cannot have read very much, because every textbook I know introduces both, the direct conversion transmitter and the direct conversion receiver (note: these go by different names, also zero-if or quadrature receivers).

Since you're already getting complex samples, we may assume you're working with a receiver that works like this

Antenna -(RF signal)->
 band-pass filter ->
  multiply with cos(2 pi f_RF) -> low-pass filter -(I signal) ->
  multiply with sin(2 pi f_RF) -> low-pass filter -(Q signal) ->
   =(IQ signal pair)=> synchronous dual ADC =(digital complex baseband signal)=>
    PC software

So, the analog part has been done for you.

Now, it's up to your PC software to do (not necessarily in this order, but all of these are typically necessary):

estimation/correction of the difference in f_RF (no two physical oscillators are perfectly identical)
estimation/correction of the difference in phase (i.e. the cos() of receiver and transmitter didn't "start" at the same time, and the signal was "in the air" for some non-zero amount of time due to limited speed of light)
reversing the pulse shaping of the transmitter
estimation/correction of the symbol rate (ie. how many digital samples is a single symbol long?)
guessing which symbols got sent by the transmitter
converting these symbols back to the bits they stood for

In reality, there's a few additional steps (channel coding, hardware imperfection corrections, equalizing, …).

Anyway, the order of these steps, and how exactly they're done depends completely on what kind of modulation scheme/waveform the transmitter produces; notice how you can't use a 3G mobile phone to get the bits your cordless landline phone sends, or use your WiFi card to receive digital TV, or use your CCTV receiver to receive satellite telephony.

Even for receivers for the same kind of signal, there's a lot of choices a SDR engineer makes – the type of timing recovery mechanism, for example, acceptable computational complexity when converting samples to symbols vs error rate when doing that, and a lot, lot more.

Thus,

  I realize this is most likely extremely simple,

Nope. It's the hardest problem, and pretty much the main problem, of  digital communications technology, the science backing one of the biggest technical sectors of modern days. Building good receivers is hard, in general. There's a lot of tools to help you deal with that. But: even these tools require solid understanding of the theory behind digital communications before you can go out and write your own SDR receiver.

But don't despair:

There's a lot of "simple" systems that are relatively easy to receive, and a tutorial alone might get you pretty far. It's an exciting field of technology, and I wouldn't have spent a solid  part of my life on that if it wasn't fun!

I've witnessed literally dozens of times how being able to actually get binary data out of something that is just a digital representation of an analog signal makes people a) happy and b) makes the understanding of the theory take giant leaps.

I'll encourage you to experiment with digital signals!

Based on your assumption, the transmitted signal was binary, I'll go ahead and take an audacious guess at what your profession is: Is it possible you are a Computer Science student?

I can see this signal world might be a bit daunting if you're coming from a background where you can safely assume that for

input---wire/file/network-->output

input==output, but that's simply not the case for digital communications.

In any case, no matter what you are, I'd recommend reading a bit on signal theory. The GNU Radio wiki has a page on recommended literature. Get one of the basic DSP books, and one of the basic Digital Comm books; some of them are freely available as PDF.

Parallely, take a 5-minute read of the introduction chapter of the GNU Radio guided Tutorials. These tutorials are meant to slowly, and carefully, and sequentially (don't skip a chapter) allow you to build a simple digital receiver. What they do is teach a tool, but they need to rely on some signal theoretic background, so they might be not that perfectly clear if you lack that.


Finding the gradient of a norm in a minimization problem I have to find the gradient of the following term with respect to $X_{1}$:

$\|\Phi\circ(X_{1}-X_{2})-u\|_F^2$ ,

where $u\in\mathbb{R}^{n}$; $X_{1}, X_{2}\in\mathbb{R}^{N\times J}$ and $\Phi\in\mathbb{R}^{n\times NJ}$. The $\circ$ operation is defined as follows:

$\Phi\circ X=\Sigma_{i=1}^{J}\Phi_{i} x_{i}$,

where $x_i$ are the columns of $X$ and $\Phi_{i}$ are $n\times N$ submatrices of $\Phi$.

I know that $\frac{\partial}{\partial X}\|X\|_{F}^{2}=2X$, but the $\circ$ operation is making it rather complicated. Can somebody help?

P.S. Also, would changing the Frobenius norm to l2 norm make any difference here (since $(\Phi\circ(X_{1}-X_{2})-u)$ is a column vector)?

This is the result of taking the derivative of $\|\Phi\circ(X_{1}-X_{2})-u\|_F^2$ w.r.t. $X_{1}$:

$[\Phi^{1'}(\Phi\circ(X_{1}-X_{2})-u), ...,\Phi^{J'}(\Phi\circ(X_{1}-X_{2})-u)]$ .

This can be got by first finding the expression for the norm, and then taking the matrix derivative (w.r.t. $X_{1}$) to get an $N \times J$ matrix. This matrix is then represented in terms of the sub matrices $\Phi^{i}$.


Implementing wavelet transform for finding transients in the power supply I am new to the concept of wavelet transforms. Can somebody please help me in understanding this ? and also how to implement it in c. Is Short term Fourier transform more efficient than Wavelet Transform for finding Transients ?

I would say that a matching Mother wavelet could be the best for detecting a transient; but both the selection and implementation would be much slower. The old adage: do you want quantity or quality :) or Bandwidth vs. noise.  Life doesn't come in our neat intellectual packets.
BTW: the simplest technique is an isolation capacitor and a threshold detector (sad to say I know a lot of ways to "cheat"  ).

Evaluation process:

I would look to octave or scilab for a simulator to use.
Construct a model of transients you think are likely or that you are looking for.
Pass the transients through various canned analysis routines and then
Find the transients in the output; the crisper the "find" the more applicable the analysis is.

Reconstruct the input and see how good the reconstruction is; i.e. apply an error criterion to achieve a faithfulness measure.

If this seems tedious it's because it can be: but at each point in the process you get an idea about the effort involved.

Selection time

Design time and complexity

Execution time

Accuracy.

And most importantly:

Knowledge!

Like the wavelet theory: start coarsely/crudely and refine the evaluation process.

I looked on google and found:
----http://www.soest.hawaii.edu/MET/Faculty/bwang/bw/paper/wang45.pdf
A real application of Waveform/Wavelet analysis of El nino and such; with all the hair:
--https://inst.eecs.berkeley.edu/~ee225b/sp14/lectures/shorterm.pdf
Which has a description of the process and history.

--http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.506.6298&rep=rep1&type=pdf
Which discusses real usage ambiguity.

--Why Wavelet developed when we already had Short-time Fourier transform
Which links to another discussion.

--https://en.wikipedia.org/wiki/Wavelet

Back to question:
Most of the techniques are "complete" in the sense that the transforms are invertable and can recover the original input.  In that sense all invertable transforms are equivalent.  But like I said they isolate different features.  Wavelet transforming is a mathematical process for utilizing Mother wavelets to isolate information in time and space. STFT is very specific whereas there are innumerable of "wavelet transforms" each having different sampling Mother wavelet shapes that are applied recursively to give many (time, freq) graphs; two dimensional pictures.
STFT and Gabor and such are (as far as I know) specific but closely correspond to normal fourier/z-transform analysis.

There are a lot of c wavelet programs in the open source domain.

An effort to clarify "disturbance" analysis.

First you have to  be able to define what a disturbance is and then the allowable type one/type two results; i.e. false positives/negatives and then do the system design around this.

Let me list a few "disturbances"; not necessarily related to power supplies.

For a nuclear containment vessel we tried to "listen" for what amounts to tearing at 100khz-1Mhz with transducers.  A time exponential increase in the power of the signal was taken as (and experimentally verified) an indication of a problem.

As far as power supplies go:
There are safety related issues such as over-voltage causing downstream failures.  For these I always put in absolute hardware; zeners, surge suppressors.  That I think will limit damage and fail-safe; shut things down.  Do the thermal calculations with respect to the calories produced by hardware shunts.

Under-voltage is also a problem in some cases, as I found out by almost smashing a co-workers hand.  Dropping the field voltage to a DC shunt drive motor causes it to speed up :)  So "pulling the plug" at random is not a great idea; and a drop in field voltage should probably be acted on.

Other expected cases are parasitic resonances/oscillations occurring in the output power stage of supplies.  These are typically produced by altered frequency response characteristics of bipolar output transistors during transitions.  Typically these are real transients and not worrisome unless some downstream equipment is sensitive.

You have specification errors where for some reason, say a passing cyclotron, the system is substantially normal but some specified tolerance is breached.  This can also happen slowly when equipment ages and drifts out of spec.

Of course in motor-generator situations transient signals might indicate mechanical wear.

Then you have fast raw supply transients that find their way through the regulators and occur at the output.  That's the reason I use hardware suppresors/zeners that can absorb some over-voltage attempts without allow ing the disturbance through and recover easily.


what is OQAM (Offset QAM) and why it is required? Currently I am studying about a modulation technique FBMC (Filter Bank Multicarrier) which is a predecessor of OFDM. The difference with OFDM is that, FBMC does not have a CP (Cyclic Prefix) which is a better use of spectrum. On the other hand, it is susceptible to ISI (Inter Symbol Interference) due to the absence of GI (Guard Interval) when transmitted over wireless channel.1. FBMC uses prototype filters to shape the out of band emission to lower the interference.
There is a concept of FBMC-OQAM (sometimes referred as staggered multitone or SMT), where the input signal is first converted to QAM symbols and further processed into OQAM symbols[2].
Can anyone explain to me in very simple words / general context:
#. What is the purpose of OQAM (ref.2 has a figure) and how it is making the signal orthogonal?

1https://ieeexplore.ieee.org/document/5753092
[2]https://ieeexplore.ieee.org/abstract/document/8255374

Actually, in theory and error performance there is no difference. But in practice, OQPSK  decrases the maximum phase alteriation from 180 to 90 degrees which is better for RF power amplifiers.


FFT-like algorithm for fast DTFT computation? Good morning!

I'm coding up a project on a microcontroller to read in some analog audio (specifically, the sound of someone whistling: a near perfect sine wave) and determine which piano note tones are present in the signal.

Although the FFT is the obvious first choice, the frequencies I wish to detect in the signal don't line up with the FFT frequencies. Furthermore, I only want to detect around 10-15 frequencies, not the whole spectrum.

So my thoughts quickly turn to the DTFT, which can be implemented with a matrix multiplication (here is my prototype MATLAB code. Details provided on request). However, surprisingly, it still takes more computation than an entire FFT!

This brings me to my question. I've made all the obvious google searches ("fast DTFT", "FFT for specific frequencies") without much luck. Does anyone here know of an efficient DTFT method?

According to Oppenheim and Schafer's "Discrete Time Signal Processing", the Goertzel algorithm will be more efficient than the FFT in computing an N point DFT if less than $2 Log_2 N $ DFT coefficients are needed.  So if you need to compute 15 frequency points, the Goertzel will be more efficient for total number of samples $N > 65536$.  Anything below that and the FFT will be more efficient than Goertzel as Marcus has implied.

As you probably discovered, one point in the DFT requires N phase rotations and N complex additions (assuming a complex input).  (As the DFT output is computed by simply multiplying your input signal through successive rotations and accumulating the result-- in this manner you are correlating to each output frequency of interest).

Each phase rotation requires 4 multiplications and 2 additions (if you have time to iterate with little multipliers, see the CORDIC algorithm as an alternative phase rotator that requires no multipliers).  The phase rotations that are +/-j are simplified, since a rotation by j is done by simply changing the sign of Q and swapping I and Q, but for a larger FFT this will be an insignificant number of the total computation.

So therefore for each point using the DFT, there are approximately 4N total multiplications and 4N total additions.  So for an N pt DFT that would be $4N^2$ of each - wow!

In comparison the FFT requires $2Nlog_2N$ multiplications and $2Nlog_2N$ additions. When considering a complete DFT, the FFT offers a dramatic reduction in required computations. (Pause here for some respectful bowing to Cooley-Tukey).

I would be interested in seeing how the Goertzel algorithm lands on the graph below but haven't muddled through those details yet.

Regarding your other challenge mentioned that your frequencies do not line up; for that I recommend windowing first to reduce the effect of distant frequencies on each tone you are measuring and then interpolate between the bins to more accurately estimate your frequency result for each tone. (A simple linear interpolation may be sufficient).  The windowing will reduce the frequency resolution of each bin but will also reduce the sidelobes from the other bins, so you would need to confirm that your frequency separation for your tones of interest are sufficiently larger than the equivalent bandwidth of your window+FFT. For a good summary of equivalent bandwidths vs window refer to fred harris' paper on windowing: https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf


Should I filter before calculating power, if I'm looking for relative changes in power? I am building a "scanner" sort of program, which reads IQ samples from an rtl-sdr stick, and keeps track of the power on certain frequencies, in certain bands.  I am interested in detecting activity on frequencies, which I do by looking at changes in power.  Since I'm simply calculating changes over time, can I get away with not filtering the data before calculating power, or is there some reason I should filter?

If the output from the SDR matches the band you're interested in, I see no point in filtering.

You want to filter when you have a scenario similar to this: you want to track the power in a band of $B$ Hz around a frequency $f_c$, but the SDR's output covers a band that is larger than $B$. In this case you want to filter, since otherwise a signal outside the band of interest could affect your calculations.


sound source localization - calculation of sound intensity vector angles I'm working on an acoustic source separation problem where I want to separate the voices of a choir based on intensity vector statistics.
I'm currently implementing a paper by Günel called “Acoustic Source Separation of Convolutive Mixtures Based on Intensity Vector Statistics” In this paper, the acoustic source separation is achieved in a couple of steps:

the signals from the microphone-array are converted to a time-frequency representation (using the modified discrete cosine transform, MDCT)
The intensity vector direction is calculated and rounded to the nearest degree. All of the obtained directions are then plotted in a histogram to find the mixture probability density. On this histogram, a von Mises distribution is fitted for every source position. (my question is related to this step)
As I understand it, this distribution is then used to define a directivity function for each sound source and a time-frequency bin and beamforming is applied to separate each source.
Finally, the inverse MDCT is calculated for each separate signal, obtaining the separate signals in the time-domain.

So this brings me to my question. After calculating the vector intensity directions for a single source in a simulated anechoic room (pyroom acoustics) and real-life measurements with a GRASS 50VI-1 Intensity probe, all obtained directions are plotted on a histogram and here, something I can’t explain happens. For a source that should be positioned at 45 degrees relative to the microphone array, the histogram shows two peaks. One peak at 45 degrees, but also one peak at 225 degrees (so, 180 degrees shifted).

This is simulated in a 3D space, but for now we are only using 4 microphones to find the source directions, so only finding the azimuth directions. These directions are found using the following formula (in python, we use the arctan2 operator):

where Pw, Py and Px are found by using our microhpone measurements as in the picture below:

I have read countless papers about DAO/separation using sound intensity, but not a single one I found shows any sign of this behaviour or talks about certain assumptions that should be made. It seems they only have a peak in the direction that the sound is actually coming from.
Finally, here is the code I used:

So this technique approximates the velocity of the particles using the pressure gradient at a time instance. This approach is only correct for microphones that are spaced closely together. The max frequency it can measure correctly with the spacer we're using is around 4.8kHz (stated in datasheet), so all angles corresponding to higher frequencies gave a wrong answer. This caused the 180 degree shifted angles.
A great explanation/visualisation from Siemens can be found here


Separate sub-figures! Any idea on separating multi-panel figures?

There are 2 different sub-figures connected each other.

Thanks!

Conceptually, it's straightforward. Let's consider a more general case, involving a grid:

First run ridge detection, box blurring, then line detection:

The purpose of this preprocessing is to make line detection more successful. Finally, partition the image vertically (only first partition shown):

then horizontally:

In Mathematica it can be done as follows:

img = Import["http://englishwalle.files.wordpress.com/2013/08/9-image-grid-no-logo1.jpg"]
{W, H} = ImageDimensions[img]
edges = ImageConvolve[Binarize@RidgeFilter[img, 0.001], BoxMatrix[2]]
lines = ImageLines[edges, 2, 0.1]
hlines = Select[lines, Norm[#[[1]] - #[[2]]] == W &]
offsets = Flatten@{0, #[[1 , 2]] & /@ hlines, ImageDimensions[img][[2]]}
hparts = Partition[H - # & /@ offsets // Sort, 2, 1]
vsections = ImageTake[img, #] & /@ hparts

(* Part 2 *)
hmask = ImageAdjust[ImageAdjust@RidgeFilter[Blur[#, {100, 0}], 5], 1] & /@ vsections
lines = ImageLines[#, 0.075, 0.25] & /@ hmask
vlines = #[[1, 1]] & /@ # & /@ lines
vlinesUnique = W*DeleteDuplicates@MeanShift[#, 1] & /@ vlines/W
vparts = Partition[Flatten[{0, Sort@#, W}], 2, 1] & /@ vlines_uniq
hsections = Outer[ImageTake, {#[[1]]}, {{0, W}}, #[[2]], 1] & /@ Transpose[{vsections, vparts}]

I leave it to you to translate this into your language of choice.

edit: With your image, I had better results with a derivative filter than a ridge filter (though that worked too, but the seam was not exactly at the right place). This is what directional blurring and derivative filtering (hmask = DerivativeFilter[Blur[img, {100, 0}], {0, 1}]) looks like on your image:

Line detection works trivially after this.


Result of inverse FFT is sometimes shifted in real space I am using the Numpy fft2, ifft2, and related functions and I am sometimes running into a strange situation where the output after performing the inverse Fourier transform the results are shifted. I know about fftshift and ifftshift, however my understanding of those is that they are used after fft2 and before ifft2 if it is needed for the DC component to be in the middle of the result instead of at one end. In particular, they are not supposed to be needed in real-space.

One example is in a demonstration of the convolution theorem. The following two things should be identical:

$$F^{-1}(F(im) \cdot circle)$$
$$im \otimes F^{-1}(circle)$$

where $F(x)$ is the Fourier transform, $F^{-1}(x)$ is the inverse transform, $\cdot$ is pointwise multiplication and $\otimes$ is convolution, $im$ is the image, and $circle$ is a simple, ideal, low-pass filter. The code I use to compare these is as follows:

# The fourier duck is a 256x256 8-bit grayscale image
im = skimage.io.imread('fourier_duck.png') / 255

# Create the ideal low-pass filter
x,y = numpy.meshgrid(range(-128,128), range(-128,128))
circle  = numpy.fft.ifftshift(x*x + y*y < 15*15)

# Perform filtering in Fourier space
f_im = numpy.fft.fft2(im)
f_im_circle = f_im * circle
finv_f_im_circle = numpy.real_if_close(numpy.fft.ifft2(f_im_circle), 10000)

# Perform filtering in real space
finv_circle = numpy.real_if_close(numpy.fft.ifft2(circle))
conv_im_finv_circle = scipy.ndimage.convolve(im, finv_circle, mode='wrap') # WARNING: takes lots of memory

The two outputs, finv_f_im_circle and conv_im_finv_circle should be quite similar but are quite different:

If the final, real-space, version is put through numpy.fft.fftshift it looks how I would expect (even though the plots don't show it, the value ranges are essentially the same as well, within FP error). The problem seems to crop up with the inverse Fourier transform of the circle which produces a kernel that also looks shifted.

I have tried doing shifts in other places and looking at the imaginary parts all to no available.

Besides the simple example I show above this has also occured with other filters designed in Fourier space, sometimes when doing inverse filtering, and even with filters designed in real space but applied in Fourier space. I can't figure out the pattern of when to shift the final result or when not to.

Ultimately my question is two fold:

Why is this happening?
How do I fix it? (which may just be how do I predict when I need to use the extra real-space numpy.fft.fftshift)

If you have a filter in the time/space domain, the location that is associated with zero linear phase in the frequency domain is the top left pixel of the filter. If the filter is centered in time/space, then you need to ifftshift prior to taking the FFT in order to remove the linear phase that centers it in the time domain. The following example Octave code (sorry, I don't use python) illustrates this:

x = zeros(10,10);   % Make an impulse in time, so that its
x(6,6) = 1.0;       % magnitude in frequency is all ones.

% Display the angle of the FFT both with
% and without the time-domain quad-swapping.
figure;
imagesc( angle( fft2(x) ) );
title( 'Phase Without Quad-Swapping In Time' );

figure;
imagesc( angle( fft2( ifftshift(x) ) ) );
title( 'Phase With Quad-Swapping In Time' );

The same thing is true in the frequency domain, which is what is happening with your filter. Specifically, you have defined the magnitude of the filter in the frequency domain, but you have not given it the phase necessary to center it in the time/space domain.

I tried to answer a related question here.


How to train and test deep neural network using MFCC features? I am working on Voice Disorder Detection problem. I have extracted 13 MFCC features, 13 delta and 13 delta-delta features from each audio file (2 to 4 secs). I extracted these features for each frame (consisting of 23 ms). So for example an audio divided into 75 frames has 75 examples in my dataset i.e each frame act as individual example in the dataset.

I trained Deep Neural Network on this dataset (70% of dataset) and evaluated my model on the rest 30% of the dataset giving me accuracy of 96% and 84% respectively.

But the problem is, while testing my model on a single unseen audio data when I calculate MFCCs frame by frame and feed every frame to the model, the model prediction for every frame is not the same. So my question is:

What should I do so that my model will give consistent prediction for all frames belonging to the same audio file?

I tried taking the majority of the class predicted by my model for that audio but that too gave me incorrect result.

One thing that I think I might have done wrong is that I haven't used every frame belonging to a particular audio file for training and so some of them are in test dataset. Is that a problem or am I missing something important here?

You may want to try training your model on 2D arrays containing entire datasets (all 75 MFCC frames of each audio file).  You can try recurrent or convolutional early layers in your DNN to handle the 2 dimensional input shape.  You may need a far larger amount of training data.


MATLAB - remove the frequency at zero in FFT I am getting a high frequency at 0Hz. How can I remove it? I am not supposed to get it.

My code is

W = 3.64*2*pi;

Fs = 50;                      % Sampling frequency in hz
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
v = 113*(1+0.40*sin(W*t));

theta0 = 12*(pi/180);         % in radians for 12 degree
theta1 = 6*(pi/180);          % in radians for 6 degree
theta = theta0-theta1*sin(W*t);

Lift = 0.5*1.225*v.^(2)*(2*pi).*theta;

NFFT = 2^nextpow2(L);         % Next power of 2 from length of Lift
Y = fft(Lift,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

subplot(2,1,1);
plot(Fs*t,Lift);

subplot(2,1,2)
plot(f, 2*abs(Y(1:NFFT/2+1)))

Frequency bin at zero is simply mean value of your signal. Just take a look on definition of DFT, for zero frequency $k$ we get:

$$\left. X[k]=\sum_{n=0}^{N-1}x[n]e^{-i 2\pi  n\cdot k} \right |_{k=0} \Rightarrow X[0]=\sum_{n=0}^{N-1}x[n]e^{-i 2\pi  n\cdot 0} = \sum_{n=0}^{N-1}x[n] $$

So after you calculate the FFT and divide by number of samples, indeed you get the mean value of your signal.
In order to remove that, please apply the following before the FFT calculation:

Lift = Lift - mean(Lift);


How to check the FFT results of a sine wave? I have been given an audio file (sine wave) 1000Hz as an input to my FFT algorithm.
I have got 8192 power spectrum samples in a array.

What is the best and easiest way to check whether my output is right or wrong?

If I give a silent audio file then the output is zero for all samples.
In a sine wav the o/p increases from 20 (0th sample) to 26059811 (743rd sample) and decreases gradually to 40.

If I get an idea of the output range then I can prove technically whether the FFT is working.

Any ideas would be helpful.

Refer to this link for any technical doubts.

You need to plot the magnitude of the output of the FFT.  I'm not familiar with your programming language, but in Python you would use something like plot(abs(fft(a))).  For a silent input, the output should be all zeros.  For a sine wave input, you should see two spikes:

For a real signal, the spikes will be symmetrical from left to right.  If you're doing a real FFT, though (which is more computationally efficient) you'll only get the left half of the plot as your output, since it ignores the redundant mirror image.

If the frequency is higher, the spikes will be closer to the center.  If the frequency is perfectly in sync with the chunk size, the spike will only be one point wide and everything else will be exactly 0.  Otherwise it will have a tapering "skirt" like above.


Facing problem in plotting histogram of DCT coefficient of an image In digital image forensic, it is assumed that the histogram of DCT coefficients of the un-compressed natural image has the smooth histogram. However, the histogram of DCT coefficients of the compressed natural image has peaks and gap in it. One such problem is discussed in this work. see this image for a better understanding .

I am using the following MATLAB code to plot the histogram of DCT coefficients of uncompressed and compressed images.

x_axis = -1000:0.01:1000;

img_1 = rgb2gray(imread('a.tif'));
img_2 = rgb2gray(imread('b.tif'));

subplot(2,2,1), imshow(img_1), title('Image');
subplot(2,2,2), imshow(img_2), title('Compressed Image');
subplot(2,2,3), plot(x_axis, histogram(dct2(double(img_1)))), title('Image');
subplot(2,2,4), plot(x_axis, histogram(dct2(double(img_2)))), title('Compressed Image');

but I am unable to get the results in the desired manner. see this figure below .

Can anyone please explain? what is the correct way of plotting the histogram of DCT coefficient of an image so that I may get the correct results?

My guess is that the assertion

  "the histogram of DCT coefficients of the un-compressed natural image
  has the smooth histogram ... However, the histogram of DCT
  coefficients of the compressed natural image has peaks and gap in it."

is valid when one uses a compression method that uses a DCT (Discrete Cosine Transform), and quantizes (sufficiently) its coefficients. This might not happen for compression methods that differ from JPEG/MPEG principles, that use a 2D DCT. Moreover, JPEG applies DCT2 on $8\times 8$ blocks, and the DC coefficient undergoes a predictive compression scheme different than that of the other AC coefficients.

So I am not surprised that you don't observe quantization when you apply the DCT on the whole image instead of blocks (I obtained the same results).

What could be more successful would be to compare the DCT coefficient histograms (DC coefficient removed) on each of the  $8\times 8$ blocks, raw and compressed, and then combining the coefficients of all the blocks. Here is an attempt, in Matlab.

As you can see, you now have holes in the right-side histograms, due to quantization.

nBin = 512;
dctCoeffLimit = 80;

img_1 = double(rgb2gray(imread('lena.bmp')));
img_2 = double(rgb2gray(imread('lena.jpg')));

% Compute block DCT with mean removed
funDCT2 = @(block_struct) dct2(block_struct.data-mean(mean((block_struct.data))));
img_1DCT= (blockproc(img_1,[8 8],funDCT2));
img_2DCT= (blockproc(img_2,[8 8],funDCT2));

% Set DC coefficients to NaN
img_1DCT(1:8:end,1:8:end) = nan;
img_2DCT(1:8:end,1:8:end) = nan;

% Convert 2D to 1D
img_1DCTflat = img_1DCT(:);
img_2DCTflat = img_2DCT(:);

% Remove DC (NaN) values
img_1DCTflat(isnan(img_1DCTflat)) = [];
img_2DCTflat(isnan(img_2DCTflat)) = [];

subplot(1,2,1)
histogram(img_1DCTflat,nBin);axis([-dctCoeffLimit dctCoeffLimit 0 Inf]);
subplot(1,2,2)
histogram(img_2DCTflat,nBin);axis([-dctCoeffLimit dctCoeffLimit 0 Inf]);


scale invariant face similarity metric I need some face similarity metric.
I know about PCA,LDA approaches that uses raw pixels data, but they are not scale invariant(I know that I can resize image but I think that's not good solution).
Maybe I need to use AAM,ASM aproaches?
But I'm interested in metric that gives me the greatest separability of faces.
At first I don't need automatic detection of metric features, for example imagine that I can mark all faces by hand.

Haven't used them first hand but I'm pretty sure that [SIFT] "Scale Invariant Feature Transform" is scale invariant. And I know that they are used, with great success, in face recognition.


how to set Equalizer's coefficient using generalized eigenvector. In that paper https://jwcn-eurasipjournals.springeropen.com/articles/10.1186/1687-1499-2012-72 .. Equation 3,

$J(w) = w^HBw/w^HCw$ ..

$B$ and $C$ are matrices, $w$ the filter coefficients vector we need to design.

The author is supposed to maximize that cost function $J(w)$. So, he used the generalized eigenvalue decomposition of the two equations $(B,C)$ then he set the eigenvector correspondent to maximum eigenvalue as equalizer filter $w$.

My question, Is my understanding for that idea correct? I mean we can use generalized eigenvalue decomposition for that purpose maximize/minimize, which means if set  setting the eigenvector correspondent to maximum eigenvalue as equalizer filter $w$ will maximize, and setting the eigenvector correspondent to minimum eigenvalue as equalizer filter $w$ will minimize.

thank you.

The generalized eigenvalue problem is given by

$$Bw=\lambda Cw\tag{1}$$

where $\lambda$ is the generalized eigenvalue of the matrices $B$ and $C$. Multiplying $(1)$ from the left with $w^H$ (with $^H$ denoting the Hermitian conjugate) and dividing both sides by $w^HCw$ (assuming that this term is non-zero), we obtain

$$\frac{w^HBw}{w^HCw}=J(w)=\lambda\tag{2}$$

This shows that the value of the objective function $J(w)$ (which is a generalized Rayleigh quotient) equals the generalized eigenvalue $\lambda$. Consequently, $J(w)$ is maximized by the eigenvector corresponding to the maximum eigenvalue, and it is minimized by the eigenvector corresponding to the minimum eigenvalue.


Least Mean Squares Algorithm confusion with Adaptive Line Enhancement I'm just a bit confused about the least mean squares algorithm to separate wideband and narrowband in an adaptive filter for voice conversation. I'm interested in the narrowband part and I'm confused about the LMS equation as follows:

Is h being iterated twice here with j and n? Wouldn't that mean I have two different sets of h values?

I'm having trouble putting this into MATLAB.

I recommend you read up on the LMS algorithm and try to understand it before you start implementing it, otherwise you won't be able to find any errors in your code. In that formula, the index $n$ is the iteration number, and the index $j$ is the vector index of the filter coefficients. So at any iteration $n$ you have $N+1$ filter coefficients (indexed by $j$) that need to be updated according to that equation.


Using Goertzel to find phase shift I have two signals with nearly the same frequency, and i would like to find the phase shift between them, since i am using micro-controller, i would like to do it in the less expensive way.

I know how to calculate Goertzel for a specific frequency bin (for each signal), but how do i calculate the angle for this bin by using Goertzel?

If you know a lower computational method, please advice.

Thanks

If the complex frequency bin outputs of the two Goertzel filters are $a$ and $b$, you can calculate the phases as $\alpha = \text{atan2}\left(\text{Im}(a),\ \text{Re}(a)\right)$ and $\beta = \text{atan2}\left(\text{Im}(b),\ \text{Re}(b)\right)$, where $\text{Re}$ returns the real component and $\text{Im}$ returns the imaginary component of the argument. Most programming languages will have the two-argument $\text{atan2}$ function. The phase difference can be calculated as $\alpha - \beta = \text{atan2}\left(\text{Im}\left(a\ \text{conj}(b)\right),\ \text{Re}\left(a\ \text{conj}(b)\right)\right)$, where $\text{conj}$ denotes taking the complex conjugate, saving you one $\text{atan2}$ evaluation and wrapping of the phase, at the cost of one complex multiplication.

There's one more trick available to calculate $a\ \text{conj}(b)$. Multiply the Goertzel output $a$ by the second input signal and lowpass filter the result. In frequency domain it looks like this:

The lowpass filtering isolates the circled spectral peak that corresponds to $a\ \text{conj}(b)$. This saves you one Goertzel filter at the cost of a lowpass filter. If you use this trick you must also take into account the phase shift due to the lowpass filter.

For fast $\text{atan2}$ see Methods of computing fixed point atan2 on FPGA.


Best method to do piano key pitch detection? For example, when playing piano, at the same time, print out the key notes by analyzing the signal of piano sound.

How to do some pre-processing to remove the noise? When calculate with FFT, noise will influence the result?

There is no easy answer to that question. Plenty of algorithms exists which are suitable to that task. Nowadays Non-negative Matrix Factorisation (NMF) is getting more and more popular in this field of research. If you have enough of resources and knowledge then you can try it. It's just a 'fancy' SVD decomposition with some 'constraints and tweaks'.

Some literature for you that might be useful in terms of choosing the approach. First one contains some vast overview of the problem as well as lot's of references:

  Benetos E., et al. - Automatic music transcription: challenges and future
  directions

  Benetos E., Dixon S. - Joint Multi-pitch Detection using Harmonic
  Envelope Estimation for Polyphonic Music
  Transcription

  Benetos E., et al. - Automatic Transcription of Pitched and
  Unpitched Sounds from Polyphonic
  Music

  Benetos E. - Automatic Transcription of Polyphonic Music Exploiting Temporal Evolution - PhD Thesis

  O'Hanlon K., Plumbley M. - Automatic Music Transcription using Row
  Weighted
  Decompositions


CMA Equaliser and FSK I'm looking into equalisation for FSK over HF channels for Digital Voice applications.  As a starting point, I am trying to understand the CMA algorithm.

I found a nice CMA Matlab example for BPSK in a previous post here on this site.

I have modified the simulation to use 2FSK but it doesn't converge to the same result.  The first plot is for BPSK modulation, the 2nd plot for FSK.  The combination of the channel and equaliser impulse response is converging to 1 for BPSK.

Any suggestions of how I can get this simulation working for FSK would be much appreciated!

Here is the source code.  the variable "tx_type" switches the tx modulation source.

N = 20000;          % # symbols
h = [1 0.45 -0.2];  % channel impulse response
h = h/norm(h);
Le = 20;            % equalizer length
mu = 1E-3;          % step size
snr = 30;           % snr in dB
M = 10;             % oversample rate

tx_type = "bpsk";   % select modulation type here "bpsk" or "fsk"

if strcmp(tx_type, "bpsk")
  s0 = round( rand(N,1) )*2 - 1;     % BPSK signal
  s0M = zeros(N*M,1);                % oversampled BPSK signal
  k = 1;
  for i=1:M:N*M
   s0M(i:i+M-1) = s0(k);
    k ++;
  end
end

if strcmp(tx_type, "fsk")
  tx_bits = round(rand(1,N));

  % continuous phase FSK modulator

  w1 = pi/4;
  w2 = pi/2;
  tx_phase = 0;
  tx = zeros(M*N,1);

  for i=1:N
    for k=1:M
      if tx_bits(i)
        tx_phase += w2;
      else
        tx_phase += w1;
      end
      tx((i-1)*M+k) = exp(j*tx_phase);
    end
  end

  s0M = real(tx);
end

s = filter(h,1,s0M);                % filtered signal

% add Gaussian noise at desired snr

n = randn(N*M,1);
vs = var(s);
vn = vs*10^(-snr/10);
n = sqrt(vn)*n;
r = s + n;          % received signal

e = zeros(N*M,1);   % error
w = zeros(Le,1);    % equalizer coefficients
w(Le)=1;            % actual filter taps are flipud(w)!

yd = zeros(N*M,1);

for i = 1:N*M-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = y^2 - 1;
    w = w - mu * e(i) * y * x;
end

np = 100;           % # sybmols to plot (last np will be plotted); np < N!

figure(1); clf;
subplot(311), plot(e.*e), title('error')
subplot(312), stem(conv(flipud(w),h)), title('equalized channel impulse response')
subplot(313);
plot(1:np,s0M(N-np+1:N),1:np,yd(N-np+1-Le+1:N-Le+1))
title('transmitted and equalized signal'), legend('transmitted','equalized'), axis([0,np,-1.5,2])

Thanks,

David

The confusion comes from the fact that what is tagged as "transmitted" isn't the real transmitted waveform but its baseband representation, which are $\left\{+1,-1\right\}$ symbols in 1 dimension. For 2-FSK, that representation would be $\left\{\left(+1,0\right), \left(0,+1\right)\right\}$ symbols in 2 dimensions.

A common convenient way to handle such 2 dimensional symbols is to represent them with complex number, mapping the coordinate $\left(+1,0\right)$ to the purely real value $+1$ and the coordinate $\left(0,+1\right)$ to the purely imaginary value $j=\sqrt{-1}$. Using that mapping, the symbol generation code for 2-FSK becomes:

if strcmp(tx_type, "fsk")
  tx_bits = round(rand(1,N));
  s0 = zeros(N,1);
  s0M = zeros(N*M,1);

  for i=1:N
    if tx_bits(i)
      s0(i) = 1i;
    else
      s0(i) = 1;
    end
  end

  k = 1;
  for i=1:M:N*M
   s0M(i:i+M-1) = s0(k);
    k ++;
  end
end

The problem then is that the CMA weight update rule from the original post is given for a purely real signal. The adjusted weight update rule for complex signals is given by:

$$
\begin{align}
  y &= \vec w_{n}^H x\\
  e_n &= \left|y\right|^2 - 1 \\
  \vec w_{n+1} &= \vec w_{n} - \mu \vec \nabla e_n \\
           &= \vec w_{n} - \mu' e_n \Re\left\{y^\star x\right\}
\end{align}
$$

This can be implemented with:

for i = 1:N*M-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = abs(y)^2 - 1;
    w = w - mu * e(i) * real(conj(y) * x);
end

Plotting the results (showing the real part and imaginary parts in separate graphs) should give you something like the following:


Intuitive understanding of scale-space extrema detection Can someone explain intuitively why local maxima and minima in the scale-space domain make for good keypoints?

I understand using LoG or DoG zero-crossing points to identify spacial variations, i.e. corners.  And I understand working at different scales to find scale-invariant features.  Both of these are pretty intuitive.

But I fail to see what it is about local extrema that makes them good features for various algorithms like SIFT.

LoG and DoG (an approximation of LoG) masks can serve as blob detectors. A blob can exist in an image at a number of locations $(x,y)$-coordinates and scales (some parameter; $t$). In some situation where scale space is divided into 3 discrete 'slices' and there are only 'small,' 'medium' and 'large' sized blobs, a 'medium' sized blob will have some response to both the 'small' and 'large' sized detectors. However, a detector whose scale is matched to the scale of the blob will have the largest response.  For this reason, we want to localize the blob by determining the maximum response in both spatial and scale coordinates.


Discrete Fourier transform - Norms of complex input signals and their transformation Given a signal $\mathbf{z} \in \mathbb{C}^n$ and its Discrete Fourier  transform $\hat{\mathbf{z} }$, does $||\mathbf{z}|| = ||\hat{\mathbf{z} }||$ hold?

The question is given to me like this with no additional details. Information about what kind of norm is also not given. Does anyone have an idea what the question might be looking for?

I may consider them as L^2 then following theorem holds.

https://en.wikipedia.org/wiki/Plancherel_theorem

https://de.wikipedia.org/wiki/Satz_von_Plancherel

https://ms.mcmaster.ca/craig/craig-pde-chap5.pdf

Also keen to learn if there is any other theorems.


Autocorrelation and Spectrum Density Wiener–Khinchin theorem states that autocorrelation function and power spectral density are a Fourier-transform pair - see Wikipidia (and lots of other resources).
That means autocorrelation should be able to be obtained by inverse Fourier transform the spectrum. The following code (run within Octave, with "pkg load signal") shows the Fourier transform of the autocorrelation DOES look like the spectrum, but the inverse Fourier transform of the spectrum does not look like the autocrrelation. What did I do wrong?
### A signal's autocorrelation and its Engergy Spectral Density are Fourier transform pairs.

###   signal
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
#plot(x);

###    autocorrelation
Rxx = xcorr(x);
#figure();
plot(Rxx); title("Rxx");

###    autocorrelation FT
RxxDftAbs = abs(fftshift(fft(Rxx)));
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
figure();
plot(freq,RxxDftAbs); title("RxxDftAbs");

###    Energy Spectral Density
xdft = abs(fftshift(fft(x)));
x_esd = xdft.^2;   # ESD is the same as autocorrelation FT. Here for visualization purpose, using absolute values.
freq = -Fs/2:Fs/length(x_esd):Fs/2-(Fs/length(x_esd));
figure();
plot(freq,x_esd); title("x esd");

### ?????????????
### is it possible to get autocorrelation from the ESD by inverse Fourier transform?
### ?????????????

###    IFT of ESD
x_esd_idft_abs = abs(ifft(fftshift(x_esd)));
figure();
plot(x_esd_idft_abs); title("x esd ift");

Thanks in advance.

Fat32's answer is correct and shows a common pitfall.
The reason that you must do this is because recall that the output of the autocorrelation of a signal of length $N$ is $2N - 1$.
You were performing the FFT with the original sample size of $N = 1000$, effectively destroying necessary information to retrieve the autocorrelation of size $2N - 1 = 1999$.