texts: How to prevent distortion after down-sampling?

texts:
I am new to audio signal processing, but in one of my machine learning projects that I am working on, I need to down-sample my audio somewhere in the middle of my neural network. Therefore, I need a down-sampling method that can apply to a tensor with the shape (batch, len). (I am using PyTorch).

Since Max pooling and Average pooling are usually used to down-sample 2D images, I am thinking about if I can try my luck with them.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample
import torch
from torch.nn import MaxPool1d, AvgPool1d, AdaptiveAvgPool1d

#-------------Generating Artifical Signal -----------------
fs = 44100
s = np.linspace(0,1,fs)
x = np.sin(2*np.pi*250*s, dtype=np.float32) + np.sin(2*np.pi*280*s, dtype=np.float32) \
    + np.sin(2*np.pi*480*s, dtype=np.float32) + np.sin(2*np.pi*1600*s, dtype=np.float32)
uprange = 800
plt.figure(figsize=(8,3))
plt.plot(x[:uprange])

#------------Comparing two resampling methods ------------
torch_downsample = AvgPool1d(2,2) # Will be used for downsampling
x_torch = torch.tensor(x)[None, None, :]
for i in range(6):
    uprange = uprange//2 # the number of samples being plotted out
    x = resample(x,len(x)//2) # downsample with fft resampling algorithm
    x_torch = torch_downsample(x_torch) #downsample with Pooling

    fig, ax = plt.subplots(1,2,figsize=(16,3)) # plot out the resampled waveforms
    fig.suptitle('Resampled {} time(s)'.format(i+1), fontsize=16)
    ax[0].plot(x[:uprange], '-b')
    ax[0].plot(x_torch.view(-1).numpy()[:uprange], '--r')

    z = np.fft.fft(x) # Analyzing the spectrum
    ax[1].plot(abs(z)[:len(x)//2])

    z = np.fft.fft(x_torch[0,0,:].numpy())
    ax[1].plot(abs(z)[:len(x)//2], '--r')
    print("len1 = {}, len2 = {}".format(x.shape, x_torch.shape))
    ax[0].legend(['fft_res', 'AvgPool'])
    ax[1].legend(['fft_res', 'AvgPool'])
    plt.show()

My results

As you can see for the first 4 resamplings, the frequency components are still correct. But starting from the 5th resampling, there's a small defect at around 200Hz. The 6th resampling results in more defects.

My Questions

Is there any method to prevent these defects?
Also, am I doing downsampling correctly with Average pooling? I have
no knowledge in audio signal processing, but I am trying hard to use
my machine learning knowledge here.

Preferably, the solution can be implemented in PyTorch. (Don't tell me to use torchaudio, it messed up my whole PyTorch installation in an attempt to install it. I don't want to go through the trouble to trouble-shoot and reinstall everything)


texts: Rate offset and jitter concealment for streaming audio

texts: Given streaming audio samples between two separated devices with separate non-synchronized but almost identical sample rates (say both 44100+-2ppm), what are some methods to best conceal any difference and/or jitter in rates from causing audible artifacts in the audio output, and with minimal required added latency for the concealment algorithm?


texts: How to calculate the displacement from the measured acceleration data?

texts:
I have an acceleration data which is collected from a vehicle. I do also have a displacement measurement.

I would like to calcuate the displacement vector from the measured acceleration data, but no success yet.

Here is the data (the previous link was not working I updated the link):
https://ufile.io/zyw3r

What I try:

filename = ('acc_time.csv');
accdata = csvread('acc_time.csv');

[mydata, myheader] = xlsread(filename);
for i = 1:length(myheader)
      % compose a command to assign each column to a variable with the same
      % name as the header
      commandExec = [myheader{i}, ' = ', 'mydata(:,', num2str(i) , ');'];
      % execute the composed command to actually create the variable
      evalin('base', commandExec );
end

vel_vector = cumtrapz(time,acc);
disp_vector= cumtrapz(time,vel_vector);
plot(time, disp_vector)

and I get:

but the displacement vector should be squiggy not like linear incremental.

Where am I wrong?


texts: Low-rank approximation of the STFT matrix : good audio compression technique?

texts:
I recently heard about this interesting math tool:

With a matrix $A$, compute the SVD ($ A = U \Sigma ^tV$ factorization, where $\Sigma$ is a diagonal matrix containing the "singular values")
Keep the $k$ biggest singular values only and discard the lowest ones (set them to $0$); this gives $\tilde{\Sigma}$.
Then $ \tilde{A} = U \tilde{\Sigma} ^tV$ will be the best rank-$k$ approximation of A.

Step #2 allows to save a lot of storage size, and thus this can be thought as a compression technique. See this example on image processing.

Then I wanted to try this on STFT. I did :

Audio sound ---> STFT ---> low-rank-approximation of the STFT matrix thanks to SVD ---> iSTFT

The good thing is that it helps to save storage size, compression!
The bad thing is that it only resulted in loosing the high frequencies :

Original STFT (listen to the audio file):

After low-rank-approximation STFT (high frequencies more or less lost) (listen to the audio file):

Is there a way to turn this "low-rank approximation of STFT" into something better? And if not, why?


texts: Is the STFT a good solution to guitar monophonic pitch detection?

texts:
Let’s assume that the audio input will be:

Guitar audio
Monophonic audio

Is the STFT overkill? Is it enough to produce a time-frequency plot and from there find the notes played at each instance? What problems could arise?


texts: best robust feature detectors for detecting similarities among archeological signs

texts:
I would like to implement edge detectors for creating a feature vector for archeological signs.
Which are best features, algorithms to retrieve? That could be loaded in a Support vector machine for classification.

If there are in the image only signs with a uniform background it's better converting an image in a binary image and calculate convexHull, solidity, eccentricity, elongation of the archeological signs.
Then, what i do? I insert these data in a vector and put every vector of every image in a Support Vector Machine for training?

How can i measure "distance", similarities among images?


texts:
I am using FFT to detect certain frequencies (DTMF). I can clearly see the frequencies in the spectrum. Is there a preferred / optimal way to decide if a certain frequency is present or not by looking at the spectrum?

I can think at least of:

Bin absolute value (ie. bin amplitude)
Difference between maximum bin and second highest
Ratio between maximum bin and second highest
Ratio between maximum bin and sum of all others


texts: DFT and multiplication/convolution equivalence

texts:
Is there a simple or potentially intuitive explanation for, with the DFT, vector multiplication in one domain being equivalent to circular convolution of the transforms of the vectors in the other domain?

Since a DFT is just multiplication by a (special) square matrix, what about this matrix and matrix multiply allows the above duality?


texts:
If we have data we need to transmit which is $i$, so we should multiply is using XOR function with spreading code $y$.

How many bits (Maximum and Minimum) should $i$ and $y$ can have ?


texts:
I am implementing phase correlation algorithm to determine shift between two images.

It generally works, but I am not sure how to interpret the resulting shift.

Pseudocode:

ft1 = FT(input1)
ft2 = FT(input2)

for each pixel p
  product = ft1[p] * conj(ft2[p])
  cross_power[p] = product / abs(product)

cross_correlation = IFT(cross_power)

where FT and IFT are Discrete Fourier Transform and its inverse, respectively.

Here are input images and resulting phase correlation taken from Wikipedia:

The peak is located at (20, 21) which corresponds to shifting of the second image.

I have obtained similar result, but it comes out flipped:

When I flip the input images the peak appears at the top left corner of the phase correlation image but how to interpret this as a negative shift (-20, -21)? Based on image quadrants?

I am using Fourier transform routine from Math.NET Numerics package. It provides only 1D FT, so I am doing transforms of rows and then columns to obtain 2D FT.

How to extract shift from the peak position allowing for negative shifts as well?

Why is the peak located differently in my result?


texts: Recover Filter Coefficients from Filtered Noise

texts:
I have a digital signal which may be represented as a pulse noise source filtered with an FIR (finite impulse response) filter. Suppose that the noise consists of discrete pulses (nonzero samples where most samples are zero), and there is an equal probability of a pulse at any sample (so for example the number of pulses in an interval follows the Poisson distribution). Further suppose the strength of each pulse is equal.

Can the coefficients of the FIR filter be recovered from the filtered signal, and how?

Now assume the noise is impulsive, but may be somewhat correlated with itself (it doesn't have perfect Poisson statistics). Also assume that the strength of each pulse may vary. Can the coefficients be recovered approximately if the exact distribution of noise is not known?

Noise source $N[i]$ where $i=0..n, N[i] = 1$ with some unknown but low probability $p$, 0 otherwise; each $N[i]$ is independent of other $N[i]$

Filter $F[k]$ where $k=0..m$, unknown, $m$ is much less than $n$

Signal $S[i]$ where $i=0..n$, $S = N * F$, where * is convolution

Given $S$, estimate $F$. $N$ and $p$ are also unknown.

To give this some practical background, the average period between noise pulses can be 20-40 samples, the FIR filter can have a few hundred nonzero coefficients (so the filtered signals resulting from each pulse overlap significantly), and the total signal is a few thousand samples.

EDIT: Some real-world data here: http://pastebin.com/N3LZU6vm  I think what you see at sample 780-810 is pretty close to the impulse response to a single impulse.  This is generated by noise fed into a complex resonator, whether the noise is impulsive as described above is not actually known but is physically plausible.

Some synthetic data here: http://pastebin.com/qTYLP5eA  This is generated from the following python code:

noise = numpy.where( numpy.random.rand(5000) < 1/20.0, numpy.ones(5000), numpy.zeros(5000))
filter = numpy.sin( numpy.arange(200) * 2 * math.pi / 11 ) * numpy.exp( - numpy.arange(200) / 30.0 )
signal = numpy.convolve( noise, filter )
for i in range(len(signal)): print int(signal[i] * (2**15)/ 3.0)

This creates a filter wavelet from a sinewave with period 11 samples and an exponentially decaying envelope, and noise where each sample has an independent 5% probability of being an impulse.

EDIT:  I think what I am looking for is very similar to what is involved in building an autoregressive model (AR).  There are slight differences in the formulation, for example AR is based on an IIR with lagged signal plus a noise term rather than FIR with lagged noise, and the noise is assumed to be Gaussian rather than impulses with Poisson distribution.  However I think AR proves that estimating a filter from filtered noise is possible, and in fact AR may just work for my problem.


texts: Defining the convolution for finite-length signals

texts:
How can the linear convolution be defined rigorously for two discrete signals $x = [x(0), …, x(N_x-1)]$ and $h = [h(0), …, h(N_h-1)]$ of different finite lengths $N_x$ and $N_h$ respectively?
Let's assume for simplicity that $N_x = 3$ and $N_h = 2$, i.e. $x = [x(0), x(1), x(2)]$ and $h = [h(0), h(1)]$. Does the sum go from $0$ to $N_x - 1$ or to $N_h - 1$? However,
$$\sum_{i = 0}^{N_x-1} x(i) h(n - i) \neq \sum_{i = 0}^{N_h-1} h(i) x(n - i)$$
i.e. the commutativity cannot hold, and also both sums will have "illegal" terms containing $h(-2)$, $h(-1)$, $h(2)$, $h(3)$, $h(4)$ and $x(-1)$, $x(3)$, $x(4)$ (which I guess can be set to zero, but that somehow seems not right since that redefines my initial signals).


texts: Is it possible to infer the characteristic curve of a recording device only from (many) recordings?

texts:
... with a certain degree of confidence / corridor of accuracy.

What do I mean?
Case in point: There is the huge archive of digitized Edison wax cylinder recordings.

They have annotations, usually year of recording, sometimes producer of recording.
Assuming that, with a good likelihood, grouping the intersection of e.g. one year of recording and a geographic area or even company yields sets of recordings very likely made with the same model / iteration of the technology within a set, given how slow-paced development was back then:

Could one infer the rough frequency response curve of an imaginary "average recording device" of those times from such a narrowed-down set of recordings?

My not exactly DSP expert naive train of thought went like this:

if I have a great enough number of as diverse as possible audio recordings, their sources (i.e. before entering the brass horn to be recorded) should have an average spectrum roughly resembling white noise. This probably means using broadly set up orchestral recordings, not so much "solo singer with background accompaniment", which is weighted too much to vocals in the spectrum.
with a good likelihood of being made on the same model of recording apparatus, an "average spectrum" calculated e.g. with FFTs over the whole set of recordings as described above, should yield the rough response curve of the imaginary recording device

Now, if premise 1. is already wrong, then I'm out of luck.
And 2. could be wrong if the spread of characteristics between individual exemplars of devices of a series is too great - although I guess that's a bit of a question with regards to "what do I need".
Anyway, this was just to show I made some effort thinking about this. If it doesn't work the naive way I am imagining:

Is there a way to do that, only from recordings, without being in possession of an actual device of that type (for measurement)?

The amount of distortion by e.g. saturation within the mechanical recording process can probably not be inferred just from recordings?


texts: User-friendly DSP for high-frequency application?

texts:
I've been doing some work with creating visual output on oscilloscopes using audio; here are some examples:
http://www.youtube.com/watch?v=UJ8-vySqEAI

This is using Max/MSP to read from a buffer of 3d geometry as audio output to an oscilloscope screen. I'm interested in doing transformations on full frames of video as well, like a digital/analog hybrid version of the old Scanimate machines from the 70s:
http://www.youtube.com/watch?v=wTO3CKM2Yy0)

But at 44100 samples/second, I don't have the bandwidth to keep a frame of video on the screen at a high enough scan rate without the signal degrading.

What I'd like to do is to use something like MiniDSP to:
1. Load a frame of pixels (monochrome) into a buffer
2. Read out the buffer to modulate the beam's brightness and move it along a raster scan path at high speed (around 60fps ideally)
3. Add LFOs to modulate the scan path to warp the image

This requires a high sampling rate for output, probably higher than normal audio rates for NTSC video, and I know MiniDSP is mainly for audio applications. I'm wondering if you can recommend a solution, ideally one that has friendly features like USB connectivity and C++ compiler, or am I dreaming? Can MiniDSP be overclocked to work at higher rates for output?

Thanks very much for any help.


texts: Is there a distinction between Tapped-Delay-Line Filters and the more general FIR (Finite Impulse Response) filters. From what I understand the Tapped-Delay-Line is a type of FIR filter since it's impulse response settles to zero in a finite duration...but I'm having trouble determining what might make these two things separate in the world of DSP.


texts: I have two .wav audio(stereo) files to be compared. The files may vary in frequency or amplitude or phase shift. How to compute the difference between the wave files in Scilab?  The expected output is that if frequency and amplitude differ between the audio files, the output must return that they are dissimilar.   Phase changes need not be compared.


texts: Is the whole data needed for the fast fourier transform with an ongoing signal?

texts:
I sample an signal with around 840k values a second, to check the spectrum around <400kHz. But i can't save the whole data as it is just to much for my microcontroller. I know that for the butterfly-calculation you need values from begining and the end.
So is the whole set of data realy needed to calculate the fourier transform or is there some other way? And how long does it take to be able to process the newest samples of data.


texts: Why are state-space filters good for time-varying systems?

texts:
"The state-space filter representation is especially powerful for time-varying linear systems."

Why is this the case?
What is the explanation behind the claim?


texts: MATLAB's $\tt filter()$ function for complex valued data

texts:
I want to use the filter() function to implement an FIR filter having complex data values as the input to it. For real, I used this function as filter(A,1,B).

Question :

Would there be a difference for the corresponding implementation when A and B are in complex domain?


texts: Finding the linear prediction coefficients for a sampled sinusoid

texts:
I am unable to work my way through this problem -

A sinusoid with random phase is given as () = 2 sin(0.25π + ). If this signal is sampled at 10
KHz, estimate the coefficients of a second order LP model. Find the magnitude and phase values of the
second order LP error filter at  = 0.5. Estimate the bias (error) in the peak frequency (frequency at the
maximum magnitude response) of ()

I am thinking of proceeding with Yule-Walker equations but unable to solve this.

Can you please provide a way to solve this.


texts: Sensitiviy and Resolution in Windowing Functions

texts:
When I read the  article Windowing Functions, there are two terms that are does not understanding them well, sensitiviy and resolution.

I looked for this two terms on the web, And I conclude that sensitivity are how much you can distinguish the components of the signal, and resolution how many components you can have for a signal. Is that correct?

If my above conclusion was correct I see that the more resolution the more components I have, yields less distinguishing of that components, which means the more resolution the less sensitivity, does that also is correct? Also, does the overlpping windows is a solution for the little resolution that some windows has, for example hanning window?


texts: For linear IIR digital filter, what happens for negative frequencies?

texts: By negative frequency, I refer to Fourier transform. Often, the frequency response of a digital filter is only displayed for positive frequencies. For a linear IIR digital filter, what happens for negative frequencies? Are frequency response for negative frequencies mirror images of what happens for positive frequencies?


texts: Doing analysis using IQ demodulation in Python

texts:
I'm trying to doing an analysis using IQ demodulation. I'm doing this by taking a noisy signal which has a desired signal in it. The desired signal has a frequency that is almost stable, which I'm simulating to be 5 hz + a small amount of random fluctuations in instantaneous frequency. I'm then mixing it with a cosine and sine function of 5 hz, then applying a lowpass filter. Then I'm taking the square root of these two filtered signals squared to retrieve the amplitude, which I simulate to be decaying exponentially. I'm not sure where my problem is, but the amplitude I extract is not stable. I'm not sure if it is a misunderstanding of something fundamentally or a code problem, but I don't think it is a code problem. I am self-teaching some DSP, and would not consider myself an expert at all. Here is the code:

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt

pi = np.pi

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

t = np.linspace(0,20,20001)
sample_freq = 1/(t[1] - t[0])

f1 = [3 + .2*np.random.rand() for x in t]
Amp1 = np.exp(-2*t)
cos1 = Amp1*np.array([np.cos(2*pi*y*x) for y,x in zip(f1,t)])

f2 = [5 + 0.1*np.random.rand() for x in t]
Amp2 = 5*np.exp(-t)
cos2 = Amp2*np.array([np.cos(2*pi*y*x) for y,x in zip(f2,t)])

f3 = [8 + .2*np.random.rand() for x in t]
Amp3 = np.exp(-t/2)
cos3 = Amp3*np.array([np.cos(2*pi*y*x) for y,x in zip(f3,t)])

f4 = [.2 + .05*np.random.rand() for x in t]
Amp4 = 0.1
cos4 = Amp4*np.array([np.cos(2*pi*y*x) for y,x in zip(f4,t)])

noisy = (cos1+cos2+cos3+cos4)

cos_mix = np.cos(5*t)
sin_mix = np.sin(5*t)

X = noisy * cos_mix

Y = noisy * sin_mix

cutoff = .5

X_low = butter_lowpass_filter(X, cutoff, sample_freq)
Y_low = butter_lowpass_filter(Y, cutoff, sample_freq)

Amplitude = 2* np.sqrt(X_low**2 + Y_low**2)

plt.plot(t, Amplitude)

Which plots an image like this:


texts:
I am fairly new to GNURadio. I am trying to generate a 2FSK with minimum frequency shift. Instead of multiplying with a 8MHz complex carrier, I directly use VCO to generate two frequencies. One at 7.5MHz and the other at 8.5MHz. At the output, QT Gui Frequency Sink shows these two frequencies. I have attached my flow graph. Since I need a continuous phase, I use the VCO.
I have two questions.

When I use 25M sampling rate, I get different amplitudes in the signal. I only want the frequency difference but amplitudes are also changed. For 18M sampling rate, I do not see two frequencies with complex output, hence I use a higher sampling rate. The amplitude change does not exists for small sampling rates (like 48k) and vco sensitivities. I use the parameters according to https://wiki.gnuradio.org/index.php/VCO Could someone please tell me what I'm doing wrong here?

VCO complex has an imaginary component where as VCO (float) does not. When I take the VCO float output, is it equivalent to the 'real' part of modulated transmitting signal? (Is it equivalent to generating IQ signals and multiplying with a complex carrier and filtering the real component?)


texts: Sliding DFT: I think my results are accurate, but I can't buy a swap-swap IDFT

texts:
I've got a sliding DFT implementation that appears to be working (judging from an output plot). I would like to be able to invert this implementation using the standard tricks of swapping the real and imaginary components or conjugating my data. However, when I attempt this my output looks nothing like the original. I also end up with a strong imaginary component where that was not expected. (I'm not entirely sure on the proper range of output for each pass through this filter.) Can you see anything obviously wrong with this code? Thanks for your time.

    [TestMethod]
    public void TestInverseDft()
    {
        var rand = new Random(42);
        var dataSize = 50000;
        var data = new List<Complex>(dataSize);
        for(int i = 0; i < dataSize; i++)
            data.Add(new Complex(rand.NextDouble() * 2.0 - 1.0, 0.0));

        var clone = new List<Complex>(data);

        var n = SlidingDFT(data, 44100, 12, 55, 13);

        Swap(data);
        SlidingDFT(data, 44100, 12, 55, 13);
        Swap(data);

        for(int i = 15000; i < dataSize - 15000; i++)
            Assert.AreEqual(clone[i].Real, data[i].Real / n, 0.00001);
    }

    public static int SlidingDFT(IList<Complex> data, int sampleRate, int binsPerOctave, double minFrequency, int stepsFromMin)
    {
        var Q = 1.0 / (Math.Pow(2.0, 1.0 / binsPerOctave) - 1.0); // denominator * Fk gives size of bin
        var k = (int)(Q + 0.5);

        // frequency at k:
        var Fk = minFrequency * Math.Pow(2.0, (double)stepsFromMin / binsPerOctave);
        var Nk = (int)(Q * sampleRate / Fk + 0.5);
        Complex exp = Complex.Exp(new Complex(0.0, 2.0 * Math.PI * k / Nk));
        var current = data[0];
        for(int x = 0; x < data.Count - 1; x++)
        {
            var lookahead = x + Nk >= data.Count ? new Complex() : data[x + Nk];
            var fx = (lookahead - current) + data[x]; //  * 0.9999999 if unstable?
            current = data[x + 1];
            data[x + 1] = fx * exp;
        }
        return Nk;
    }

    public static void Conjugate(IList<Complex> data)
    {
        for(int i = 0; i < data.Count; i++)
            data[i] = Complex.Conjugate(data[i]);
    }

    public static void Swap(IList<Complex> data)
    {
        for(int i = 0; i < data.Count; i++)
            data[i] = new Complex(data[i].Imaginary, data[i].Real);
    }


texts: Isolate High Variance vs Low Variance Sections of Signals

texts:
I have the following one-dimensional signals:

The goal is to identify the sections with smaller variances compared to the rest (For Python, that is to get their index value in the array). In the first signal, there are two such sections (one at the left end and one at the middle). In the second signal, there is one such section at the middle.

What I've tried so far

I've tried computing the local stdev in each signals with window size 25 and produced the following plots:

From here on, I can use threshold to select sections with small stdev. However, one can see that both have different y-axis, and need different threshold values to isolate the low variance sections. Furthermore, applying a rolling window is notoriously computational expensive (I have to apply this on ~5,000,000 similar one-dimensional signals). So the question is, what is a quick method to identify sections that can be applied to both signals?

EDIT

I've put two sample signals at here and here.


texts: Testing user's pronunciation of Arabic against correct pronunciation

texts:
I'd like to be able to compare two audio streams: one containing the correct pronunciation of Arabic text (from the Quran), and the other containing one from a user learning Arabic. Here is a sample of what Arabic sounds like in case it helps: https://youtu.be/ZYaZ6Odbx_Y
Questions I'd like to be able to answer:

Did the user pronounce everything correctly? Or close enough.
Which words did the user pronounce correctly? This will require splitting the audio streams into words.
Which letters did the user pronunce correctly? This will require further breaking down the word into letters, which I'd imagine should be similar to #2. Perhaps a little trickier to identify some letters, but I have the textual source to help and I might be able to do a large preprocessing of the data to learn the profile of each letter.
Going deeper into the intricate rules of Arabic. I have the domain expertise in Arabic that would help here, and it's also something I'd rather postpone solving to later.

While my general algorithm knowledge is strong, my knowledge of signal processing is fairly limited.

Do you have any suggestions of how to solve this problem? Do any existing libraries implement such analysis?
Where can I start to read up on this topic?
I hope this allowed here, but if not tell me and I'll remove: I'm looking to hire someone to build a prototype of this. If successful, I hope to develop this into a shippable Android app. I have a design doc if anyone is interested.


texts: Practical insightful time domain functions

texts:
This question might be a little bit different, but i am having hard time to find some sort of a list/book/website/advice, with very practical time domain functions that provides useful insgihts for a series of numbers in time domain (frequencies of less than 1Hz)

Examples - FFT to find periods in a noisy signal, Autocorrelation  for finding an interesting correlations between 2 curves.

Very curious to know about more, maybe more functionality about prediction and information that is useful and critical in real life.


texts: What is the difference between feature detectors and feature descriptors? Which among these are detectors and which are descriptors: Harris, SURF, Min Eigen, FAST, SIFT, BRISK


texts:
We have this signal:

$$x(t)= 2 + 2\cos(2\pi f_0 t) + \frac{2}{T_0} \operatorname{sinc}\left(\frac{2 t}{T_0}\right)e^{j2 \pi 4/ T_0} + \operatorname{sinc}\left(\frac{2 t}{T_0}\right)e^{-j2 \pi 4/ T_0}.$$

I have calculate Fourier trasformed and have found:
$$
X(f)=\delta(f)+ \delta(f-f_0) + \delta(f+f_0) + \operatorname{rect}\left(\frac{f-\frac 4{T_0}}{\frac 2{T_0}}\right) + \operatorname{rect}\left(\frac{f+\frac 4{T_0}}{\frac 2{T_0}})\right)
$$

I must find the average power. If use Parseval's have:
\begin{align}
P_x&=\frac{1}{T_0} \int_{-\infty}^{\infty}\lvert X(f)\rvert^2df \\
&= \frac 1{T_0}\left(\int\lvert\delta(f)+\delta(f+f_0)+\delta(f−f_0)+\operatorname{rect}+\operatorname{rect}\rvert^2df\right).
\end{align}
How can I solve this?


texts: Difference between jpg compression and mp3 compression

texts:
My understanding of jpg compression is that, if done repeatedly, any image will eventually be reduced to a single flat color, due to the way it deals with neighboring pixels. However, years ago I was told that mp3 does not have this problem, a 192kbps song will stay at the same audio quality with little file size variation no matter how many times it is compressed at the same setting.

Why is that? What's the differentiating factor?

Thanks


texts: Obtaining the FFT of a data set collected through an accelerometer

texts:
I'm using MATLAB to analyse a very large sample of acceleration data collected using the MPU9150 IC. My aim is to plot the FFT for this large sample of data.  The daya is of three axes; $x,y,z$ and I have successfully imported the data into matlab as vectors.

As I am new to MATLAB I have taken over month to plot the FFT graph but I believe its highly dependent on the smoothing factor. I wish to first illustrate my approach and then showcase my full code for your kind reference. Sampling frequency of of the accelerometer is 10 Hz.

Approach in exact order.

Import the raw acceleration vector with $x,y,z$ axes Determine
butterworth Filter parameters
Apply butter-worth filter using filtfilt()
use smooth() function with a smooth factor. ( NOTE: smoothing factor
effects FFT shape significantly)
obtain FFT using abs(fft(acceleration_DC_shifted_Filtered_Vector))
smooth the FFT line using smooth function (this smoothing has no
effect on FFT there fore can be skipped if required.)
plot FFT

My code is below,

%importing data from workspace.
     accelerationd2 = Xaccelday2;
% Setting basic Parameters. Number of samples should be same in all day files.         number_of_samples = length (accelerationd2); sampling_frequency = 10;

%Setting filtering paramters.
          cut_off_frequency = 5;
          wn = cut_off_frequency/sampling_frequency;

           order_of_butterworth_filter = 1;
              [b1,a1]= butter(order_of_butterworth_filter,wn,'low');

%Filtering data from all ten days.
                  accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

%DC shift removal`

              accelerationd2_smoothned= smooth(accelerationd2_filtered,36.5);%acceleration 1filtered was passed into all smoothfunctions
          accelerationd2_DC_shifted = accelerationd2_filtered-accelerationd2_smoothned;

%Computing FFT
          accelerationd2_FFT = smooth(abs(fft(accelerationd2_DC_shifted)),8000); %1000 is ballpark
          bin_vals = (0 : number_of_samples-1);
          fax_Hz = bin_vals*sampling_frequency/number_of_samples;
          N_2 = ceil(number_of_samples/2);

%Detecting peaks
          figure;
          plot((fax_Hz(1:N_2)),(accelerationd2_FFT(1:N_2)),'r');
          hold on
          legend('Day-02');

%[maxtabDay1, mintabDay1] = peakdet(accelerationd1_FFT(1:N_2), 15, fax_Hz(1:N_2));
          [maxtabDay2, mintabDay2] = peakdet(accelerationd2_FFT(1:N_2), 15, fax_Hz(1:N_2));

%This is where peaks are dressed with markers. Markers must have the same color of the line. Markers can be different symbols. Types of markers that can be used : ^ , *
           plot(maxtabDay2(:,1), maxtabDay2(:,2), 'r*');
          xlabel('Frequency (Hz)');
          ylabel('Magnitude (dB)');
          title('Single-sided Magnitude Spectrum (Hz)');
          %hold on
          hold off
%axis tight

My problem is that the FFT changes highly as the smoothing factor in

 accelerationd2_smoothned= smooth(accelerationd2_filtered,36.5)

is changed. Any help would be so much appreciated.

Is it my thinking that there is only one correct FFT and the FFT should not change depending on the process I obtain it from?

Another important reason to note is that all samples have a DC shift. Which is why the smooth function has been used. If it is not used, the time domain is not brought about $y=0$-axis before inserting to the FFT function.

I would appreciate if a better procedure can be explained that can be used to obtain the FFT through MATLAB.

P.S.
It would be nice to have a method I can obtain the FFT for my data set. I collected this data about the motion of an particular animal. And I need to generate the FFT for its movement. I have many records but not a proper method. My data goes like,
$X_{\rm acceleration}, Y_{\rm acceleration}, Z_{\rm acceleration}$.
But there can be DC shifts seen on the time domain because the animal because the collar which has the logger was not perfectly in the middle. Therefore I must remove these DC components caused by gravity.


texts: Why do we quantize high frequencies coarsely, shouldn't it be the opposite?

texts: Also, in Huffman Coding, why do we use fewer bits for high frequency data. I don't get the intuition. Shouldn't it be the opposite or does this have to do something with the psychoacoustic model of the human hearing which says that our cochlea is sensitive to lower frequencies more than higher, hence we use up higher bits for lower frequencies?


texts:
I am working on a simulation of QPSK data using the code below (Matlab):

j = sqrt(-1);
N = 10^6; % number of symbols
Es_N0_dB = 12; % ratio of signal energy to noise energy in dB
QP = (2*(rand(1,N)>0.5)-1)+j*(2*(rand(1,N)>0.5)-1); % actual QPSK data
s = (1/sqrt(2))*QP; % normalization of energy to 1
n = 1/sqrt(2)*(randn(1,N)+j*randn(1,N)); % white Gaussian noise, 0 dB variance
y = s + 10^(-Es_N0_dB/20)*n; % recovered noisy data
plot(y, 'b.')

This is my resulting constellation diagram:

To identify the errors (and calculating the bit-error-rate), I have divided the diagram into 4 quadrants and used the following to find out which elements have changed sign (i.e. crossed the decision boundaries):

sr=real(s); si=imag(s);
yr=real(y); yi=imag(y);

for i = 1:N
    if sr(i)*yr(i)<0
        errorI(i)=yr(i);
        if si(i)*yi(i)<0
            errorQ(i)=yi(i);
        end
    end
end

err = errorI+j*errorQ;
BER = numel(err)/numel(y);

However, this does not return any data for the in-phase and quadrature errors. What is the problem with this approach?

Any explanation would be greatly appreciated.


texts: Using DFT to determine control system bandwidth from user input

texts:
I am currently working on a race car project where steering and braking is to be done by a fully electronic x-by-wire system. To roughly determine the minimum bandwidth of the controllers, I took user + car data from an existing race car. I though I would do an FFT of the input signals to see what the 'user' bandwidth is, and let the controller sample time be at least 10 times this value.

As an example, this is what I did with the steering input:

I took the sampled time-series of the steering wheel angle, s(t), for 2 laps, and removed the mean value from the series.
I determined |s(t)|^2 to obtain a power-like quantity from the user input
An FFT was performed on the signal from step 2, and the bandwidth of this signal was determined. At first I was going to determine the half power bandwidth, but I quickly noticed the spectral components above 2 Hz were negligible so I took this as the 'user signal bandwidth'
From step 3, the min sample frequency for the steering input is 20 Hz as rule of thumb.

My questions are:

was I correct in thinking that an FFT of the input signal represents
how much input 'energy' is present at each frequency in the FFT?

does step 2 makes sense or should I skip it?
Am I correct in thinking the system will steer nicely if I sample the input signal at at least the frequency from step 4?
My grad prof asked me 'so an FFT of the steering input, what am I seeing here? What does this graph actually mean? Is the driver
actually turning the wheel at 1/3 Hz, what have you plotted?' and I
had no idea how to clearly answer this question. Any thoughts?

I am new to this forum so hope my questions are clear enough for you guys.

Thanks in advance for helping me out!


texts: extracting ultra accurate phase and frequency with FFT with real world data

texts:
I was doing some experiments in Matlab, and write this code to experiment with it.

I want to know if it's possible to extract accurate within -+0.01° phase info from samples of an ADC, the matlab FFT function says that we can.

But it's weird somehow:   It would get me very accurate results, is it doable in a real-world scenario?

I assumed a 64 Sample per Cycle in my Wave under test.

Here is the matlab code

t=0:2*pi/64:2*pi-2*pi/64;
R=sin(t+pi*178.03/180);
FFT=fft(R);
angle(FFT(2))*180/pi

After this code I start to do some measurements on my real hardware, and I have increased my sample rate to 128 samples per cycle.

I have made some measurements with my system, which is a Cortex M4@200MHz device and an external ADC, I have a precision source injecting voltages and currents to my ADC which is connected to  CT and PT's, the first 4 channels are CT’s and the last 4 channels are PT's.

I have injected some sine waves to the system from 40Hz to 400Hz and in the 50Hz zone I have incremented the frequency by 0.1Hz, also I have tested the system with various phase angles with the channels, the sampling frequency of my ADC is 6400 sample/sec and my ADC is a true simultaneous one.

Here is the results, when I extract the phases respect to the V0 channel and in the 50Hz region, with different phases the extracted phase is within good ranges, but as the frequency changes the error would increase.

I have only 128 samples in a full 50Hz cycle and I cannot increase my sample rate higher due to my system's overhead. I have to extract all the phase and frequency of the waves from this 128 samples, roughly 5ms, aside of many other things to calculate.

Feel free to play with these data to see if you could suggest a good enough solution to extract phase and frequency info with under 0.1 degree for phases and 0.01Hz for frequency.

And all the extracted waves are ready to be injected into matlab.

Here you can get the wave files

https://ufile.io/5b7y5


texts: Generating a noisy sine wave in Python, efficiently

texts:
(Questions after code block.)

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

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

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

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

def writeWav(w):
    # audiophilic magic

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

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

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


texts: Direct form 1 filter equation for a biquad filter

texts:
In general, we all know that the direct form 1 equation for a biquad filter is given by
$$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$
But when it comes to practical applications like for audio processing applications, I have seen many people using
$$y(n) = b_0x(n) + 2b_1x(n-1) + b_2x(n-2) - 2a_1y(n-1) - a_2y(n-2)$$
as the filter equation in their codes.
So, I just want to know the reason why we need to multiply $a_1$ and $b_1$ coefficients by 2?


texts:
I am trying to find the autocorrelation function on the intensity of a pixel over time stored in an array "Intensity". Usually one would expect the highest point to be at t = 0 and then decay. I am getting a random mess.

I checked my code with multiple posts on here regarding autocorrelation such as Autocorrelation: numpy versus FFT I am working in python

correlation = np.correlate(Intensity, Intensity, mode='full')/Intensity.size
correlation = correlation[correlation.size//2:]
plt.plot(range(0,len(correlation)),correlation)
plt.show()

You might say well maybe your data looks very similar all around, but then either way why would the highest value not be at t = 0? its not even close. I feel like I must be doing something wrong


texts:
It seems that you can decompose it as such:

$f(n) = a^n u(n)  +  a^{-n} u(-n-1)$

But I already have issue here,

is it basically saying that $ u(n) + u(-n-1) = 1$?

this is the plot of u(n) and u(-n-1):

if you sum them, wouldn't you just basically combine them? So you get 1 from $- \infty$ to -1, and from 0 to $\infty$

between n = -1 and 0 it's 0

no?

I know that the answer to z-transform $2^k$ is $\frac{z}{z-2}$ but how do you arrive at this?

general formula of z-transform is:

$F(z) = \sum_{k=\infty}^{0} f[k] z^{-k} $

applying formula:

$F(z) = \sum_{k=\infty}^{0} 2^k z^{-k} $

it would look like this, no?


texts: suppose I have this OFDM signal the pink part is cyclic prefix. I should sample at point b. Because of lack of Timing I started to sample at point a, but I get wrong detection. How can I compensate this problem. I appreciate any help.


texts: Accuracy of Goertzel algorithm for measure harmonic

texts:
If I use Goertzel algorithm for calculate the magnitude of current or voltage, use for measure current harmonic or voltage harmonic or fundamental component. I not use for detect harmonic. Do you know the accuracy for measure harmonics is good or not?,How much about relative error?

Another one question is The Accuracy relative with number of data sample?
Ex. If I have measure voltage harmonic,
 I measure waveform of current by use sample rate = 4000 Hz, target frequency for measure fundamental harmonic is 50Hz, so if I use number of sample (bin width) for Goertzel algorithm is 80 sample this is minimum, If I use more then 80 samples like 4000 sample the accuracy will more accurate right? how much better?


texts: Convolution of signals sampled on a logarithmic grid

texts:
Is there a practical accelerated algorithm or a theoretical discrete (Fourier) transform based method to convolve discrete-time signals sampled on a logarithmic grid? What I mean is representing a discrete-time signal by the sequence:

$$[1, 1, 1],$$

as a shorthand notation for samples:

$$\big[\big(\log(1), 1\big), \big(\log(2), 1\big), \big(\log(3), 1\big)\big],$$

with autoconvolution:

$$[1, 1, 1]*[1, 1, 1]$$

$$=\big[\big(\log(1)+\log(1), 1\big), \big(\log(1)+\log(2), 1\big), \big(\log(1)+\log(3), 1\big)\big]+\\
\big[\big(\log(2)+\log(1), 1\big), \big(\log(2)+\log(2), 1\big), \big(\log(2)+\log(3), 1\big)\big]+\\
\big[\big(\log(3)+\log(1), 1\big), \big(\log(3)+\log(2), 1\big), \big(\log(3)+\log(3), 1\big)\big]$$

$$\begin{align}=&\big[\big(\log(1\times1), 1\big), \big(\log(1\times2), 1\big), \big(\log(1\times3), 1\big)\big]\\
+&\big[\big(\log(2\times1), 1\big), \big(\log(2\times2), 1\big), \big(\log(2\times3), 1\big)\big]\\
+&\big[\big(\log(3\times1), 1\big), \big(\log(3\times2), 1\big), \big(\log(3\times3), 1\big)\big]\end{align}$$

$$\begin{align}=&\big[\big(\log(1), 1\big), \big(\log(2), 1\big), \big(\log(3), 1\big)\big]\\
+&\big[\big(\log(2), 1\big), \big(\log(4), 1\big), \big(\log(6), 1\big)\big]\\
+&\big[\big(\log(3), 1\big), \big(\log(6), 1\big), \big(\log(9), 1\big)\big]\end{align}$$

$$\begin{align}=&[1, 1, 1, 0, 0, 0, 0, 0, 0]\\
+&[0, 1, 0, 1, 0, 1, 0, 0, 0]\\
+&[0, 0, 1, 0, 0, 1, 0, 0, 1]\\
\rule{0pt}{4ex}=&[1, 2, 2, 1, 0, 2, 0, 0, 1]\end{align}$$

The sequence notation is not affected by the choice of base of the logarithm. In an alternative representation, a sample $(\log(0),1) $ $=$ $(-\infty,1)$ could be added to the signal (appended to the beginning of the sequence), giving the pleasures of zero-based indexing without affecting the rest of the convolution result:

$$[1,1,1,1]*[1,1,1,1]$$
$$\begin{align}=&[4, 0, 0, 0, 0, 0, 0, 0, 0, 0]\\
+&[1, 1, 1, 1, 0, 0, 0, 0, 0, 0]\\
+&[1, 0, 1, 0, 1, 0, 1, 0, 0, 0]\\
+&[1, 0, 0, 1, 0, 0, 1, 0, 0, 1]\\
\rule{0pt}{4ex}=&[7, 1, 2, 2, 1, 0, 2, 0, 0, 1]\end{align}$$

I don't have a real application mind, but encountered such a convolution in the context of multiplication of discrete random variables. It is also interesting that if an infinite sequence of 1's would be log-sampling-convolved by itself, then the result would be the sequence OEIS A000005, the number of divisors of positive integers.


texts: Recovering signal statistics from non-uniform sampled signal

texts:
I'm interested in estimating the mean and standard deviation of a signal that was sampled non-uniformly. Assuming I have an estimate of the signal bandwidth, what algorithms would provide estimates of these quantities?

I imagine simply weighting the values by their relative frequencies could work. But is this actually a good approach? Does it have some theoretical guarantee like how the sample mean is an unbiased estimator of the population mean?

I'd appreciate any insight or references I could turn to. Thank you.


texts:
I'm trying to reduce my magnitude spectra, which I need in octave bands. I have 5 measurements of a sound source, recorded with the same microphone at different times. I want to get an average model of my data. Therefore I do the following for each of my recordings:

Divide my data in time windows
Compute STFT
Average STFT over time frames
Sum bins that fall in a specific range, corresponding to different octave bands

Once I have my octave band spectra for my 5 files I do the average of all the 5 magnitudes.
I can't find any information about averaging different magnitude spectra.
Am I doing this wrong? Any help would be appreciated.
Thanks!


texts: Decomposition of **3D** structuring elements for morphological operations

texts:
I am struggling to implement a mathematical morphology toolset in an image processing package. I base my implementation on what I saw in MATLAB, and on several papers on the subject.

There seems to be abundant literature on morphological operation optimization through structuring element (strel) decomposition. For instance, one can get a tremendous speed bonus by using two orthogonal lines instead of a square as structuring element for dilation. Several papers give methods for optimization through strel decomposition:

Rolf Adams, "Radial Decomposition of Discs and Spheres," CVGIP:
         Graphical Models and Image Processing, vol. 55, no. 5, September
         1993, pp. 325-332.
Rein van den Boomgard and Richard van Balen, Methods for
       Fast Morphological Image Transforms Using Bitmapped Binary
       Images, CVGIP: Models and Image Processing, vol. 54, no. 3, May
       1992, pp. 252-254

etc...

However, all these publications are about 2D structuring elements. I could not find much on 3D decompositions.

Do you have any clues on how to decompose:

a 3D sphere. Not a ball that is used in 2D grayscale morphology, but an actual flat 3D sphere;
a 3D diamond?


texts: The sampling frequency for band pass signal which is having frequency from 4 kHz to 6 kHz(rectangular shape from 4 to 6 kHz), can I prefer the sampling frequency of 12 kHz(2*the highest frequency present is message signal) or 6 kHz  (2*bandwidth of message signal(6-4)kHz), can u please suggest which one to prefer?


texts: Simple problem regarding rectangular pulse shifting

texts:
So i have the following rectangular pulse:

$$h(t)=\begin{cases}
      1 &  0<t< T \\
       0 & \text{otherwise}
   \end{cases}$$

Now, i am supposed to create $h(T-t)$, and the problem is that i do it two different ways, and get two different results, and it is obvious that this is not possible. Anyway, here is what i have.

My first attempt is step by step building of the $h(T-t)$ from the $h(t)$, firstly, i "flipped" variable $t$.
so i got
$$h(-t)=\begin{cases}
      1 &  0<-t< T \\
       0 & \text{otherwise}
   \end{cases}=\begin{cases}
      1 &  -T<t< 0 \\
       0 & \text{otherwise}
   \end{cases}$$

Now, i add T to this, and i have
$$h(T-t)=\begin{cases}
      1 &  -T<T+t< 0 \\
       0 & \text{otherwise }
   \end{cases}=\begin{cases}
      1 &  -2T<t< -T \\
       0 & \text{otherwise}
   \end{cases}$$

However, i've done this directly:

$$h(T-t)=\begin{cases}
      1 &  0<T-t< T \\
       0 & \text{otherwise}
   \end{cases}=\begin{cases}
      1 &  -T<-t< 0\\
       0 & \text{otherwise }
   \end{cases}=\begin{cases}
      1 &  0<t< T \\
       0 & \text{otherwise}
   \end{cases}$$

Now, there is something wrong here, but i don't see any mistakes in any of the two cases. Any help appreciated!


texts: Questions about the stability (and stationarity) of a system and state space representations

texts:
i'm pretty new to the topic and I'm trying to understand how to determine the stability of a process. I'm giving this discrete-time stochastic system:

$$
\cases{ s_t = 2s_{t-2} + 3w_{t-2} \\
y_t = s_t + v_t \\
w_t = wn(0, q) \\
v_t = wn(0, r) \\
w_t \perp v_t
}
$$
where $w_t = wn(0, q)$ means that the process $w_t$ is white noise with 0 mean and variance $\sigma^2 = q$ and $w_t \perp v_t$ means that the processes are uncorrelated.

I'm asked if the process $y_t$ is stationary and I'm finding some difficulties in deriving the answer.

I (think I)'ve written $s_t$ in terms of his transfer function (interpreting $z$ as the lag operator):

$$
s_t = \frac{3z^{-2}}{1 - 2z^{-2}} w_t
$$

This transfer function has poles in $\pm \sqrt{2}$, outside the unit circle, hence I'd say that this process is not stable and that $y_t$ cannot be stable, either.

Regarding the stationarity of $y_t$: do I have to solve the discrete SDE and evaluate mean and variance of $y_t$ or is a better and smarter way to answer the question?

Can I give a state space representation of $s_t$ and work it out in that domain?

Thank you for your help!

Edit: @mark-leeds thank you very much, you have been very clear. I was unsure about this interpretation of an AR model but now I'm OK with it.

One last question regarding the state space representation: is it possible to express the process $s_t$ in a state space model? Something like

$$
x(t+1) = A(t)x(t) + B(t)e(t) \\
s(t) = C(t)x(t) + D(t)e(t)
$$

Thank you very much again!


texts: How can an SDR recover a high-frequency signal?

texts:
How can software-defined radios operate at high frequency?

The Nyquist rate dictates that you need to sample at twice the frequency to fully recover the signal. If my signal of interest is modulated at $2.5 \,\rm{GHz}$, doesn't my SDR need to sample at $5 \,\rm{GHz}$ to fully recover and demodulate the signal of interest?


texts: Reducing noise from the same frequency band as signal

texts:
Sorry for the vague question (as I'm not even quite sure what I want to do is possible), just asking for some general direction to take my research.

For a brief description, my signal resembles exponentially distributed noise, bandlimited, with a triangle shaped spectrum around the baseband. I'm looking to remove additive white Gaussian noise within the same frequency band as the signal, ideally without affecting the signal bandwidth.

Below are some simulated examples of a section of my signal and noise:

and their bandwidth

So the question is whether I can, somehow, attempt to improve the SNR, given that they share the same frequency band, and my signal is not known beforehand (though it differs from the noise in distribution and spectral shape).

EDIT Averaging successive acquisitions is not an option, as what we want to measure is a local time delay of the blue trace. As such, just averaging would interfere with the measured delay.

EDIT 2 I see people commenting that the SNR is fairly high in this case, just want to clarify: This is the current use case and we are currently happy with the performance.
In the case of our system, without going too much into irrelevant details, if we increase the spatial resolution the SNR will deteriorate (up to 10 or 20dB), so maybe the representation isn't accurate of what I was trying to convey.
I will post new figures as soon as I am back in the office, or if I get an opportunity tomorrow of a noisier case.


texts: Is it possible to extract peak locations in the time domain using help from fourier/wavelet analysis?

texts:
The signal I'm studying has fundamental frequencies of 20 and 60 cycles per minute (shown in the Periodogram graph). It is straight forward to extract the peaks in the time domain belonging to the 20 cycles/minute frequency (both are circled), however the peaks associated with the 60 cycles per minute frequency can be noisy.

6 seconds of the time domain signal

Both signals in the time domain will ALWAYS overlap so low pass filtering isn't completely necessary but I'm just curious to know if there's a better way of extracting the 60 cycles per minute peaks (aside from peak-to-trough analysis) using frequency domain information?

Thanks!


texts: What are the options for turning a bit stream into a byte stream?

texts:
I am creating a system (modulator and demodulator) for broadcasting binary data through an analog signal. I'm using BPSK with great success – at least, I'm generating great graphs.

Green is input signal, red is carrier phase, pink is digital zeroes and ones.

My question is – now that I have that pink stream of binary data, how to delineate bytes? I have been reading about clock recovery, 8b/10b, Manchester, NRZI... And I'm a little confused about what those offer me.

I imagine I could just come up with a system like Manchester, where each 0 is just 01 and each 1 is 10 – and then use a pattern like 0000 to mean "start of next byte" – but this doesn't seem right. Under that system, how would you interpret 10000010? It could mean the last byte ended with a 1 and the next byte starts with a 1 – or it could mean the exact opposite – the last byte ends with a 0 and the next byte starts with a 0.

I know there must be some basic DSP concepts that could help me here, please point me to some resources!


texts: Causal downsampling of a small buffer of GPS positions

texts:
Not sure this is the right place to ask, but I hope so.

Context

A GPS receptor is installed in a moving vehicle, and saves the measured positions with a resolution of at most $r_0$ to a small (a few dozen positions) buffer. The buffer is small due to memory limitations.

The buffer is kept as empty as possible while the vehicle is connected to a network. However, the network coverage has several gaps large enough to allow the buffer to be filled. As a result, positions measured after the buffer is filled are currently thrown away.

This causes our main use case, which is to display and analyse routes, to present large jumps between the edges of a network shadow zone. It would be better if the positions in or close to the shadow zone are saved in a smaller resolution than not at all.

Specifics

So the current approach is to downsample the buffer when it's full, in ordere to make more room for new positions. Since the sampling rate across the buffer is not constant, we just find the chunk which has the highest resolution and cut an element of it.

Since this happens in real-time, the downsampling is causal (we only know the past points), and thus selecting the best point in this buffer is not looking trivial (to me at least).

Question: are there general methods for this problem or should I continue via the ad hoc approach?


texts: What are the drawbacks of Harris corner detector?

texts: What are drawbacks of Harris corner detection when trying to match features in a stereo pair? I notice that whenever I use Harris corner detector to match features, images with strong occluded areas are matched incorrectly.


texts: Does the 'FFTW' have a cost of k*N*log_2(N)?

texts:
@Michael C. Grant wrote that "the cost of FFTW isn't an easy formula based on log_2 anymore."

But the Wikipedia article says that the FFTW "can compute transforms of real and complex-valued arrays of arbitrary size and dimension in O(n log n) time."

So can it compute 1 forward FFT with the size 1920*1080? in kNlog(N)/log(2) =

5*1920*1080*log(1920*1080)/log(2)=
217559066 real operations?

Or should I perphaps take a different value for k?


texts:
In the documentation for the [XAudio2 filter parameters] it mentions

Filter radian frequency calculated as (2*sin(pi*(desired filter cutoff frequency)/sampleRate))

Where does this formual come from? What is "Filter radian frequency" and "filter cutoff frequency" and what is the difference. Can someone translate this terminology to mathematics of DFT?


texts:
I am beginner and I want a help to know how to detect hand in Matlab.

First, I need code which detects skin and then draw a rectangle around a region so it show to me face and hand.

When I try to make rectangle only around hand and the hand is away from the face, the the rectangle is drawn around the face instead.

What do I need to do to discount the face and only draw the rectangle around the hand?


texts: What is facebook's image compression algorithm?

texts:
For a project I'm working on (master thesis), I need to figure out what is facebook's image compression algorithm.
The goal: to be able to mimic facebook's image compression locally (using python).
The images I use have the following exif: Image exif (Q = 90)
After facebook compression, I have the following exif : Compressed image exif
Note: If the images I uploaded with the link seem unavailable, please refresh the page! :)
After searching for a while, I could not find a recent and clear answer to my question.
I tried to run the compressed image through ImageMagick's identify. It tells me that it is a JPEG compression (DCT progressive + Huffman) with a quality of 85. But when I try to do my own compression using python, the result is not the same.
Here is my code:
from PIL import Image
img = Image.open('original.jpg')
img.save('output.jpg', optimize=True, quality=85, progressive=True)

Do you have any ideas or resources that could help me? I'm running out of ideas.
Thank you!


texts:
I've been studying image processing methods (without signal processing background) and many things are interesting to me. I was think about noise. In many definitions noise is unwanted  and without meaning  electrical or electromagnetic energy that degrades the quality of signals and data. Considering noise a signal type. So my question is

Noise is undesirable signal, or incomprehensible signal?

Note: Sorry by the philosophical question.


texts: Analogue Filters (Picking the right filter!)

texts:
I'm new to Analogue and I'm trying to find the Ideal filter for the following specifications:

$0.1\ \rm dB$ ripple
Cutoff Frequency $F_c = 80\ \rm{kHz}$
An attenuation of $40\ \rm dB$ at the stop-band frequency $F_s = 130
\ \rm kHz$

I'm not sure if I'm on the right path but till now with all the research I've done I think I have to use a Low pass, 7th order Chebyshev Filter (Type 1) using Sallen key topology (not sure why but I was thinking since it have less noise at a gain of 1 I should go with a Sallen-key).
Any help and explanation would be great!
Thanks in advance :)!
Pix


texts: How is the formula of the cutoff frequency of a 2nd order Butterworth low-pass filter derived?

texts:
According to articles and notes, the cutoff frequency of a 2nd order Butterworth lowpass filter:

is expressed by the formula:
$\omega _c = {1\over\sqrt{R_1R_2C_1C_2}}$

and since $\omega _c = 2 \pi f_c$, alternatively

$f_c={1\over2\pi\sqrt{R_1R_2C_1C_2}}$

Where does this formula come from? I'm asked to write the equations and formulas which eventually lead to the formula of the cutoff frequency, but I've been very confused and tutorials and notes don't help at all. Thanks in advance.


texts:
I have generated 16 chirps in experimental study:

and the signal was picked up by some electronic sensor, and then I have calculated the spectrum of the received signal:

Following that I have calculated the cross spectrum density between the sound emitter and receiver and did some post processing, i.e.:

x=voltage1(:,2);
y1=vertgeo1(:,2);

Fs=1000;
window=hamming(512);
noverlap=512/2
nfft=512;Fs=1000;

[Pxy1,f]=cpsd(x,y1,window,noverlap,nfft,Fs);

From this code the number of the points that shows the frequency resolution, essentially the number of the point that f contains is $1\times 256$.

I would like to know how I can increase this to $1025$ points.


texts: Getting HRTF from 3-D scan + acoustic simulation

texts:
This question follows from Modelling propagation of sound wave by particle simulation

This last fortnight I have started experimenting with Binaural audio; it really blows everything else out of the water.

The main reason it hasn't taken off, I think, is that every individual needs to get their own HRTF calculated, which consists of sitting in a soundproof room with microphones in your ear canals, while something makes popping noises at thousands of locations around you.

I'm only aware of one service provider that may calculate your HRTF: http://www.physiol.usyd.edu.au/~simonc/hrtf_rec.htm  and they are in Australia!

I'm wondering whether a more practical method may emerge, which would involve taking a 3-D scan of one's head (maybe using something like http://www.david-3d.com/) and shipping it off for a heavy dose of distributed computing which would return HRTFs for that individual.

I wonder if we may see a day where people have their own HRTF data stored in the cloud, and can enjoy binaural sound on their iPods.

How far away is such a technology?  And other any other contenders for measuring HRTF?

π

PS could we possibly have some more tags like "binaural", "hrtf/hrir"


texts: Signal values we will 'miss' between sampling instances during sampling of band limited signals

texts:
According to the Nyquist–Shannon sampling theorem, any continuous time signal with a bandwidth $B$ smaller than Nyquist frequency $f_N=f_s/2$ (with $f_s$ the sampling frequency), which is sampled at sampling frequency $f_s$ can be perfectly reconstructed by sinc interpolation (i.e. the Whittaker–Shannon interpolation formula).

Assume we sample an unknown, bounded in magnitude, continuous time signal with constant sampling time $T=1/f_s$ at sample instances $kT$, ($k\in\mathbb{Z}$), without sampling jitter or quantization. We add the constraint that $B=\alpha f_N$, with $0\leq\alpha\leq 1$.

What I would like to figure out is the following:
At sample instant $k$, I want to determine for each $\alpha$ a worst case fractional 'overshoot' of any continuous time signal between samples $k-1$ and $k$, that I could have had. I.e. how much the continuous time signal was higher than the highest (absolute) sampled values at sample instants $k$ and $k-1$. The continuous signal, or reconstruction (since the sinc interpolation is perfect!!), that we have 'missed' by sampling.

Example:
We set $\alpha=1$ and assume discrete time signal [1,0,1,0,1,1,0,1,0,1] (notice the double 1 near the middle, and does this signal even have $\alpha=1$?). Its sinc reconstruction (blue line) from the samples (black impulses) looks as follows (I have plotted the sincs belonging to each sample in gray):

The 'overshoot' between samples $k=0$ and $k=1$, is $\approx 0.7$ or $70\%$. So we missed a peak of value 1.7 in our original band-limited continuous time, or the 'perfectly band limited reconstructed', signal. If I would have put 3 or more consecutive 1's the overshoot would have been less (the Gibbs phenomenon is in the end much smaller). Therefore, 2 consecutive continuous samples like this is 'worst case'.

Extending the signal in both directions will make the overshoot grow:

Which shows relative overshoot of $\approx 1.1$ to a value of almost 2.1.

For any sequence length $2m$, this 'overshoot' $o(m)$ will grow indefinitely, $o(m)\propto\ln{(m)}$, which goes to $\infty$ when $m\to\infty$. This is because each sample the sincs will create constructive 'interference', and the sum of $1/\pi n$ (the contributions of all the envelopes of unit sinc's) for $n\to\infty$ does not converge.

This (I think) simlar to the following: if constantly sampling a value 0, I could also reconstruct a continuous time signal with infinite amplitude that is only sampled in the nodes at values of 0, e.g. $\sin{\pi f_s t}$. This tells me the same thing: that if I allow a signal to be at the Nyquist frequency the worst overshoot I could 'miss' is infinite.

We can now state that $o(m\to\infty)|_{\alpha=1}=\infty$. And we can reason that that $o(m\to\infty)|_{\alpha=0}=0$ (sampling a constant signal of which you know that it is band limited has a unique constant reconstruction).

What if $\alpha<1$?

If we now assume we do this same sinc interpolation, but know for sure $\alpha<1$, like $\alpha=0.5$. Then, (my gut feeling says) this effect should go down and should even remain finite (when $m\to \infty$)!. Since for any signal brick-wall limited to bandwidth $\alpha f_s/2$, we get a filter impulse response of $h(t)\propto \text{sinc}{\left(\frac{t-kT}{\alpha T}\right)}$ (right?). Therefore, signal transitions cannot be as fast as for the changing impulse train example above, and therefore contributions of each sinc function during reconstruction cannot create infinite constructive interference.

My problem:
I don't know how to proceed from here; how to form a 'proof' of the worst case overshoot I could ever have find between 2 consecutive samples, knowing that $\alpha<1$, for $any$ signal (not necessarily these unit impulse-train like examples). A given value for $\alpha$ gives me a slope $\frac{dh(t)}{dt}$ of the band-limiting convolution kernel $h(t)$, which should tell me something about how much consecutive samples need to be different, but I do not see the steps to take from there to reach any generic conclusion.


texts: How to find energy of sinc function $x(t) = \frac{\sin \pi t}{\pi t}$ with out the help of Fourier transform or Parseval's theorem?


texts: What is the difference between Remez exchange and Parks-McClellan filter design?

texts:
A question recently came up regarding Parks-McClellan and some of the comments point out that the wikipedia article on Parks-McClellan states,

  ...the Parks-McClellan algorithm is a variation of the Remez algorithm or Remez exchange algorithm.

Is there a difference in performance?  Is the remez function in octave equivalent to the Matlab function firpm?


texts:
Given this image:

How can  I extract the 'rectangle' my arm is holding  up?
I know about closing, dilation and erosion etc but I still manage to get artefacts of my hand included after filtering.

The size of that rectangle can vary due to varying distance if the camera.

I am using emgucv:

CvInvoke.Threshold(MyImage, MyImage, 10d, 255, ThresholdType.Binary);
var element = CvInvoke.GetStructuringElement( ElementShape.Cross, new Size( 3, 3 ), new Point( -1, -1 ) );
CvInvoke.Erode( Threshold, Threshold, element, new Point( -1, -1 ), 1, BorderType.Reflect, default(MCvScalar) );

... and I would like image to 'roughly' look like this:

Any ideas?


texts: Natural, Forced and Total System Response - Time domain Analysis of DT

texts:
I am a self-learner of DSP by reading related books.

While reading Time domain Analysis of DT (Discrete Time) system, I suddenly got introduced to difference equation representing a system. It then further stated -

The response of any DT system can be decomposed as :

Total Response = Zero State Response + Zero Input response

Where Zero state response is system response due to input alone and zero input response is when input is zero.

It goes further on to explain Natural, Forced and Total System Response - but unfortunately I am unable to map it to a concept.

For me its like response of a car engine (example per say or anything - a system) is like response when engine is not started + when engine is started - but this does not makes any  sense for me. The engine system why should it consider zero input response at all. Its just a change of its state.

Similarly I am trying to understand Natural / Forced response but do not get the concept quiet clear. Also What is the role of studying them and where can we apply these concepts?

What is meant by recursive and non-recursive system - though my understanding is when system output is fed again into system as its input, it would then make the system call as recursive system - but why does a system need to do that?


texts: Example exercise - Shannon's law on PCM-digitized signal

texts:
I am studying for a Digital Communications course in my university, and I am presented with the following example exercise.

  You have a sound signal with spectral components in the 300-3000 Hz range. The signal is going to be digitised using Pulse Code Modulation.

  a. For SNRdb == 30dB, what is the number of uniform quantum levels needed?

  b. What is the minimum data transfer rate needed for the transmission of the digitized-with-PCM signal?

I know that for both of those questions I need to use Shannon's law, C = B*log_2(1+SNR) with SNRdb = 10*log10(SNR).

For task a. I need to find M, for task b. I need to find C.

My question: What do I use as B in both of the questions? How do I utilize the digitized-with-PCM info?

Do I take a median 1650Hz in a?

In b, do I write that PCM uses from 8 to 20something bits per sample, and use B = 20*3000Hz? Or do I just use the maximum B of the range, 3000?

Any help would be very much appreciated.


texts: zero-phasing frequency components while keeping the same magnitude, in Matlab

texts:
How is it possible?

I was thinking of taking just the real part of the DFT of my signal (isn't that zeroing out phases?) with

real(fft(X))

but the magnitudes of that result don't match :

abs(real(fft(X))) != abs(fft(X))

Thanks in advance.


texts: Why would we simulate baseband rather than passband

texts:
My question is, Why do we usually use baseband in simulation rather than passband? however, if we want to test our simulated idea in real environment, it must be sometimes in passband !!

So, should simulation in baseband give the same results if we made it in passband? for term of BER, MSE and so on

Here is an example as asked in the first comments,

In baseband, Suppose you have ofdm system, with $s$ is the $n^{th}$ block with N symbol, when performing N-by-N FFT to get $x = FFT^Hs $, after adding a CP as guard interval with length P, we will have $M = x + P$ is the transmitted signal as baseband.

But in Passband, we will not transmit the $M$ as it's, we will have additional process when having $M$, this process is summarized  conversion --> matches filter --> carrier modulation by multiplying (.*exp(-j*2*pifct), fc is the carrier frequency t is the time period, and them transmit real part of this results, then in the receiver side, we follow the inverse of these steps till getting the baseband signal $M$ then, both passband and baseband will follow the same steps.

Thank you


texts: Virtual Microphone using Mic Impulse Response

texts:
I am working on a virtual microphone project which tries to emulate the output of an audio signal that would be obtained using another microphone. Let me know if this sounds reasonable as a solution or not.

Let's suppose, we have two microphones: Mic1 and Mic2. A sound source s is played through a loudspeaker using Mic1. Let h1 and h2 be the impulse response of Mic1 and Mic2 obtained using the same loudspeaker respectively. Now, if we record using Mic1, we obtain:

s1(t) = h1(t) * s(t) % (*) operator is convolution.

If we want to achieve s2 (the sound produced by Mic2), we can do the following:

% Take fft of h1(t) & h2(t)
H1f = fft(h1); H2f = fft(h2);

% Obtain S(f) using division in fft
S1f = fft(s1);
Sf = S1f/H1f;

% Obtain s2(t) using multiplication in freq domain & ifft
S2f = Sf .* H2f;
s2 = ifft(S2f);

Does this sound feasible as a solution to the problem? A more practical implementation (suggestion) would be appreciated.


texts: What language or framework would you use to develop and test natively algorithms meant to be ported on an embedded hardware later?

texts:
I am trying to find a nice clean solution to this problem since some days now, and couldn't make myself to stick to a solution.

I need to develop some stream and block processing for audio, including cascaded biquad, FT and iFT, FIR, dynamic correction. The whole audio chain will be embedded either in a DSP or on a ARM running Xenomai Linux.

What I am looking for is a way to quickly compile and test bits of codes, getting audio in and out, and in a language that will allow me to target a broad range of devices. (Would be C or Cpp surely).

I don't want to choose the definitive platform now because I don't know yet how much processing power will be required to make my full project runs.

Here are some tought:

I am using Gnu Octave but afaik it will not compile a real-time pipeline and the code is hardly portable.
I tried Audio Weaver from DSP Concepts (which is great) but the targets DSP are limited and the engine is proprietary.
I looked at Faust which is a very promising framework for audio (even able to generate Cpp code), but unfortunately it was designed with stream processing in mind and not meant to process blocks of data
I looked at some LabView + FPGA solution but the price and licensing is just WOW !
I looked at some IDE from Ti and AD but there also the price is too high.
Using a simple VST or ALSA framework might be ok but I don't really know where to start.


texts: Independent Subspace Analysis and Topographic ICA

texts:
Can anyone recommend the best references for independent subspace analysis (ISA). I have read 'Multidimensional Independent Component Analysis' by Cardoso and would like to know of any modern, solid (i.e. generally accepted to work) approaches.

If you have implemented ISA before how did it perform?


texts: How to detect continuous noise in audio call?

texts:
We want to detect bad audio quality in audio call.

Should we extract this information from amplitude or frequency spectrum ?
Is there any technique/transform/filter available for this ?

Example:

I am attaching two signal graph(Clean Signal vs Noisy One) of Amplitude vs time.

Clean Signal Graph:

Noisy Signal Graph:

As shown in above graph, there is a continuous noise in background in Noisy signal.

How can it be differentiated between this two signal ?


texts: Using MATLAB Function `mpiir_l2()`: Results Are Not IIR but FIR Filters - Why and How to Avoid This?

texts:
I'm using mpiir_l2() from user Matt L's PhD thesis to design IIR filters. I set the number of numerator and denominator coefficients both to the same value (between $100$ and $300$). The maximum pole raduis is $0.98$. My observation is that the result is effectively a FIR filter, with the denominator coefficients allways being $[1, 0, 0, ...]$. Basically, the result is equal to that of the lslevin() function for least squares FIR design from the same source.

Clearly something in my filter spec makes the FIR solution have better errors than any IIR filter. The specification is the result of numerical optimization and looks pretty random to a humans eye, so I'm not sure what properties to look for.

Is there a way to get a true IIR filter from this function? I want to compare IIR anf FIR filters for the same specification and see what's better. Right now, I can't really compare when the results are basically the same.

EDIT: Here is a sample script and corresponding desired response. It seems like I have to use a third party hosting site for this, sorry.

close all;
clear all;

% loads the specification, weights and frequencies as three vectors
% the sepcification has a relevant band that is padded with zeros to both
% sides. There is a "don't care" gap between the relevant band and the
% zeros.
load 'desiredResponse.mat' % loads Ws, freqs, desiredResponse

% search area
allowedOrder = 1:1:3;       % numerator and denominator order
allowedShifts = 0:0.125:2;      % linear phase term

% other variables
Fs = 48000;
normalFreqs = freqs / Fs * 2 * pi; % normalized frequencies
poleRadius = 0.98;

% preallcoate for the search
lowestError = inf;
bestNum = [];
bestDenom = [];
errors = zeros(length(allowedOrder), length(allowedShifts));
orders = zeros(length(allowedOrder), length(allowedShifts));
shifts = zeros(length(allowedOrder), length(allowedShifts));

% search for best aproximation
for nl = 1:length(allowedOrder)
    ord = allowedOrder(nl);
    for ns = 1:length(allowedShifts)
        orders(nl, ns) = ord;
        shifts(nl, ns) = allowedShifts(ns);
        % add the linear phase shift
        phaseShiftTerm = exp(-1i * shifts(nl, ns) * (ord-1)/2/Fs * 2 * pi * freqs);
        D = desiredResponse .* phaseShiftTerm;

        % design the filter (change denominator order here, if you like)
        [num,denom] = mpiir_l2(ord, ord, normalFreqs, D, Ws, poleRadius);

        % calculate filter response
        H = freqz(num, denom, freqs, Fs);
        % calculate the weighted error
        errors(nl, ns) = sum(Ws .* abs(H - D).^2);

        if errors(nl, ns) < lowestError
            lowestError = errors(nl, ns);
            bestNum = num;
            bestDenom = denom;
        end
    end
end


texts: GPS CDMA signal acquisition (BPSK modulation)

texts:
Sorry, this may be a rookie question: I have a CDMA gps communication in BPSK. Codes have a duration equal to $1\ \text{ms}$ and their length is 1023. At the carrier frequency ($f_c = 1.57524\ \text{GHz}$), Doppler frequency $f_d \in [-5, +5]\ \text{kHz}$.

The receiver is depicted here. The mixer shifts the signal to intermediate frequency $f_i = 1.3\ \text{MHz}$.

Sampling of the A/D converter is at $\frac{1}{T_s} = 8\ \text{MHz}$ and samples are quantized to 1 bit for I and 1 bit for Q. As a result, I have 24000 samples for each channel (so they should "cover" $3\ \text{ms}$).

The purpose, of course, is to recognize the transmitter and its Doppler frequency.

In my understanding, it should be sufficient to recollect the transmitted symbols with $\hat{x_k} = I_k\ \cos(2\pi f_i k T_s) - Q_k\sin(2\pi f_i k T_s)$ and correlate the sequence obtained with each of the codes (better, their "BPSK" equivalent $2\ b_k - 1$) and find the one maximizing this correlation.

I understood that once obtained the data, I can estimate the phase shift between the signal and the local oscillator (thus the doppler) having $Z_k = I_k \cdot\ Q_k = \frac{1}{8}\ x_k \cdot\ \sin(2 \Delta \theta)$

Hoping that what I just wrote is correct:

are these considerations still correct now that I do have the quantized samples for I and Q and not their actual value?
since I/Q samples (and then the $\hat{x_k}$) and the codes have a different sampling rate, what's the best option to compute their correlation? I had in mind to run

ifft( fft(code, Nfft) .* conj(fft(xk, Nfft)) )

but having the two signal different sampling, I am pretty sure this will fail.


texts: How to Estimate Carrier Frequency & Amplitude for Several Overlapped Template Signal

texts:
As is shown in the picture,
a template signal is delayed and overlapped in time domain,
the shape of  template signal is already known
The positions of each template are already known.

Each template has its own carrier frequency, and amplitude.
Or the carrier has a Gaussian distribution $\mathcal{N}$($\mu$,$\sigma^2$).

How to estimate them？ Please provide some recommendations, thank you.


texts:
How does an equalizer prevent a sample and hold effect?

Does a reconstruction filter only smooth out the output of an DAC (+S&H) or does it have another function?


texts: GNU Radio Companion - Constellation Modulator block - how to generate single constellation point only?

texts:
I'm beginning to learn about the Constellation Modulator block, and just trying it out by setting up a 4-QAM system, with the constellation points arbitrarily set as : [-1-1j, 1-1j, 1+1j, -1+1j]
And the corresponding symbol map set as : [0, 1, 2, 3]
I tested this block by using a Vector Source with a value of '2', which is the third value of the symbol map. The corresponding constellation point (the third value in the list) --- which is 1+1j -- is indeed plotted.
Although, I also notice that another constellation point (which happens to be the first point in the list) - namely -1-1j is also plotted.
Is there a way to get the simulation to focus only plotting the 1+1j point? (that is - to plot a single point that corresponds to the the '2' value (ie. the 1+1j constellation point).
I've attached a screenshot showing the point that I assumed to see only (marked in green) ----- while the point marked in red seems to be plotted as extra. Any recommendations about what I can do to stop that extra constellation point (at -1-1j) from coming into that constellation plot?
As a result of that extra constellation point coming into the picture ----- I think that the constellation plot also picks up points in-between the two constellation points ----- which comes up as points along a straight line. Basically - picking the trajectory in-between -1-j and 1+1j. I'm just trying to get a plot with 1+1j appearing in the constellation plot only.
Thanks very much in advance for any recommendations/help!

=======================================
UPDATE - I have added waveform plots (obtained from a QT GUI Time Sink block) - which shows some interesting waveforms. These waveforms aren't what what we expect to see, right? Thanks in advance for more comments can recommendations about the Constellation Modulator block.


texts: Linear Time-Invariant system frequency response and impulse response given input and output in time domain

texts:
I am having difficulty with this question specifically with the $\sin()$ part. So here goes nothing:

Consider a Linear Time-Invariant system with input $x[n]$ and output $y[n]$. When the input to the system is
$$x[n]=5 \sin(0.4\pi n)/(\pi n) + 10 \cos(0.5 \pi n)$$
and corresponding output is $$y[n] = 10 \sin(0.3 \pi (n-10))/(\pi (n-10))$$

Determine the frequency response $H(e^{jw})$ and the impulse response $h[n]$ for the Linear Time-Invariant system.

Thanks in advance.


texts:
Structure Tensor is a matrix in form:

$S=\begin{pmatrix} W \ast I_x^2 & W \ast (I_xI_y)\\ W \ast (I_xI_y) & W \ast I_y^2 \end{pmatrix}$

where $W$ is a smoothing kernel (e.g a Gaussian kernel) and $I_x$ is gradient in the direction of $x$ and so on.

Therefore,size of structure tensor is $2N \times 2M$ (were $N$ is the image height and $M$ is its width). However it is supposed to be $2\times2$ matrix to decompose eigenvalues such that we obtain $\lambda_1$ and $\lambda_2$ as it is mentioned in many papers.

So, how to calculate $S$ matrix?


texts: Auto-correlation of the sum of two generic signals

texts:
Be $x[n]$ and $y[n]$ two generic discrete-time signals. Given $s[n] = x[n] + y[n]$ I want to evaluate its autocorrelation $R_s[l]$.

By definition (https://en.wikipedia.org/wiki/Cross-correlation):
$$R_s[l] = \sum^{+\infty}_{m = -\infty}s^*[m]s[m+l]$$
Which expands to:
$$\dots = \sum^{+\infty}_{m = -\infty}(x[m]+y[m])^*\,(x[m+l]+y[m+l])$$
Since $(a + b)^* = a^* + b^*$ :
$$\dots = \sum^{+\infty}_{m = -\infty}x^*[m]x[m+l]+y^*[m]y[m+l]+x^*
[m]y[m+l]+y^*[m]x[m+l]$$
Using the definition again:
$$\dots = R_x[l] + R_y[l] + R_{xy}[l] + R_{yx}[l]$$

Is this result correct?


texts: Why does Phase-locking value show spurious effects at the beginning & end of a signal?

texts:
I am trying to understand the behaviour of the phase-locking-value, as e.g. defined here, by some simple examples.
Basically, I am creating two signals with the same frequency & amplitude, where the phase-offset between the two is random. Now, as expected, at most time points the PLV gives me a low value indicating that there is not a consistent pattern of phase onset. However, at the beginning & end of the signal, there are some strange effects. I suppose that these might be somehow introduced by the Hilbert transform, or is this a problem with my implementation?
Here is how I defined the function:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

def plv(x,y,trial_axis=0):

    hilbert_x = signal.hilbert(x)
    hilbert_y = signal.hilbert(y)

    normed_x = hilbert_x/np.absolute(hilbert_x)
    normed_y = hilbert_y/np.absolute(hilbert_y)

    p = np.absolute(np.mean(normed_x*normed_y.conj(),axis=trial_axis))

    return p

Here we have a plot of the behaviour:
t = np.linspace(0,10,1000)
N = 50
s_1 = np.array([np.sin(10*2*np.pi*t) for i in range(N)])
s_2 = np.array([np.sin(10*(2*np.pi*t)+np.random.uniform(low=0,high=2*np.pi) ) for i in range(N)])
p = plv(s_1,s_2)
plt.plot(p)


texts: What's the correct graphical interpretation of a series of MFCC vectors?

texts:
I'm studying speech-recognition, in particular the use of MFCC for feature extraction. All examples I've found online tend to graph a series of MFCC extracted from a particular utterance as follows (graph generated by me from the software I'm writing):

As you can see in the graph above:

the x-axis is used for each of the MFC coefficients (from 1 to 12 in this example)
the y-axis is used for the values of the coefficients (ranging aprox from -12 to 42 in this example)
you have as many lines as frames or feature vectors you have extracted (140 in this example).

Now, this doesn't make too much sense to me, because what we are seeing here, is the superposition of all the feature vectors at once, losing any time information. I'm having a hard time to understand how this representation is useful.

In my mind, I would represent the extracted vectors as follows (again, graph generated by me):

In the graph above:

the x-axis is the frame or vector number (1 to 140)
the y-axis are the coefficient values (again, from -12 to 42 aprox)
you have one line for each feature (12).

To me, this representation should be more useful because you can see the evolution in time of each particular feature, and in my mind that should have a stronger impact on how to apply comparison algorithms on spoken words.

Maybe the two representations are equally valid and useful for different purposes, very much like when you need to study a signal in the time domain or in the frequency domain, but in the case of speech recognition I would expect the evolution in time of each individual feature to be more meaningful than the density of values for each feature (and perhaps I'm completely wrong :P).

So, two questions in fact:

Why is the first representation the one than seems to be widely used and not the second one?
When you want to compare two sets of extracted MFCCs, for example by using Dynamic Time Warping - DTW, and related to this topic, do you compare the feature vectors (i.e. 140 vectors of 12 features) or the frames (12 vectors of 140 frames)? (in other words, MxN or NxM?)

Thanks!


texts: High frequency noise out of a dummy phase vocoder

texts:
I am trying to understand where does the high-frequency noise comes from when I use a phase vocoder with no modification to the spectrum or speed.
Here is the MATLAB code I used. I changed it to my wave file, and turn hopout to 121 so that the vocoder is neither changing the spectrum nor it is changing duration.
Still, it appears that the vocoder is introducing a lot of noise, especially so in the high-frequency region.
Here is how the input looks like:
The Wave-form:

and the spectrogram

The output looks like this:

and the spectrogram

The concerning part is the spectrogram of the output. I don't understand why it got all red.
The code re-normalize the overall amplitude when writing this output, so I wonder the effect might be caused by it. However, when the input is re-normalized using the same formula, it yields the same spectrogram, so normalization is not the reason.


texts: How to calculate fourier coefficients of sum of two discrete time with different fundamental periods

texts:
Assume that we have two discrete-time signals named x[n] with fundamental period 3 and fourier coefficients ak (k from 1 to 3), and y[n] with fundamental period 5 and fourier coefficients bk (k from 1 to 5). Now we want to calculate fourier coefficients for z[n] which is defined as below:
z[n]=x[n]+y[n]
How should we do this knowing nothing else about z[n]?


texts: Comparing multiple signals for similarity

texts:
I have multiple (between 2 and 100) signals and need to determine when a significant number diverge from the rest. We're exploring machine learning techniques, but we also want tackle this as a signal processing problem and see where we get the best results.

This very informative post suggests that best results come from a weighted ensemble of techniques, including:

Similarity in time domain (static): Multiply in place and sum.
Similarity in time domain (with shift*): Take FFT of each signal, multiply, and IFFT. (matlab's xcorr)
Similarity in frequency domain (static**): Take FFT of each signal, multiply, and sum.
Similarity in frequency domain (with shift*): Multiply the two signals and take FFT. This will show if the signals share similar spectral shapes.
Similarity in energy (or power if different lengths)

But this is a fairly high-level outline. Can anyone point me to a more thorough discussion of these techniques, preferably with some python code or in lieu of that, some code in R?


texts: Interpolation formula for two dimensional signal reconstruction in the frequency domain from polar samples

texts:
In the book, Advanced Topics in Shannon Sampling and Interpolation Theory by Robert J. Marks II, one may find an interpolation formula for reconstructing a two dimensional signal from regular polar samples in the frequency domain:

  Let $f(x,y)$ be space-limited to $2A$. Let its Fourier transform in polar coordinates $F(\rho,\phi)$ be angularly band-limited to $K$. Then $F(\rho,\phi)$ can be reconstructed from its polar samples via

  $$F(\rho,\phi)=\sum_{n=-\infty}^\infty\sum_{k=0}^{N-1}\widetilde{F}\left(\frac{n}{2A},\frac{2\pi k}{N}\right)\operatorname{sinc}\left[\frac{2A(\rho-n)}{2A}\right]\frac{\sin\left[\frac{1}{2}(N-1)\left(\phi-\frac{2\pi k}{N}\right)\right]}{N\sin\left[\frac{1}{2}\left(\phi-\frac{2\pi k}{N}\right)\right]},$$
  where $N$ is assumed even and
  $$\widetilde{F}\left(\frac{n}{2A},\frac{2\pi k}{N}\right)=\begin{cases}F\left(\frac{n}{2A},\frac{2\pi k}{N}\right),&n\ge 0,
\\
F\left(-\frac{n}{2A},\frac{2\pi k}{N}+\pi\right),&n<0.\end{cases}$$

I figured that the argument of the cardinal sine function contained a typo, since the $2A$ in the numerator and denominator would appear to cancel. Thus I would guess that the author actually had

$$\operatorname{sinc}\left[\frac{\rho-2An}{2A}\right]$$

in mind. Correct me if I'm wrong?

Then the last quotient with the sine functions looks very similar to the Dirichlet kernel:

$$D_N\left(\phi-\frac{2\pi k}{N}\right)=\frac{\sin\left[\left(N+\frac{1}{2}\right)\left(\phi-\frac{2\pi k}{N}\right)\right]}{\sin\left[\frac{\phi-\frac{2\pi k}{N}}{2}\right]}$$

but is not quite the same.

Incidentally, my attempts to implement the interpolation formula in the theorem via Matlab for specific examples have failed completely with regards to reconstructing a signal. Perhaps someone can point me in the direction of a correct formula for signal reconstruction in the frequency domain from polar samples?


texts: Designing a fast linear operator with $\pm 1$ entries with low condition number and low Hamming distance between consecutive rows

texts:
I need to design a matrix for compressive imaging where each row represents an $N$-pixel filter in a focal plane through which light is masked, summed, and measured (think of Rice's single-pixel camera). The requirements are the following.

There is cost associated with switching elements of the mask. Thus, the desire to limit the Hamming distance between consecutive rows/measurements (e.g., $25\%$ sign flips).
A low condition number is desirable from an isometry perspective.
It should be constructed, at least in part, with fast operators (e.g., Hadamard-Walsh transforms, circulant matrices) to speed iterative reconstruction.

Note that the constraints of minimizing the condition number and the Hamming distance are conflicting. For example, any set of orthogonal rows, as one would have with a Hadamard matrix,  minimizes the condition number ($=1$), but all rows necessarily have a Hamming distance of $\frac N2$ (exactly half the signs flip between any pair of rows).

Does anyone have any ideas about prior art or insight how to approach this problem from a principled perspective?

Update: circulant matrix approach

As Thomas Arildsen notes in the comments,  if you create a circulant matrix from a length-$N$ vector $v$ that has $K$ sign changes, each row will vary from the previous row in exactly $\frac KN$ spots.

Furthermore, the matrix spectrum is equal to $\mbox{FFT}(v)$. This allows quick evaluation of condition number.

The problem is thus how to create the vector $\vec{v}$ with $K$ sign changes that minimizes the objective

$$\vec{v}^* = {\text{argmin} \\ \vec{v} \in \Gamma_K} \frac {\max_i
|[\mathbf{F}\vec{v}]_i|^2} {\min_j |[{\mathbf{F}\vec{v}]_j|^2 }}$$

where

$$\Gamma_K \triangleq \{\vec{v} \text{ has K sign flips}\} \subset \{-1,+1\}^N\\ \mathbf{F} \text{ is FFT}$$


texts:
We have an expression shown in red bellow for non ideal DAC ZOH,which links the sampling frequency to the amplitude loss of the input signal.
I took $f_b=1kHz$ and $f_{noise}=4kHz$
$f_{in}=f_b+f_{noise}.$
Our ADC is a sinc function with $x=10kHz$ so I did a convolution for out_put_ADC=conv(fin,sinc).
And basically that is the point i got stuck with the MATLAB code shown bellow.
If everything was correct till this point how to implement the DAC so we will witness the amplitude deterioration formula of ZOH as a function of sampling frequency?
Thanks.

x=0.0001:0.0001:1; %10KHz ADC sampling

ADC_response = sinc(x);

fb=1e3;  %1Khz frequency

f_noise=4e3;  %1Khz frequency

yb=sin(2*pi*fb*x); %1Khz frequency sine input signal

yf_noise=sin(2*pi*4*fb*x); %4Khz frequency sine input noise

y_in=yb+yf_noise;%total signal going into the ADC

ADC_out=conv(y_in,ADC_response);%signal after ADC


texts: Transfer function and difference equations: why does $H(z)$ numerator polynomial not correspond to $Y(z)$?

texts:
For a discrete time LTI system, I understand that from a difference equation description of the system in the form
$$
\sum\limits_{k=0}^N{a_k y[n-k]}=\sum\limits_{k=0}^M{b_k x[n-k]}
$$
I can determine the transfer function of the system in the form
$$
H(z)=\frac{b_Mz^M+b_{M-1}z^{M-1}+\ldots+b_1z+b_0}{a_Nz^N+a_{N-1}z^{N-1}+\ldots+a_1z+a_0}
$$
Given a rationanal transfer function, I can easily go back to the difference equation. However, since the transfer function can also be expressed as $\frac{Y(z)}{X(z)}$,
$$
H(z)=\frac{b_Mz^m+b_{M-1}z^{M-1}+...+b_1z+b_0}{a_Nz^N+a_{N-1}z^{N-1}+...+a_1z+a_0}=\frac{Y(z)}{X(z)}
$$

Doesn't this imply that the polynomial in the numerator corresponds to the output and the denominator corresponds to the input?

But we know that the numerator polynomial derived from the input side and the denominator derived from the output side of the difference equation.

I'm obviously missing something in my understanding, but I don't see how to resolve this conundrum.  What am I missing?


texts:
I have a pressure platform that provides 64x128 pixel grayscale images (1px aprox 5mm) of a subject walking on top of it.
If I capture the whole step, I can obtain the average distribution of pressure for all the frames. Just by averaging (in fact the operator used is the max of each pixel) I obtain what I need. That is because at each timestep, each point in the foot always touches the same sensor.
The difficult case is when I put this pressure platform below a treadmill. From the point of view of the pressure platform, the user is "sliding", i.e. each sensor captures information from different parts of the foot at different times.
What this means is that there is an offset $o_i$ between $image_i$ and $image_{i-1}$.
Right now I have two algorithms: one local (that works in real time, providing a rough estimation) another global (very slow, though unoptimized).
The local one is just to check the offset (with subpixel accuracy of, say, 0.1px, so linearly interpolating the images) that minimizes the SSD between both images, that is:
$$
o_i = argmin_{o} ||translate(image_i, o) - image_{i-1}||
$$
which results in images like this:

The global method is to maximize a focus measure over the aggregated image, that is, if we take $O$ as the vector of all the alignments for a whole step consisting of $N$ images then:
$$
O^* = argmax_O = focus(\sum_i^N translate(image_i, o_i))
$$
Here I use the variance of the laplacian as focus measure, and it results in images like this:

As you can see, the alignment is much better, particularly noticeable in the front of the foot (the left part of the image).
To implement it, I initialize $O$ with the alignment from the local algorithm, then iterate through each element assigning the value that maximizes the global focus of the resulting image (leaving the rest fixed). I perform multiple passes until a maximum is reached (no change of global focus respect last iteration).
The following link has a video of a whole sequence, aligned using the local algorithm: https://imgur.com/a/CMolyUt
The displacement is only in the X direction, so only a positive real number is estimated. In the case of the video, that step is 182 images, so the parameters are 182 positive real numbers.
For the global method, I just do a brute force search over all the parameters, and reduce the search space by optimizing over a finite list of parameters (e.g. $0.3, 0.35, 0.4, 0.45, 0.5, 0.55$) which I know are in the range of possible offsets (just take min/max of the local algorithm output and add some tolerance).
The question is two-fold:

Is there any other algorithm/measure that I could try in order to solve the alignment problem
Is there another way of optimizing my current global algorithm instead of the brute force approach I am following?