texts: Fixed-point Exponential Moving Average filter with given cutoff frequency

texts:
I want to implement in fixed point arithmetic an exponential moving average filter, with a specific cutoff frequency for a given sampling rate.
The formula for the filter is
$$
y_n = \alpha x_n + (1 - \alpha)  y_{n-1}
$$
Looking at Exponential moving average cut-off frequency one can derive the value for alpha for a given sampling rate and a wanted cutoff. For instance, for a sampling rate of 1000 Hz and a cutoff of 50 Hz, the value of alpha is 0.267730531659313:

Fs = 1000
f3db = 50 %This is the cutoff frequency that I want
format long;
omega3db = f3db * pi/(Fs/2)
alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313

Now I tried to implement the exponential moving average filter in fixed point arithmetic as described here:  https://www.embeddedrelated.com/showarticle/779.php but I am getting an offset that I can't explain.

Here my Arduino code:

#define PERIOD_MICROSECS 1000

static unsigned long lastRead = 0;

int analog_pin = 0;   // sensor connected to analog pin 0
int analog_input0 = 0;  // variable to store the read value;
int analog_input0_filtered_fixed_point = 0;

/*
 * Formula to compute the alpha parameter of the EMA given a desired cutoff frequency.
 * Taken from: https://dsp.stackexchange.com/questions/40462/exponential-moving-average-cut-off-frequency
 *
 * Fs = 1000
 * f3db = 50 %This is the cutoff frequency that I want
 * format long;
 * omega3db = f3db * pi/(Fs/2)
 * alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313
 *
*/

//I calculate in advance the coefficients for the filter
int Q = 12; // I can choose Q freely
//I want a cutoff of 50Hz at a sampling frequency of 1000 Hz. So I calculate the corresponding alpha with formula above
//For cutoff 50 Hz, alpha = 0.267730531659313; then alpha_scaled = round(2^Q*alpha) = round(2^12* 0.267730531659313) = round(4096* 0.267730531659313) = 1097
int alpha_scaled = 1097;
// N ~= -log_2(alpha) = -log_2(0.267730531659313) = 1.901146423365515 ~ 2
int N = 2;

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N){
  y += (alpha_scaled * (x-y)) >> (Q-N);
  return y >> N;
}

// ************************* //

void setup() {
  Serial.begin(115200);
  analog_input0_filtered_fixed_point = analogRead(analog_pin); //set EMA y for t=1
}

void loop() {
  if (micros() - lastRead >= PERIOD_MICROSECS) {
        lastRead += PERIOD_MICROSECS;

        analog_input0 = analogRead(analog_pin);
        analog_input0_filtered_fixed_point = low_pass_EMA_fixed_point(analog_input0, analog_input0_filtered_fixed_point, alpha_scaled, Q, N);

        //Check the original and filtered signals with the serial plotter
        Serial.print(analog_input0);
        Serial.print(" ");
        Serial.println(analog_input0_filtered_fixed_point);
  }
}

See above a plot with the serial plotter.

Any idea about how to solve the offset problem? What am I doing wrong? Any suggestion about how to solve the general problem of a fixed point exponential moving average filter with a specific cutoff frequency?


texts: How can I modify a sound signal based on it's MFCC

texts:
I'm working on an NN algorithm with which I want to filter noise from sound signals.

I have some clean speech sound samples and noise samples, by combining them I get the signal that I want to filter.

I get the 24 MFCC values which I feed in a NN and the output is another 24 MFCC values.

How do I modify the original sound sample using the MFCC values to actually see the difference in it?

I've read about this topic here and he does something similar, but doesn't quite explains how does he modify the sound signal using those MFCC values.
By the way, im using keras and tensorflow in jupyter notebook.


texts: Fourier Transform of sampled reflectance measurements

texts:
Goal:
Calculate the complex index of refraction ($\hat n = n +jk)$ from reflectance measurements.

Data:

Reflectance measurements for various materials were taken with a FTIR (Nicolet™ iS50 FTIR, from Thermofisher). The spectral range of the reflected measurements are from 400nm - 15um, they are reported in both wavenumber and wavelength (relationship: $\hat\nu = \frac{1}{\lambda} [cm^{-1}]) $. Note, that  the full wavelength range was captured with three different detectors.

To simplify the matter (and give access to data) I have been using reflectance measurements form the following NASA database https://speclib.jpl.nasa.gov/library. I have chosen SiO2 from the database (2um-15um) which can be navigated to by selecting:

Select Spectral Type-> Minerals

Filter by class-> Silicate

Wavelength-> VSWIR+TIR

Quartz SiO_2

Brief Theory:

Here's a short summary: https://shimadzu.com.au/sites/default/files/Appl_FTIR_Polymer_specular_reflectance_055_en.pdf

The imaginary part of an analytical signal can be found from the real part alone through the Hilbert transform, e.g. $\tilde {x_c}(t) = x_r(t) +j x_i(t)$ meaning $x_i(t) = HT{x_r(t)} $. https://epdf.pub/hilbert-transforms-volume-1-encyclopedia-of-mathematics-and-its-applications.html

Changing mindset to the optical domain, and the question at hand, noting that most optical literature calls the Hilbert transform the Kramer Kronig. The FTIR measures spectral reflectance (i.e. Intensity).

$$\tilde r = re^{j\phi} = \sqrt{R}e^{j\phi}$$

$$R = |\tilde r|^2$$

Where r is the reflectivity
R is the reflectance (This is what is measurable, intensity)
$\phi$ is the phase change at the surface caused by the absorption of the material.

Through the Fresnel equations the complex refractive index can be calculated from the below equations.

Dispersion:
$$n(\nu) = \frac{1-R(\nu)}{1 + R(\nu) - 2\sqrt{R(\nu)}cos(\phi(\nu))}$$

Absorptive index:
$$k(\nu) = \frac{-2\sqrt{R(\nu)}sin(\phi(\nu))}{1 + R(\nu) - 2\sqrt{R(\nu)}cos(\phi(\nu))}$$

So, the goal is to take the measured reflectance data and calculate the phase through the Hilbert Transform, as shown below.

$$\phi(\nu_g) = \frac{2\nu_g}{\pi} \int_0^\infty \frac{ln\sqrt{R(\nu)}}{\nu^2 - \nu_g^2}$$

where $\nu$ is the wavenumber

Since I have discrete data (over a finite range) the Hilbert transform can't be calculated directly. It's usually calculated through Maclaurin's method or the double Fourier Transform (Discrete Fourier Transform (DFT) in this case).

The approximation of the Hilbert transform using the double FT is given shown in the below equation.

$$\phi(\nu_g) = 4 \int_0^\infty cos(2\pi \nu_gt)dt \int_0^\infty ln\sqrt{R(\nu)}sin(2\pi \nu_gt)dv$$

My Question:

My question pertains to the Fourier transform (FT) of the sampled data. All the equations I see on line show the FT of evenly spaced sample data with respect to "time", which I don't know. My data is evenly sampled with respect to wavenumber $\nu$, can I still use the FT but replace the sampling time (1/$\Delta$t) with (1/$\Delta\nu$)? If yes, is there anything extra I need to consider with the respect to the reluctance data.

The last equation showing the double Fourier transform has a intergal with respect to time, but I don't have any sampling time information, is there a way around this, or some conversion?

I found the following relationships online but can't quite see how to apply them.

$$f = \frac{1}{t} ....||||.... \nu = \frac{1}{\lambda}$$

$$w = 2\pi f ....||||.... k = \frac{2\pi}{\lambda}$$


texts: Carrier and timing not properly locking?

texts:
I have designed a simple BPSK system which is operating in UHF band. My transmit data rate is chosen to be 500 kbps. My receiver has AGC, timing recovery and carrier recovery bocks.

I am using Gardner timing recovery with 2 samples/symbol. The loop bandwidth for 2nd order loop filter for timing and carrier recovery is chosen to be

My signal level varies from -30 to -10 dBm at the receiver. I am observing the lock in carrier and timing is not constant. Is it because of the loop parameters or signal level or the header bits chosen for the carrier and timing recovery?


texts: How do I perform the sifting process in empirical mode decomposition?

texts:
I am programming a voice activity detection algorithm and I found a paper that recommends the use of HHT. I have been trying to program it and understand how to calculate $m1(x)$, but the shifting process confuses me.

Could someone explain what exactly I need to do after I calculate the mean value between the two splines of the original data?

Here is what I found at this website:
https://www.clear.rice.edu/elec301/Projects02/empiricalMode/process.html

The first component h1 is computed:
$h_1=X(t)-m_1$
In the second sifting process, $h_1$ is treated as the data, and m11 is the mean of $h_1$'s upper and lower envelopes:
$h_{11}=h_1-m_{11}$
This sifting procedure is repeated $k$ times, until $h_{1k}$ is an IMF, that is:
$h_{1(k-1)}-m_{1k}=h_{1k}$

Why does it go from $h_1$ to $h_{11}$? I am really confused by this.


texts: A priori spectral efficiency of OFDM scheme

texts:
Of course I can compute the spectral efficiency as the ratio of the measured bit-rate (of data bits, no header, no preable, nor redundancy bits included) and the measured spectral bandwidth. This  would be an a posteriori approach.

However, I wonder if there exists an analytical expression for a priori determination of the spectral efficiency of an OFDM communication system.

The system uses preable, header, forward error correction (in and out), cyclic redundancy check, cyclic prefix and taper window. And the modulation by subcarrier can also be changed, as well as the number of subcarriers, and which of them are nulled or pilots.

Is there some kind of formula that includes all of these parameters and returns the effective bit rate and spectral efficiency?


texts: Decimation vs Mean in downsampling operation

texts:
Let's suppose I have a discrete signal $x[n]$ with a certain frequency $f_x$.

Now, due to various unimportant reason, I want to obtain a new signal $y[m]$ with a new frequency $f_y = \frac{f_x}{k}$, with, for example $k = 10$.

Now, there are various ways to do it, I am for the moment interested in two:

Decimation:

$y_d[m] = x[m* k + (k - 1)]$

Mean:

$y_a[m] = \frac{1}{k}\sum_{i=m*k}^{m*k+k-1}{m[i]}$

What are the differences between the two, from a theoretical standpoint?
What effect on the spectrum of $y_d$ and $y_a$?

I have limited knowledge of signal processing, so sorry for if something confuses you. A brief explanation with some references for further reading would be perfect.


texts: "Reflected added to self" - word for it?

texts:
x += x[::-1]

I thought "reflect-superposed", or "reflect-added" ... any formal term for this? e.g. with linear chirp:


texts:
I was looking for a way to measure how dirty is the document in image processing. Example images below:

a. Dirty Backgrounds

b. Dirty Backgrounds with Text

Can someone give me an initial idea to start my analysis?


texts: What is the meaning of channel capacity?

texts:
If the transmitted information signal takes values from a modulation constellation, then is there a way to know what the channel capacity is?

Is there a connection between the modulation kind and the channel capacity?
Why and when do we do $\log_2(M)$ where $M$ is the number of modulation symbols? For example, if the modulation is 16 QAM, then what does $\log_2(16)$ tell us?
In equalization, when pilot symbols are sent, we get better estimation quality. This is known as non-blind equalization. But in many research articles often it is mentioned that non-blind technique wastes a lot of bandwidth due to periodically sending the training symbols. What is the meaning of bandwidth wastage and how does it occur? Does capacity decreases with using training/ pilot symbols?

These specific information is extremely hard to extract from text books and other online tutorials. Hence I have posted these problems here, hoping to get an intuitive answer. Shall be really obliged.


texts: DT unit impulse function properties proof

texts:
How to prove the following properties of DT unit impulse function. If anyone got link to a proof please mention it. I have search on the web. But only found the properties, not a proper method of proof.

Thank you.


texts: How many 8x8 JPEGs exist with distinct DCT given a quality level?

texts:
For a given quality of image (or quality of compression), Multiple 8x8 JPEG images can have the same DCT (Discrete cosine transform) matrix, for a given quality of JPEG compression. I want to find the number of those 8 by 8 JPEG files which have different DCT codes.

I write the question in a different way:

  Assume a set $A$ of images, which is initially empty. We consider all possible DCT's for 8 by 8 Jpegs one by one, and convert each one to an image. If the image already exists in the set $A$, we throw it away, otherwise, we put it in set $A$. After we are done with all possible DCT's, how many images are in the set $A$? (what is the cardinality of $A$?)

Wikipedia has the formula for the DCT conversion. But I cannot realize from that how to find what I need. I am wondering how to solve it or if that is already solved in the literature. Any idea?


texts: How is the ROC of a transform function determined?

texts:
Suppose I have $x[n]$ and $y[n]$, and I calculate their respective Z-transforms $X(z)$ and $Y(z)$ as well as their respective ROCs.

Calculating $H(z)$ is as simple as calculating the quotient of $\frac{Y(z)}{X(z)}$, but how do I get $H(z)$'s ROC?

An example:


texts:
I am studying decision trees and machine learning. I have encountered with this word several time. Although, it seems that it not an important concept, I cannot understand this word in machine learning literature. For example:

  if the test does not significantly satisfy the criteria, then it creates a leaf and computes the prototype of the instances belonging to the leaf"
  .

what does prototype in this sentence mean?

According to Wikipedia:

  A prototype is an early sample, model, or release of a product built to test a concept or process or to act as a thing to be replicated or learned from.

But it does not make sense.


texts: How to find the normalized coefficients?

texts:
I use the below matlab function to design the unquantised coefficient set:

  [b,a]=ellip (4,2,60,2*[0.2,0.3]);

and then i use the matlab finction "tf2sos" to convert the direct form coefficients to a set of a factorized second-order sections:

  [sos,g]=tf2sos(b,a);

My result is:

The normalized version of the difference equation is:

How i can find the b and the c1,c2,...

Is the actual result in this form:

or factorized means normalized and the result is in this form:

The gain plays any special role?


texts: Can we use ICA (Independent component analysis) on a noisy signal

texts: If a receiver received a noisy signal can we use ICA to separate the noise from the signal ?


texts: Determining Causality From Discrete Impulse Response

texts:
Trying to wrap my mind around the concepts of this one...

  Consider the following impulse response $h[n]$ for a linear, time-invariant system:
  $$
h[n] =\left\{\underline{1} , -2, 2, -1\right\}
$$
  where n=0 starts at 1 for h[n]

  (a) Is the system causal? Why or why not?
  (b) Is the system stable? Why or why not?
  (c) Is the system FIR or IIR? Give a reason for your answer.
  (d) Does the system have memory? why or why not?

Based on my understanding I've concluded...

(a):  ?
(b): Yes the system is stable because the summation of |h[n]| is less than infinity.
(c): The system is an FIR (Finite Impulse System) because the duration of h[x] is not infinite. It has a finite number of values.
(d): This system is not memoryless, therefore has memory. This is because in order to be memoryless (to my understanding) an impulse response should not have values at any times other than n=0. h[n] has multiple values

Updated the question to reflect how the values relate to time. The underlined 1 in h[n] denotes the origin n=0


texts: discrete-time (slew) rate limiting algorithm

texts:
Consider uniformly-spaced samples of smooth, bandlimited signal observed in noise and subject to some artifacts (small jumps).  Physical restrictions impose a constraint on the maximum magnitude of the derivative of the noiseless, "artifactless" signal.  In an effort to smooth the signal and remove artifacts, I could impose this constraint by simply limiting the size of the sample-to-sample change (e.g., matlab, but this seems like it will introduce other effects to compromise the smoothness of the signal.

The only other option I came across are nonlinear feedback loops (e.g., see this paper.  I suppose that one could also impose derivative inequality constraints on, say, some local polynomial representation of the signal as well.  Any other ideas, references, etc?

Thanks!


texts:
Im struggling to implement / build convolutional decoder .
I've done implementation of FEC encoder which it's as shown in the photo below:

So this is the encoder that I've implemented, and I want to do for it the decoder, from where could I start and how can I implement for this encoder the decoder? (decoder is working opposite to encoder and I know that but from where do I start and how to implement that decoder as what I did in my encoder photo)
thanks for any help/assistance.


texts: Calculating the SNR of Audio Signal (Recommended Libraries)

texts:
Are there any open source packages or libraries available which can be useful in calculating the SNR(signal to noise ratio) of an audio signal.

The input will be just an audio signal and I have to calculate the SNR of that signal. Can anyone suggest any library in python/C/C++ or any other language which does this or can be useful?


texts: Time duration of a signal to determine the peak frequency

texts:
I have two questions regarding signal processing:

1) In my analysis I generate a signal in the time domain with a frequency of 0.1 Hz. I apply the FFT and then determine the peak frequency of the signal, which should be 0.1 Hz. But I notice that when I simulate my signal for 200s, the peak frequency is always close to the 0.1Hz, but when I simulate my signal for 20s, then the peak frequency is not always 0.1Hz. Is there a certain rule of thumb for the duration of the signal to have the peak frequency occur in the frequency domain?

2) Currently I am just determining the peak frequency of my signal, as the highest peak in the frequency domain. But sometimes I have a lot of peaks near each other and at times they might overlap. What are some other recommended techniques to determine the peak frequency? I am using Matlab

Thanks!


texts: Self studying, getting a quality spectrogram

texts:
So I've been self studying signal processing, and been trying to use Python to get quality spectrograms, quality spectrogram like the one that would be produced with audacity etc. The steps I have done are:

Me saying 'computer', was recorded on my mac in .m4a, then I used ffmpeg converting it into .wav

$ ffmpeg -i computer.m4a -acodec pcm_s16le -ac 1 -ar 16000 computer.wav
Imported into Python via:

sample_rate, samples = wavfile.read('/path/to/computer.wav')
Now I just do a quick plot of the samples:

plt.plot(samples)

and I get this kinda graph:

Why are my peaks capped?
Anyway, I carry on and as I understand it, I perform a Fourier transform that will take my signal, decompose it into the sine waves of different amplitudes and frequencies, and then get the real part of the complex numbers that are generated:samples_fft = np.fft.rfft(samples)
frequencies = np.abs(samples_fft)
I now plot a graph of frequencies:

This is what I guessed would be like, as vocal frequencies are typically around the range of 100 to a thousand, right?
A spectrogram is when I divide my samples into a certain amount of windows, and then take the Fourier transform of each of the windows, resulting in a time frequency graph.
last_wind = 0
windows = np.linspace(0,len(samples),1000)
spectrogram = []
for window in windows:
    if window != 0:
        samples_in_range = samples[last_wind:int(window)]
        samples_in_range_fft = np.fft.rfft(samples_in_range)
        frequencies_in_range = np.abs(samples_in_range_fft)
        spectrogram.append(list(frequencies_in_range))
        last_wind = int(window)
spectrogram2 = np.asarray(spectrogram).T
Plot the spectrogram:audio_length = len(samples)/sample_rate
audio_time_intervals = np.linspace(0,audio_length,999)
tempfreq = np.linspace(0,1,19)
plt.pcolormesh(audio_time_intervals, tempfreq, spectrogram2)

So now I have a spectrogram, but of very poor quality. but I've been watching youtube videos such as this one, how do I get better resolution and clarity so as to read stuff such as formants in vowels, i tried changing the size of my partitions, but still get very similar results, more windows get less data of frequencies in each window, and less windows get poor time resolution.

I am very new to signal processing, so any help would be greatly appreciated, such as a point in the right direction or pointers as to what I'm doing wrong, in order to get clear spectrograms so as to study.

Note, I have tried python packages such signal.spectrogram, but the results are very poor aswell. Plus I'd rather get to know the nitty gritties so I really understand what is going on.


texts: FFT vs DFT Run Time Comparison (Complexity Analysis) in MATLAB

texts:
So we were given an assignment to plot the time taken by the FFT algorithm by MATLAB and a DFT algorithm written by me in MATLAB.G

The expected output should have been the DFT algorithm following the O(n^2) complexity, however when I plot the output it was coming out to be linear.

Is there some problem with my implementation of DFT or some other reason that is causing it?

The first graph is time taken to compute dft vs N(N point fft) whereas the other is for my implementation of DFT vs N.

timedft = zeros(1,length(2:1024));
timefft = zeros(1,length(2:1024));
%DFT
sig = rand(1,1024);
for N = 2:1024
    tic;
    opdft = zeros(1,N);
    for i = 0:N-1
        for l = 0:length(sig)-1
            temp = sig(1,l+1).*exp(-2*pi*l*i*j/N);
            opdft(1,i+1)=opdft(1,i+1)+temp;
        end
    end
    timedft(1,N-1) = toc;
    tic;
    fft(rand,N);
    timefft(1,N-1) = toc;
end
figure;
subplot(211);
plot(2:1024,timedft);
title("fftvsdft");
xlabel("N");
ylabel("Time")
subplot(212);
plot(2:1024,timefft);
xlabel("N");
ylabel("Time");


texts: camera calibration change orientation of axis

texts:
I want to solve for the extrinsics by using direct linear transformation on corresponding 3D LIDAR points and 2D camera points. I already have the intrinsics.

Problem is, points behind the camera gets re-projected as well (see picture below).

So I constrain to only points "in front of the camera", i.e z > 0. The problem is, on different trials where different sets of points are used, the produced extrinsic matrix produces differing axes. Sometimes, constraining z > 0 gives the right results (centre part of image), whereas other times I need z < 0. So the question is, how do I constrain the Z axes of the camera to be sticking out of the camera?

def with_intrinsic(points2d, points3d, intrinsic):
    cam1_K_inverse = np.linalg.inv(intrinsic)

    #direct linear transformation calibration, assumes no intrinsic matrix
    assert points2d.shape[0] >= 3
    assert points3d.shape[0] == points2d.shape[0]

    A = []

    points2d_homo = []
    for u,v in points2d:
        points2d_homo.append([u, v, 1])

    points2d_homo = np.array(points2d_homo).T #columns to be data points

    points2d_inv = np.dot(cam1_K_inverse, points2d_homo).T
    assert points2d_inv.shape == (points2d.shape[0], 3)
    assert points2d_inv[0, 2] == 1

    for idx in range(points2d.shape[0]):
        x3d, y3d, z3d = points3d[idx]
        u, v, _ = points2d_inv[idx]

        A.append([x3d, y3d, z3d, 1, 0, 0, 0, 0, -u * x3d, -u * y3d, -u * z3d, -u])
        A.append([0, 0, 0, 0, x3d, y3d, z3d, 1, -v * x3d, -v * y3d, -v * z3d, -v])

    A = np.array(A)

    U, D, VT = np.linalg.svd(A)
    M = VT.T[:, -1].reshape((3, 4))

    error = get_reprojection_error(points2d, points3d, intrinsic, M)
    logging.debug("error with_intrinsic: %s", error)

    return M


texts:
I have a signal,

fs = 100;                                % Sample frequency (Hz)
t = 0:1/fs:10-1/fs;                      % 10 sec sample
x = (1.3)*sin(2*pi*15*t)                 % 15 Hz component
+ (1.7)*sin(2*pi*40*(t-2)) ;             % 40 Hz component

when I take fft of x then its amplitude and phase spectrum are:

And now if deliberately I will put the leading 10 values to zeros

x(1:10) = 0;

then as per my understanding the frequency of the complete data set is not changed so the fft would not be changed. Which is supported by the amplitude fft response.

But the phase plot is showing some strange result.

Can some help me why the response is changed. Although they both are following the same equation. Also If I have a equation cos(2*pift) and my t is from 3:Tp and before time =3sec the response is zero, then will it not have the same response as t:0:Tp.


texts:
I need to define an encoding to beacon IDs. The codes are going to be periodic but in the receiver side I can have phase shifts so I' wondering the best way to decode:
- 8 bits code groups that represent the same ID:

Group1: 00000001, 00000010, 00000100 .....  10000000
Group2: 00000011, 00000110, ......11000000, 10000001
Group3: 00000101, 00001010, .............., 01000001
...

How can I verify that a code belongs to a group? maybe the Real Part of a 8 bit DFT?
Is there a group/codes generator?


texts:
To introduce my situation: I'm developing a digital synthesizer in a form of a C++ library, working with low level APIs like WASAPI, ASIO, ALSA etc. It's probably not very practical and I'm mostly "reinventing the wheel" but my intention is to learn about digital synthesis in depth. So far I have successfully implemented basic concepts like oscillators and modulation of their properties. The next logical step is a filter.

So my question is: How does a digital filter work on this low level? How exactly does it modify the individual samples?

I understand, that this involves a lot of math. That's not a problem for me. I only need a good starting point (some sources to learn) and an intuitive explanation as all I was able to find were either analogue explanations or just formulas explained with a lot of advanced terminology that I'm not familiar with.


texts:
Recently I've found a formula relating mean average deviation (MAD) SNR to root mean square (RMS) SNR:

$$\text{SNR}_\text{MAD}  = \sqrt{\frac{2}{\pi}} \text{SNR}_\text{RMS} $$

which assumes Gaussian statistics. This is for a definition of MAD something like

$$ \text{MovingAverage}[\ |x - \text{MovingAverage}[x]|\ ]$$

for some signal $x$ measured in the time domain. My question is where does the above SNR formula come from? I'm pretty new to signal processing so maybe there's a very basic reference that describes this. Is there some intuition for why MAD SNR is smaller than RMS SNR?


texts: MUSIC algorithm for DOA with floating point input signals

texts:
I am writing a Direction of Arrival (DOA) estimator using the MUSIC algorithm in python and am having trouble getting the correct spectrogram output.
Specifically, my resulting graph comes out aliased as such:
On the left is the SYNTHESIZED signal, using the same source azimuth and elevation as the real dataset. On the right is the real-world signal.

These graphs looks like this when I am looking for 2 signals. If I am looking for only 1, they look like this:

My guess is that the reason this is happening is that I am trying to use a Hermitian, which should be reserved exclusively for complex numbers, when my signals are represented using only real values, without an imaginary one. The the real world dataset is presented as an NxM .wav where N is the number of channels (8) and M is the length of the signals (174800 or so) sampled at 44.1KHz. The data arrives in int16 format, which I then convert into floating point.
The eigenvectors I obtain from the covariance matrix look realistic, with 2 peaks, 1 twice as high as the other. The synthesized signal has only 1 peak. As seen in the graphs, if I limit to using all but the highest peak eigenvector, the spectrogram doesn't really show anything, while taking M-2 does provide peaks, but theyre aliased. I am certain this is aliasing and not a reverberation of some sort, because setting an extremely low wavelength variable for the steering matrix makes the real-world signal spectrogram perfectly symmetrical. If I shuffle microphone positions, I also never get more or less than 2 peaks.
The following are some code snippets to get an idea of how I implemented the algorithm:
The synthesized signal is created using these two functions as the core:
    def emit(shape,C):
        x = np.random.randn(shape)
        x += np.random.randn(shape)
        x *= np.sqrt(0.5)
        x *= sqrt(C)
        return x

generate(sources):
    A = array.steering_matrix(sources, wavelength)
    S = source_signal.emit((channel_count,n_snapshots),1.0)

    Y = A @ S
    N = noise_signal.emit((channel_count,n_snapshots),1.0)
    Y += N
    return Y

Once Y is returned, I do not tamper with it in any other way.
The full music function is as follows:
def music(A,x,k):
    R = (x @ x.conj().T) / x.shape[1]
    Ev, E = np.linalg.eigh(R)
    En = E[:,:-k]
    v = En.T.conj() @ A
    return np.reciprocal(np.sum(v * v.conj(), axis=0).real),En

Steering matrix:
def steering_matrix(sensors,sources,wavelength):
    sources = np.deg2rad(sources)

    s = 2 * np.pi / wavelength
    cos_el = np.cos(sources[:, 1])
    cc = cos_el * np.cos(sources[:, 0])
    cs = cos_el * np.sin(sources[:, 0])

    D = s * (np.outer(sensors[:, 0], cc) +
            np.outer(sensors[:, 1], cs) +
            np.outer(sensors[:, 2], np.sin(sources[:, 1])))

    A = np.exp(1j * D)

    return A

I am nearly certain that this is related to complex numbers. If I synthesize a signal using the steering matrix, I get a signal that has mixed values on both the real and the imag portions, but my real world data has no imag portion. So if the implementation I am using for MUSIC requires some of the information to be retained on the imag portion, everything would mess up. I am specifically questioning my use of conj().


texts: How to perform spectral inversion in the frequency domain to convert a low-pass filter into a high-pass filter?

texts:
To convert a linear-phase FIR low-pass filter into a high-pass filter with the same cut-off frequency, we can invert the sign of the low-pass filter's impulse response $h(n)$ and then add one to the center point: $\delta(n)-h(n)$, where the $\delta(n)$ is the unit impulse function.

I am trying to verify the result in the frequency domain. If the DFT of $h(n)$ is $H(f)$, then the DFT of $\delta(n)-h(n)$ is $1-H(f)$. However, in MATLAB it seems that I cannot generate the desired high-pass filter kernel by applying the ifft() function to $1-H(f)$. The unit impulse function is added to the first point of the high-pass filter kernel instead of the center point. What is the reason for this and how can I correct my code?

MATLAB code to generate this figure:

%% Generate window-sinc filter coefficients
% CutOffFreq and TransBandWidth should be expressed as fractions of sampling rate between 0 and 0.5.
CutOffFreq = 0.1;
TransBandWidth = 0.01;

% Approxiately calculate the filter order (must be a even integer).
M = 4 / TransBandWidth;
M = round((M-2)/2)*2+2; % Round to the nearest even integer

% Generate the sinc function. Filter length = filter order + 1.
idx = 0:1:M;
Sinc = sin(2*pi*CutOffFreq*(idx-M/2)) ./ (idx-M/2);
Sinc(M/2+1) = 2 * pi * CutOffFreq;  % Solve the divided by zero problem

% Generate the Blackman window.
BlackMan = 0.42 - 0.5*cos(2*pi*idx/M) + 0.08*cos(4*pi*idx/M);

% Generate the filter coefficients (the impulse response function).
b = Sinc .* BlackMan;

% Normalized the filter coefficients to have unity gain at DC.
b = b / sum(b);

%% Spectral inversion in the frequency domain
bfft = fft(b);
b1 = ifft(1-bfft);

figure('Color','w');
plot(b,'Color','b','LineWidth',1);
hold on;
plot(b1,'Color','m','LineWidth',1);
xlim([-1,length(b)+1]);
legend({'Low-pass filter kernel','Wrong high-pass kernel'});


texts: How can I synchronize signal from an LED and TTL pulses better?

texts:
I would like to synchronize my video and voltage data using Python. (I obtain these data using Bonsai and Open Ephys for an open field exploration experiment.)

I sent $150\text{ ms}$ long TTL pulses to a computer recording the voltage data, at the same time as to an LED in the field of view of my camera to turn it on. I sent the pulses at random times with 20-60 second gaps. I would like to align the LED pulses with the TTLs to synchronize my datasets.

The two recordings did not necessarily start and end at the same time, so it was possible for one of the datasets to have fewer pulses than the other one.

The voltage data was recorded at $30000\text{ Hz}$ sampling rate, and the video at about $30\text{ Hz}$.

I tried to downsample the voltage data, correlate the two arrays and find the maximum correlation to calculate the lag and then shift the video data by the offset.

My first problem is downsampling the voltage data. I calculated the real sampling rates from the average step size, and found that the camera is about $30.08207\text{ Hz}$.

I tried to use resampy's resample:

resampy.resample(voltage_data, voltage_data_sampling_rate, video_sampling_rate)

where voltage_data_sampling rate was $30000\text{ Hz}$ and video_sampling_rate $30.08207\text{ Hz}$. This returns an array that is the same size as the original voltage array.

1. How could I get a signal from this that has the same sampling
    rate as my video?

To resample my data without using this I generated an array with multiples of the rate of sampling rates ($\frac{30000}{30.082074} = 997.271$)), converted these to integers and took these indices. I understand that this is very basic, but I thought it would allow me to check if I can write the rest of this code.

I calculated the correlation using numpy.correlate:

correlation = np.correlate(video, voltage, "full")

I set a threshold based on standard deviation and the median and replaced all values with $0s$ below threshold to remove noise.

I checked which array was shorter and padded the end with $0s$.

short_array = np.pad(short_array, (0, len(long_array)-len(short_array)), 'constant')

I calculated the lag between the two signals based on the maximum correlation:

lag = (np.argmax(corr) - (corr.size + 1)/2)/avg_sampling_rate_video

This result is about right. It reduces my 2-5 second lags to about $100\text{ ms}$.

plt.plot(downsampled_voltage_time, downsampled_voltage_TTL)
plt.plot(video_time, video_sync_pulses)

I need this to be better than that, I would like the lag to be near $30\text{ ms}$.

2. Is the downsampling causing this problem, or are there other problems with this approach?


texts:
I am looking for a framework for extracting hand region from the video feed from the webcam(analogous to Viola Jones HAAR detection of face).
This region of interest will serve as input for a CNN for hand gesture recognition, which I have already coded using python and CAFFE.
There are way too many frameworks and I am confused on which one to use,
priority would be to minimise the processing time. I already tried HAAR cascades for hand detection but they are pretty unreliable.
It would be most helpful if implementation is readily available in python.

Thanks in advance


texts: Calculate Horizontal Pixel Shift for Super Resolution

texts: I have a sample that is moving at a constant rate which is known.  Since I know the distance traveled I can equate that with the number of pixels for the pixels shift.  The only way I can think to do this is to do a camera calibration but I am working on a small scale and was wondering if there is another way to calculate the pixel shift for each image?


texts: Why does this video's audio change depending on the playing device?

texts:
Yesterday I noticed that if I played this video from my smartphone on speaker, one could barely hear anything. On the other hand, if I used earphones I could hear it perfectly.

Then, I tried with several devices (other phones, my PC's speakers, a Bluetooth speaker) and on some of them one could listen to the audio with a really high quality, but on others it was very low and like 'echoish'.

I believe there should be a signal processing explanation to this phenomenon. I don't really understand what's going on, but I suspect it has something to do with stereo and mono audio. Sound quality was good on those devices that were stereo (speakers from PC, earphones, headsets) but on those that were mono (speaker of a cellphone, Bluetooth speaker) I could barely hear anything.

What is really going on? Is this phenomenon a mistake from the video uploader or could it be an intended effect?


texts: What does it mean for a function to have frequencies?

texts: In a lecture, my professor mentioned that $\cos$ has two frequencies. I see that using the inverse Euler's formula we can express $\cos$ as a some linear combination of complex exponentials, each with a $\omega$ term.  But, I don't really have an intuitive understanding of what it means for a function to "have" a frequency. What other functions have a frequency, and why?


texts: Dominant Frequency Peak Decreases with Increasing Window Size

texts:
I have a signal that looks like this.

I analyse it using fast Fourier transforms to identify the frequency with the largest peak, which is always close to zero. (There are no other clear peaks.) If I use windows of different sizes (by "window" I mean size of the snippet I use – I don't know if this is the correct use of the terminology, however), the result shows that the dominant peak gets closer and closer to zero as the size of the windows increase. That is to say, if I plot 1/frequency against window size, it is linear.

My initial thought is that if I had an infinitely long window, then the peak would be located at zero implying a DC offset, but the signal is centred on zero in an effort to avoid such an offset. Can someone explain why I see this phenomenon and how I can avoid it?

PS I even tried this using sliding windows so that I can get frequencies based on different parts of the data, but the issue persists.


texts: detecting edge direction using Julia-language

texts:
Let say we have following image:

made with following Julia-code:

using Images
using PyPlot
using ImageFiltering

function chessboardimg(m::Int, n::Int, blkwidth::Int, blkheight::Int, )
  img = zeros(m, n)

  togglex = true
  toggley = !togglex
  for i=1:m
    togglex = true
    for j=1:n
      img[i, j] = abs(togglex - toggley)
      if j % blkwidth <= 0
        togglex = !togglex
      end
    end
    if i % blkheight <= 0
      toggley = !toggley
    end
  end
  return img
end

### creating image ####
img1 = transpose(triu(ones(20, 20)))
img2 = chessboardimg(20, 20, 10, 10)
img = [img1 img2]

Now i want to get the edge directions in this image as below:

gx, gy = imgradients(img)
mag, grad_angle = magnitude_phase(gx, gy)
grad_angle *= 180 / pi
imshow(grad_angle, "gray", interpolation="none")

These are the images of each calculated step:

as you can see in the "Gradient Direction"-image twoh vertical lines are missing.
which should actually indicate the transition from black to white region in the y-direction.

1) what am i doing wrong?


texts: Moving average filter output magnitude in simulink

texts:
I was using the Moving average filter provided by simulink. I set the window length equal to 31 samples and i was using a fixed step solver with a step size of
$\frac{1}{(50\times24)}$. I used a unit amplitude sine signal with a frequency of 50 Hz as the input to this filter (This signal is sampled at the same rate as the solver step size i.e. $1/(50\times 24)$ or $ \sin(\omega_{0} \frac{n}{(50 \times 24)})$. The filter output magnitude was around (approximately) 0.2 p.u. Now, by changing my solver step size to $\frac{1}{(50\times24\times31)}$ (which again is the same as my signal sampling rate i.e $ \sin(\omega_{0} \frac{n}{(50 \times 24 \times 31)})$),  where 31 pertains to the filter window length, the resulting filtered signal had the correct amplitude (approximately 0.9971) 1 p.u.
My questions are as follows, 1) why is the amplitude very low in the filter output when i do not incorporate the filter length in the step size of the solver. 2) How do i correct for this amplitude decrease without changing my solver step size i.e. in the first case.


texts: Training custom classifiers using trainCascadeObjectDetector

texts:
Is there any good documentation about training your own classifiers(with examples included) different from this one provided by MathWorks.
Thanks in Advance


texts:
I have a digital audio signal (float samples). I need to apply the Vibrato effect to a signal. "Vibrato" is no common term, but it is the same as the frequency modulation.
So, I need to modulate my audio signal by a low-frequency sinusoid. I need to set two parameters for the modulation: "amplitude" and period. In this case, "amplitude" of freq modulation is a maximum shift of signal frequency, for example 10 Hz. Period is just a period of the sinusoid.
So, I need a C/C++ function for such modulation.
The related article, but I cannot quite understand it: https://christianfloisand.wordpress.com/2012/05/01/shaking-it-up-with-vibrato-2/


texts: Upsampling and downsampling signals as a preprocessing step for a neural network

texts:
I have audio data acquired from a 4 channels sensors array.

As a preprocessing step for a neural network, I want to beamform and focus on the sound source. For higher resolution in the beamforming process, I have upsampled the signal with a factor of 5 using cubic interpolation, from 48 kHz to 240 kHz.

I am aware that now my network will train much slower due to the longer vectors. Nevertheless, my network is not performing well (it did without preprocessing step). Is there any connection? Is there any reason to downsample the frequency back for machine learning purposes specifically and for classical signal processing method in general?


texts: Derivation of Circular Mean Square Error

texts:
I would like to understand how Eq. (36) in [2] was derived:

The rationale behind the definition of circular sample mean in Eq. (37) is clear, but there is no motiviation for the CMSE definition in Eq. (36).

[2] Lovell, Brian C., and Robert C. Williamson. "The statistical performance of some instantaneous frequency estimators." IEEE Transactions on Signal Processing 40.7 (1992): 1708-1723.


texts:
I am stuck on this exercise. I don't know how to deal with this squared y in the denominator. What am I supposed to do to obtain state space equations?

ex.1
Simplified dynamic model of the steel ball kept in the magnetic field has the following form:
$$my'' = mg - C{u^2 \over y^2}.$$

u – voltage of the magnet, y – vertical ball position, g – gravitational constant, m – mass, c – constant. We
consider normalization m = g = c = 1.

Write down state space model equations.
What is the operational point for y = 1?
Write down the linear model equations around that point.
Write down the respective transfer function.


texts: For an LTI system, why does the Fourier transform of the impulse response give the frequency response?

texts:
I know that for a given system, the Fourier transform of its impulse response gives its frequency response.
I want to find where this property comes from, but haven't been able to find if it's a definition or if there's a mathematical proof available, for both continuous and discrete-time systems.


texts:
The aim is to identify problems in an ADC FFT simulation in an early phase of the simulation using the output bitstream of the ADC. The reason is to avoid long simulation times while using wrong data.
Please if you are aware of papers, documentation or books on the topic write it in the answers.


texts: Where does the imaginary part of a singnal come from?

texts:
So I recently watched this lovely overview of DSP/SDR explaining how demodulation of AM and FM becomes trivial by turning a signal into a complex number. (AM being the magnitude of the number and FM being the phase angle difference between two samples).
But I'm still confused. It's simply transferring the problem to "How do you get the envelope of a signal" to "how do you represent the signal in 2d". Where does this imaginary component come from?
The guy says to get the envelope of an AM waveform you need only one sample. And while you already are in the complex plane, yes. Sure. Get the magnitude of that one complex number. But a sample, as far as I understood, is a real value. For example a power measurement from a sine wave generator. So you get "0.7V". What's the envelope? It could be anything. Could be a rising edge of the sine, falling edge, a peak, even a valley if you have DC offset. So I call BS you need only one sample.
The person in the video simply glosses over how you actually get the imaginary part? Do you just add a sine wave imaginary component? (sin(t)*i)? Do you phase shift your entire signal by 90° (f(t-tau/4) * i) and add that? Or something else entirely?
Heck I don't even know how to tag this question.


texts: Seeing spurious peaks in acquired vibration data; peak frequency changes with sampling rate

texts:
I am acquiring vibration data from an accelerometer using a DAQ. Due to some reason I am seeing spurious peaks/impacts in the FFT spectra at specific frequency and its harmonics all the way till the Nyquist frequency. This frequency changes with the sampling rate as follows:

Sampling rate = 48 KHz, peaks at 100 Hz and its harmonics
Sampling rate = 8 KHz, peaks at 83.33 Hz and its harmonics
Sampling rate = 24 KHz, peaks at 93.75 Hz and its harmonics

These peaks don't seem to correspond to any real mechanical/vibration source. Initially I suspected aliasing, but it seems the DAQ has anti-aliasing built-in. I also tried doing anti-aliasing in software using python/scipy.signal.decimate(), but it didn't help.
I am attaching below couple of FFT spectra plots for the 24 KHz case.
Can you please advise what could be possible reason for such spurious FFT peaks?


texts: How Is Mixed Norm ($ {L}_{1, 2 }$) Better than $ {L}_{1} $ Norm for Sparse Representation?

texts:
Using $ {l}_{1} $-norm regularization for the purpose of achieving sparsity of the solution has been successfully applied a lot. But many times, I found the paper using mixed-norm instead of $l_1$-norm.

Considering the mixed norm $ {L}_{p,q}$ norm defined as:

$$ {\left\| a \right\|}_{p,q} = \left( \sum_{i} \left( \sum_{j} {\left| {a}_{i,j} \right|}^{p} \right)^{
\frac{q}{p} } \right)^{ \frac{1}{q} } $$

How is mixed-norm better than $l1$-norm for sparse representation? Most of what I have seen is $l_{1,2}$.

Reference, for example,

Sparse Nonnegative Matrix Factorization for Clustering eq.12
Compressive Sensing of Jointly Sparse Signals as a Method for
Dimensionality reduction of Mass Spectrometry Data eq.7
L1/2 Sparsity Constrained Nonnegative Matrix Factorization for Hyperspectral Unmixing


texts: Are time series data always contain noise?

texts: I wonder if every time series data should contain noise or not. For example I am taking the price of a ticker, say Yahoo, every hour and noting the values. Does this data contains noise or not?


texts: Relationship between Frequency response, magnitude spectrum and gain

texts:
I have been reading quite a lot of DSP related materials online and offline (books etc.) lately. I come across the terms frequency response and gain. I have also noticed that both are expressed in dB.

I do understand the fact that expressing the magnitude spectrum in dB is more informative and intuitive and helps identify the presence of different frequencies. But as seen here (or here) and here, the plots basically have the same/identical parameters on X-axis as well as Y-axis i.e, its dB vs frequency. This makes me a bit confused.

What is the relationship between frequency response and gain. Are they the same thing? Are they used according to the context (in that case on what occasions it is apt to use either of them?)?


texts: Can someone explain me the statistical relationship of random signals with a linear time-invariant and deterministic channel to the output of the channel when having random signal as an input?


texts: Calculate lines of the Inverse of Autocorrelation Matrix

texts:
I need to calculate the inverse of a autocorrelation matrix

$$\mathbf R_{xx} = E\left\{\mathbf x \mathbf x^T\right\}$$

Where the samples $\mathbf x$ are $266000\times 1$ vectors, which means I'll have a $266000\times 266000$ matrix and need to invert it. But I only need the first 729 lines. The problem is, I have to find a classifier $\mathbf h$, where
$$\mathbf h = \mathbf R_{xx}^{-1}\mathbf u$$

where $\mathbf R_{xx}$ is the autocorrelation matrix and $\mathbf u$ is a $266000\times 1$ vector.

I have $\mathbf u$ and I can calculate $\mathbf R_{xx}$, but I can't save it on memory (it's a matrix larger than 400GB and I only have 380GB available). I'm only interested in the first 729 lines of $\mathbf h$


texts: Higher-order moment of output of LTI system

texts:
Assume a very simple LTI system. Assume $x$ is white Gaussian i.i.d. with variance $\sigma^2$.

The output variance is straightforward to obtain. For example, for a continuous-time system:

$$\mbox{Var}(y) = \mathbb{E}(y^2) = \sigma^2 \int_0^{\infty}|H(s)|^2 \,\mathrm d s$$

However, I am interested in higher order moments, for example, $\mathbb{E}(y^4)$. How can this be obtained (is there a closed form solution in the first place)?

PS: Note that $\mathbb{E}(y^4)\neq\sigma^4 \int_0^{\infty}|H(s)|^4 \, \mathrm d s$.


texts: When to use the DTFT vs the DFT (and their inverses) in analysis?

texts:
In many of my readings, whenever some author speaks about working in the frequency (transform) domain (of a digital signal), they often times take the DFT, or the DTFT, (and of course their corresponding inverses). Different authors will tend to work with one or the other.

I have not been able to really ascertain a particular pattern regarding this. In that, why would you pick the DTFT over the DFT or vice versa in explaining algorithms? Where does one help you over the other?


texts:
I am recording a 17khz signal at a sample freq of 44.1khz. I want to perform cross-correlation between the received and transmitted signals for calculating TDOA. But when I do xcorr, the results are inconsistent. Is it because my sample freq is too low for my signal freq? If so, how can I overcome this? Will upsampling help?

Edit:
My received signal also consists of noise. I am passing my received signal through a bandpass filter with limits 16khz to 18khz. I am using butterworth filter. Is it better to use FIR filter instead?


texts:
We currently have C/C++ FFT code with Hann Window function support. For a large data (over 10,000) we want to split this into 1024 samples with 50% overlap. For each 1024 sample, we window and fft.

What is the recommended way to combine all the results from these smaller overlapping results to obtain a presentable results for the over 10,000 data?

Update: For the poor explanation.
The process is best illustrated by the following links

https://spectraplus.com/DT_help/overlap_percentage.htm

https://www.youtube.com/watch?v=LlyH6YciDhw


texts: How do I manually plot the frequency response of a bandpass Butterworth filter in MATLAB without freqz function?

texts:
I have code like below that applies a bandpass filter onto a signal. I am quite a noob at DSP and I want to understand what is going on behind the scenes before I proceed.

To do this, I want to know how to plot the frequency response of the filter without using freqz.

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)

Given the outputs [b, a] how would I do this? This seems like it would be a simple task, but I'm having a hard time finding what I need in the documentation or online.

I'd also like to understand how to do this as quickly as possible, e.g. using an fft or other fast algorithm.


texts: Entropy when having symbols with same probability

texts:
I am trying to work out this problem:

Assume that we are transmitting data from a dictionary of 15 symbols. Let the
probability of sending one symbol be the same as that of sending any other
symbol. Calculate how often we need to fetch a symbol to transmit information at
a rate of 5 bits/s.

How can this be found?


texts: What is the state of the art for vocal elimination/reduction?

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

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


texts:
What specifically is PCM? I've seen it refer to:

the sampling and quantisation of an analogue domain signal into digital samples

those sample quantisations alone in the digital domain i.e. in memory or storage i.e. .pcm files

the transmission of quantisations or indeed any bitstream as a digital signal with a digital line code through a DAC into the analogue domain

the process of 1+2

the process of 1+3


texts: Playing samples in a loop, without clicks or pops

texts:
I'm trying to build an audio player that has this specific feature: when you scroll to x sample time, [x to x + 4096] samples are played in a loop until you move the mouse to another place.

The samples are taken from the backing audio. The problem is, I fail to play them back as a continuous note: the looping is not seamless.

So far I tried to start and end the loop at zero crossings. This slightly improved it, but I can still hear the seamless looping.

One idea I had was to go to the frequency domain, remove all, but the max magnitude bin and do an inverse FFT afterwards. This, sort of works, the looping is better, not perfect, but better, however, removing all the other freq components from the loop leaves it really "dry", not resembling the original (missing the harmonics)

Is there a way to handle this in the time domain? Where should I look? Are there any filters/algorithms that could help me?


texts: 2D Digital Signal Processing - Applying High Pass on and Image

texts:
I have been doing my research on 2D digital signal processing. What I am doing is about noise removal of digital image:

Type of noise that I want to remove from an images is like the Gaussian noise and 'salt and pepper'. For your information, I intend to use the butterworth lowpass filter to remove the noise. This is my code that I used in matlab.

Function Highpass

function [ out ] = butterhp( im,d,n )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
h=size(im,1);
w=size(im,2);
[x,y]=meshgrid(-floor(w/2):floor(w-1)/2); -floor(h/2):floor(h-1)/2;
out= 1./(1.+(d./(x.^2+y.^2).^0.5).^(2*n));
 end

Function Lowpass

function [ out] = butterlp( im,d,n)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
out=1-butterhp(im,d,n);
end

matlab code

a = imread('penang.jpg')
f = rgb2gray(a);
k = imnoise(f,'gaussian')
lp=butterlp(k,15,1);
af=fftshift(fft2(k));
aflp=af.*lp;
aflpi=ifft2(aflp)
noise = uint8(real(aflpi));
imshow(noise)

This is the image that i used.

And running the code, this is what I get.

So,my problem is when I try to filter it to make it looks more clear by adjusting the number of 'n-th order' and 'cutoff frequency', it seems it does not work. It still gives me blur picture just like that.


texts: Acceleration time series to velocity conversion using DFTs

texts:
I have 1000 Hz time series data for acceleration (512 data points), which I want to convert to velocity. I am trying to use the omega arithmetic method to achieve this. Following are the steps I am taking:

Converting the acceleration-time data to its DFT.
Converting the acceleration DFT to velocity DFT using the omega arithmetic formula(dividing each DFT value by i * omega (i.e. 2 * pi * i * freq))
Taking the inverse FFT of the velocity DFT to get velocity-time data.

These steps can be seen in the picture above. I am using numpy's FFT functions, and not normalizing the DFT to avoid any confusion.
My question is: Why are the velocity values so low in the velocity-time plot? Am I making any mistake in my steps?


texts: Inverse $\mathcal Z$-transform when region of convergence goes outwards from the inner pole?

texts:
I am looking for the inverse $\mathcal Z$-transform of the following:

$$
\frac{1}{1-\frac 12 z^{-1}}+\frac{1}{1+\frac 13 z^{-1}}
$$

When the region of convergence is $z > 1/3$. I have found the $\mathcal Z$-transform for when $z > 1/2$ as a right sided signal:

$$
\left(\frac{1}{2}\right)^n u[n] + \left(-\frac{1}{3}\right)^n u[n]
$$

I can't seem to find a inverse $\mathcal Z$-transform when the ROC goes outwards from the inner pole. Does the inverse $\mathcal Z$-transform not exist for $z>1/3$?


texts: Emulating tape distortion of low-level audio signals

texts:
I'm working on a tape echo emulation in Reaktor, and I've based the frequency response, tape speed warble, high-amplitude saturation and tape hiss on a paper I found about the Roland Space Echo. However, I'm now looking for information about the distortion of low-amplitude signals by audio tape. When I let an echo trail die out to just above the tape hiss level, and then turn up the feedback level to make it audible again, the sound is too clean, whereas a real tape delay would produce a heavily distorted sound.
Which simple curve or mathematical transformation of a signal would emulate the distortion caused by recording audio on tape at very low levels? And how should I adapt it to emulate different tape speeds?
It doesn't have to be very realistic or high quality, just something that will add some tape-like dirt to low-level signals only, and be mostly transparent at higher levels.


texts:
I am not sure if this question will be a little off-topic on this forum, that I will give it a try anyway, since it implies signal process arguments.

By using MATLAB, I'm on my way to generate a turbulent wind field. This wind field is generally a Gaussian process with all the time series having 0 mean value and unitary standard deviation. I will shortly introduce you to the problem:

1) Discretize the space in Y-Z plane in a certain number of points Ny, Nz;

2) Discretize the time scale (generally 600s) with a proper time sample, dt = 0.15;

3) Calculate distance among points;

4) Generate only positive frequency vector f = (0:length(t)/2-1)*df;

5) Generate power spectra

$$ S_{ii}(f)=\frac{4\sigma_{ii}^{2}\left(L_{ii}/U\right)}{\left(1+6f\left(L_{ii}/U\right)\right)^{5/3}} $$
where
$$ ii = u,v,w\,\, (speed\,\,components\,in\,x\,y\,z\,direction)$$
$$ U = const $$
$$ L_{ii} = const $$

6) Generate coherence matix

$$ Coh_{u}=\exp\left(-12\sqrt{\left(\frac{f*dist}{U}\right)^{2}+\left(0.12\frac{dist}{U}\right)^{2}}\right) $$
$$ Coh_{v,w}= \exp\left(-12\, dist\,\frac{f}{U}\right) $$

7) Generate random numbers with 0 mean and 1 std

nn_u = complex(normrnd(0,1,size(H_u)),normrnd(0,1,size(H_u)));
nn_v = complex(normrnd(0,1,size(H_v)),normrnd(0,1,size(H_v)));
nn_w = complex(normrnd(0,1,size(H_w)),normrnd(0,1,size(H_w)));

8) Cholesky factorization of the Coherence matrix (retrieves lower triangular factor)

9) Generate time series

$$ U = [real(fft(Coh_{u}*nn_{u}))\,\,imag(fft(Coh_{u}*nn_{u}))] $$
$$ V = [real(fft(Coh_{v}*nn_{v}))\,\,imag(fft(Coh_{v}*nn_{v}))] $$
$$ W = [real(fft(Coh_{w}*nn_{w}))\,\,imag(fft(Coh_{w}*nn_{w}))] $$

Now, after this consistent theorethical introduction, I can move towards the real matter: when calculating the standard deviations for eac time series (u,v,w components), they don't match the expected values
$$ \sigma_{u} = 1 $$
$$ \sigma_{v} = 0.8 $$
$$ \sigma_{w} = 0.5 $$
because of several reasons:

space discretization introducing smoothing effects
time sample period
power spectrum only within a short frequency range

Btw, I would like to, let's say, scale the time series in order they to be still Gaussian variables with zero mean and a standard deviation better approaching the expected values.

Do you have any hint?


texts: What can be done with I/Q data which would otherwise be impossible?

texts:
As I understand it, I/Q data is the complex representation of any waveform. It enables us to understand two things at any arbitrary point in time:

The projected amplitude of the waveform.
The direction of the waveform.

However, it's still unclear to me what makes I/Q data so special.
Logically, I can discern the direction of a waveform implicitly: If I'm talking (or transmitting), I know a signal is travelling away from me. On the other hand, if I'm listening, I know any signal I hear must have been travelling toward me.
And if I want to know the amplitude of a waveform, I can just wait until I see a peak.
This in mind, what sorts of modulations/demodulations would be impossible with real data instead of I/Q data? Respectfully, I'd appreciate an example based on reality as I understand that I/Q data makes the result of a FFT transform deterministic, but I don't really understand why that's relevant.


texts: What is the phase and magnitude response of white noise?

texts:
I would like to create white noise in the frequency domain, and then transform it to the time domain using python. To understand the problem, I simply generated white noise in the time domain, and transformed it to freq domain:

import scipy.signal as sg
import numpy as np
import matplotlib.pyplot as plt

e = np.random.normal(0,1,1e3)
E = sg.fft(e)

plt.figure("Bode plot")
plt.subplot(211)
plt.title("Magitude")
plt.plot(abs(E))
plt.subplot(212)
plt.title("Phase")
plt.plot(np.angle(E))
plt.show()

I does not look at all as I expected:

Questions:

Isn't white noise supposed to have a flat magnitude response? (equal amounts for all frequencies)
What is the relationship between the standard deviation (1 in my example) and the magnitude and phase?

Thank you in advance!


texts: Why a Large Diaphragm Microphone is limited to lesser frequency range

texts:
Here Above are the frequency responses of different size microphone(condenser)

what actually is the relation between the microphone capsule size and frequency response ?

  The conclusion is that a large diaphragm microphone will have a more limited frequency range than a small diaphragm. This is illustrated in the graphs

i found some sites explaining the reason but i did not get the reason well, can any one illustrate for me plz


texts:
I have this block-diagram of Karplus-Strong basic algorithm.

Where $$H_{d}(z) = \frac{1+z^{-1}}{2}$$
The problem occurs when I try to get the difference equation. As I can see from the diagram, that would be
$$y(n) = x(n-N)+\frac{y(n)+y(n-1)}{2}$$
but I've found online this one(which I know how to implement in code, unlike the first one)$$y(n) = x(n)+\frac{y(n-N)+y(n-(N+1))}{2}$$
Are these two equation somehow equivalent?


texts: FFT returns large low-frequency power - probably because signal is interpreted as containing cycles that do not fit. How is this called?

texts:
So I have a large number of signals like the one in the first picture below and I would like to extract and compare the frequencies within them.
I applied a Fourier transform which resulted in results like those shown in the second picture. The large power of very low frequencies suggests to me that the signals include something that is interpreted as being of such low frequency that one cycle does not fit into the time window (which is 1.2 seconds).
I hence applied a High pass filter with different cut-offs and the pattern is gone after using a HPF with 2Hz (last picture).

I have two questions about this:
1) Does my reasoning make sense?
2) How is such an effect / pattern called?

I need to include this in a paper I am currently drafting and hence am especially unsure about question 2). My supervisor's first reaction was that these are edge-effects, but looking into them I feel like that is not really what is at work here.

Thank you!!


texts: What does the operator 'Range Detector' do in this block diagram? How can it be modeled?

texts:
I'm trying to implement a 'Compressor/Expander/Noise gate' for automatic gain control of audio signals. The block diagram of the system is shown below:

I don't understand how the 'Range Detector' selects the appropriate output, given that for some values of xRMS(n), both the output of the expander and noise gate are non-zero.


texts: How to use ImageNet Database in Matlab for object classification or recognition? I didn't find any tutorial which could practically tell me how to train images on ImageNet database and classify the objects.


texts: Using Total Variation Denoising to Clean Accelerometer Data

texts:
I know this is maybe a very basic question but I am doing this as a hobby and I can't find a solution to this problem.
Basically I am trying to remove some noise from data I am reading from an accelerometer. This is what I want to achieve (Taken from Total Variation Denoising (An MM algorithm)):

I read in Picking the correct filter for accelerometer data that Total Variaton Denoising would fit my needs. So I read Wikipedia - Total Variation Denoising article from Wikipedia and I think I have to use one of this equations:

But I don't understand how I apply this to my signal. Suppose I have a set of x,y points like in the plots above, how I apply the equation to that data? I implemented some simple low-pass and high-pass filters like this:

gravity[0] = alpha * gravity[0] + (1 - alpha) * event.values[0];

But this is maybe too complex and I don't know where to start or how. I want to implement this in Java or C so Matlab is not an option (I have seen a lot of MatLab implementing this). I will appreciate any help to guide me in the right direction!


texts: I am working on calibrating my recording system. I have a signal generator playing a 500mVpp sine wave signal into my utrack24 recorder. I am trying to determine what the units of the recorded utrack .wav file is so I can compare that to the input signal. I have read the .wav file into matlab with audioread. I can use either the normalized version of the file or the 'native' version; I'm not sure which is better.


texts:
How to recognize a zebra crossing from top view using OpenCV?
It is not a straight zebra crossing - it has some twists and turns.
I have some ideas:

parallel line detection, but the available techniques only work to find straight parallel lines, not curved
template matching to match a template of back to back white and black stripes, but it becomes tedious since I can't find any pattern matching techniques with scaling and rotation


texts: Complex analysis book recommendation for electrical. eng

texts: I need recommendation about complex analysis book. As I'm electrical engineering student, it should cover everything one engineer need to know about that mathematical field, but without strict mathematical formalism.


texts: Why should one choose the maximum as reference for power-to-decibel conversion?

texts:
I have seen that one converts the amplitude to decibels through:
$\text{dB(S)}=10 \text{log}_{10} \big(\text{S/ref}\big)$, where $\text{S}$ is the output of a STFT and ref the reference value for the logarithmic scale, i.e. the value that should be mapped to decibel, right?

However, I've seen implementations with $\text{ref = max(S)}$
and I do not understand why. Can someone clarify on that?


texts: How to deal with a negative pole (unstable) in the pre-filter of a control system?

texts:
So while answering how to design a PI controller for a first order time delayed system   (Question Here )

Here is the closed loop equation to a control system:

$$
G_C(s) = \frac{\frac{K}{T}(1-sT)(s)} { s^3 + (\frac{1}{T} + a - KK_p)s^2 + (\frac{a}{T} + \frac{KK_P}{T} +K_I)s+\frac{KK_I}{T}}
$$

Question: How do you deal with normalizing the numerator in your closed loop transfer function when the filter is unstable? (Pole on RH of plane)

Typically you introduce a filter before your controller that does:

$$ \frac{1} {\frac{K}{T} (1-sT)(s)} $$

to normalize the numerator

But the filter itself is unstable because of the term:

$$ \frac{1}{(1-sT)}$$ is unstable for a step response which would create an issue realizing the system at all.

One way I've thought about dealing with this is multiplying it by its complex conjugate
$$
\frac{(1+sT)} {(1+sT)}$$

but im not really sure about the merits of it.


texts:
The goal is to create an LTI filter which is exactly, or approximates, damping of harmonic modes.

The equation of course is:

$$\frac{d^2 x}{dt^2} + 2 \xi \omega \frac{dx}{dt}+\omega^2x=0$$

This can be done with very costly convolution, I suppose.

I was thinking that instead of treating the velocity term as a part of a time dependent potential (no longer time invariant), I can treat this as a two dimentional state and therefore eliminate second derivatives and have it a $z^{-1}$ feedback filter. It is working weird.

A continous alternative of something close is:

$$e^{-\xi t} \cos \omega t$$

How is this approximated efficiently?

Links are enough.


texts: Determine filter type using recurrence relation

texts:
Given the recurrence relation:
$y[n] = x[n] + 0.5y[n-1]$
I want to determine the filter type (i.e. LPF, HPF etc.)

I try to use Z transform, and get that the the transfer function is
$H(z) = \frac{2z}{2z-1}$
So I derive that the filter is HPF (High pass). Is it correct?


texts: Why do we get different imaginary parts of a zero centered Gaussian for the the same number of data points N?

texts:
Suppose we have a total number N= 2048 points in a data and we wish to have zero centered Gaussian. There are two possibilities that we use the x-axis as

x1=[-1023:1:1024]; % x axis spans from -1023 to 1024 with 1 unit steps

x2=[-1024:1:1023]; % x axis spans from -1024 to 1023 with 1 unit steps

and if we make two zero centered Gaussians using these x values:

Gauss1=normpdf(x1,0, 10); % The syntax is normpdf(x, mean, standard deviation)

Gauss2=normpdf(x2,0,10);

and obtain their FTs as follows in MATLAB.

FFTGauss1=fft(Gauss1);
FFTGauss2=fft(Gauss2);

The real parts are identical and their magnitudes exactly match. For some reasons, the imaginary parts vary drastically. Why do we see large imaginary parts in one case and almost non-existent imaginary parts in the other? Thanks.


texts:
Learning the cepstrum analysis for speech recognition, I have met two different definitions of cepstrum (for discrete signals):

$F^{-1}(ln|F(x[n])|)$ . That is, the cepstrum is the inverse Fourier transform of the logarithm of the magnitude of the Fourier transform of the signal $x[n]$. So, cepstrum coefficients have the dimension of time.
$F(ln|F(x[n])|)$. That is, the cepstrum coefficients are the Fourier transform of the logarithm of the magnitude of the Fourier transform of the signal $x[n]$. Hence, the spectrum of the spectrum. From here, the name of quefrency to indicate the variable of the cepstrum coefficients.

The first equation seems the most used, and intuitively makes more sense to me w.r.t the source separation for speech signal processing. However, I met the second definition quite few times, causing me a great confusion.

How the conceptually differ these two definitions? And which is the most correct one? Always in regards to speech signal processing.


texts: Conversion from laplace transform to z-transform

texts:
I would like to know if

$$
  \text {Z-Transform ( }  G(s)H(s) \text{ )} = \text {Z-Transform (}G(s) \text{)} \text { Z-Transform (} H(s) \text{) } = G(z)H(z)
$$

where G(s), H(s) are the Laplace transform representations of g and h, and G(z) and H(z) are the Z-transform representation of g and h.

Is this relation true?


texts: What exactly is meant by "stochastic sampling" and is it profoundly different from the regular Nyquist-Shannon sampling theorem? Is it related to sampling a stochastic process?


texts: What is different between Euclidean distance and the Geodesic distance?

texts:
We have two kinds of distance in image: Euclidean distance and the Geodesic distance. What is different between them? Could you show mathematic formula or visualization to make it clear? The reference said

    Difference between the Euclidean distance and the Geodesic distance calculated with the Fast Marching Method. The Geodesic distance is the distance of the minimum length inside the figure path and the Euclidean distance is the straight line distance

This is example


texts: Detecting changes between voice and music

texts:
I have a bunch of audio files that have both speech and music parts (think radio). For instance, a given file might consist of the following sections:

music --- speech --- music --- speech --- music --- speech --- music

My goal is to detect when the file changes between music and speech, and cut away the music parts leaving only the speech. However, I don't know what method to use to distinguish between speech and music. I was considering detecting the tempo, but some of the music is very slow or has a varying tempo.

(While we're at it, so that I don't have to write my own, are there any existing software libraries out there that can distinguish between speech and music?)


texts:
I am looking to design a set of FIR filters to implement a low pass filter. I am also trying to reduce the latency of the signal through the filter so I am wondering what the minimum number of taps I can use might be.

I know that more taps can lead to a sharper cutoff of the frequency and better stop band rejection etc. However what I'm interested in is more fundamental - if I want to implement a low pass filter with cutoff at $\frac{f_s}{100}$ say does that mean that I need at least 100 taps in order to attenuate the lower frequency signals? Or can I get away with less taps and if so is there some theoretical lower limit?


texts: importance of using histogram equalization

texts:
I'm reading opencv tutorials and I'm diving into histogram equalization. i have looked in wikipedia, there is a nice example that sums up exactly the problem:

original:

equalized:

but for getting this result i would take a different approach:

find the minimum and maximum in the original
normalize (remap) everything upon it.

no histogram, no  cumulative distribution function. for sure a more stupid approach but i can't see the difference.. why use histogram equalization? can someone can help me in make the reason out?


texts:
I am not sure I understand what the Mel Scale is.

Googling doesn't give me various answers. I seem to be getting the same response again and again.

Which would be something like:

"The mel scale reflects how people hear musical tone"

First of all, does the concert pitch has a specific frequency? So why do we even have a scale stating how do people hear it?

I am currently trying to understand how MFCC uses them, and it seems like it maps the frequency to the Mel scale, or there are some filter banks which are applied on the given frames


texts: Design Hilbert filter at matlab (pure 90 phase and magnitude 0db)

texts:
I need to design a filter  able to  take a measurement of alternate current and shift the phase of the input by 90°.

The filter should work between 40 and 60Hz and at this range should have unity gain. Outside this range I don`t care about the gain or phase.

I research about Hilbert filter and sounds that this kind of filter is what I need, but I don`t know how to use.

Let me paste my matlab code (fixed afther Matt recomendation):

%     close all
clear all
clc

%Period
Fs = 1000; %sample frequency

fl = 40;   %lower frequency

f_min=fl/(Fs/2);

b = firpm(40,[f_min (1-f_min)],[1 1],'h');   % Bandpass Hilbert
%     fvtool(b)

t = 0:1/Fs:1.2;
t_window=0.2;

%input signal Frequency

f=40;
%Generate AC signal (input)- at 40hz (lower limit)
y1 = sin(2*pi*f*t);

%Generate filtered signal-40hz
y1f=filter(b,1,y1);
y1 = [zeros(1,(length(b)-1)/2), y1(1:end-(length(b)-1)/2)];

f=60;
%Generate AC signal (input)- at 60Hz (upper limit)
y2 = sin(2*pi*f*t);

%Generate filtered signal-60hz
y2f=filter(b,1,y2);
y2 = [zeros(1,(length(b)-1)/2), y2(1:end-(length(b)-1)/2)];

figure
hold on

plot(t,y1,'*-')
plot(t,y1f)
plot(t,y2,'*-')
plot(t,y2f)

hold off
grid

legend('input','filtered')
axis([max(t)-t_window max(t) -1.1 1.1])

range(y1f(end-(length(t)/10):end))/range(y2f(end-(length(t)/10):end))

I was expecting that my signal yf is 90° shifted, but this doesn`t work.

Someone have some idea about how to solve this issue and what is wrong at my code?

Part 2 - Thanks Hilmar and Richard.
Let me fix my code and use your code and show what is happening. Done. But this don`t work yet...

Objective 1- 90 phase shift -ok! (thanks guys!)
Objective 2- Magnitude without change- Not OK. The signal was atenuate by 5% if I change the input frequency from 60 to 40Hz.

Could you help me to explain why this is happen at magnitude?and What I should do to fix this?.... I need magnitude as flat as possible between this frequencies.
I tried to increase the order of the filter(40->60) but don`t solve my issue yet.

Part 3 - Thanks Matt.
Let me fix my code and use your suggestion. Now work!

Now my filtered signal is shifted properly and with low attenuation between 40 to 60hz. The attenuation is less than 0.4% at this range.

Thank again guys(Matt, Richard and Hilmar)


texts: How to find out the transfer function of a FIR filter?

texts:
$$h[n]=\begin{cases}a^n & \text{if } 0 \le n < N \\
0 & \text{otherwise}\end{cases}$$
And for which values of $a$ the filter is stable

I know that the transfer function will be

$$H(z)=\frac{z}{z-a}~,~|z|>a$$

how to find out the values of $a$ for it's stability ?


texts: gnuradio: How can I change the value of the tag exiting the "Correlate Access Code - Tag" block?

texts:
Background: I have built a conventional demodulator in gnuradio and am trying to send payload data out of the flowgraph to an external Python script. I'm using gnuradio v3.7.12.

It seems to me that the optimal method would be to use the flowgraph to:

find the preamble and tag it
convert the data following the payload to a pdu
send the pdu as a message to the external Python code

I've built a flowgraph that successfully does this using the following blocks:
input bitstream->"Correlate Access Code - Tag Stream"->"Repack Bits"->"Tagged Stream to PDU"->"ZMQ PUSH Message Sink".

The problem is that the "Correlate Access Code - Tag Stream" block expects a very specific header between the preamble and the payload: 32 bits containing the payload size (16 bits repeated twice). This is fine if I am the one controlling the input signal. Any other signal I try to receive, however, is not going to have this header field.

So instead of the "Correlate Access Code - Tag Stream" block, I tried to use the "Correlate Access Code - Tag" block, but this block sets the length tag value to 0, which prevents the rest of the signal chain from working correctly. It seems like all I need to do is change the value of the length tag to my expected payload length, but unfortunately, after many hours of scouring the documentation regarding tags, messages and pdus, I have been unable to figure out how to do this.

Am I on the right track? I've attached my working flowgraph (with the "Tag Stream" block) as well as the non-working flowgraph (with the "Tag" block). I tried a number of things (Tagged Stream Align, Stream to Tagged Stream), but I couldn't get anything to work.

Working GRC: https://www.dropbox.com/s/ilv0ry1x4nvqelj/tag_stream_test.grc?dl=0

Non-working GRC (can't set length tag!):
https://www.dropbox.com/s/aoueiovuvj9ss8z/tag_test.grc?dl=0

Thanks for any help you can provide.
PC


texts: The two types of stability and "Why exponential"

texts:
When I was learning about LTI systems,
I noticed that
LTI system is said to be BIBO stable
if and only if its square sum of impulse response is finite.
This expression is found on many textbook of signal processing.

I continued study and found this expression
LTI system is stable if and only if
its impulse response is exponentially bounded.
I could not understand this at first,
but by further research,
I elaborately acquire this interpretation.
However, I could not find explicit description about this.
So my question is "is these things true ?",

There are two types of stability (of the system), internal stability and external stability.
BIBO stability is equivalent to external stability.
If the system is controllable and observable, internal and external stability is equivalent.
If the system is linear, "exponentially bounded" is equivalent to internal stability.

and if 4. is true,
why the function must be "exponential decaying" ?
How is this restricted to be "exponential" ?


texts:
I have to do 'inverse engineering' with piano sound and find sequence of frequencies in it. I have to use spectrogram() and adjust parameters but my problem is that I don't know how. I set window length as 512 to see more frequencies and overlap on 256 (as default), is this correct? Also in my task there's: "FFT length (zero padding)" so I used nextpow2. Can somebody take a look at this?
Here's what I've already done:

[y,fs] = audioread('piano.wav');
%sound(y,fs)
f = 0:1:0.5*fs;
Nx = length(y);
NFFT=2^nextpow2(Nx);
n= 0:Nx-1;
dt = 1/fs;
t = dt*n;
figure(1)
plot(t,y)
figure(2)
spectrogram(y,512,256,NFFT,fs); %


texts: Why correlation coefficient template matching perfoms better than cross-correlation and square-differences method

texts:
I am using this template in OpenCv's match template function by converting all images to grayscale.

Here are my result using sum of square differences,cross correlation and correlation coefficient.Both SSD and Cros Correlation method fails matching but correlation coefficient performs well.

Not in this case but in my prior experience I always find cross correlation method the most robust one.I didn't understand why here correlation coefficient method works better than others ?


texts:
Last week I've been trying to implement the Order Analysis in MATLAB in vain. I've read a lot of docs about it like this one ,  but I still can't figure it out.
This will be my last shot.

Say I have recorded 10sec audio file of a rotating machine, the sampling frequency is 96Khz. In this 10sec interval, the speed of the machine will continuously increase, say from 800rpm to 1500rpm. Is there any way to calculate the Order Analysis for this machine?
I'll be grateful for any hint!


texts: Is digital up and down sampling linear/causal/time-invariant?

texts:
So I am trying to determine whether the process of up sampling (zero padding) and downsampling are linear, causal and/or time-invariant.

Based on some resources I found online, I am under the impression that it is linear, but not time-invariant. However I am unsure about causality (since it depends neither on future nor past values?).


texts: How to calculate SNR of dead-band quantizer?

texts: Suppose we have a dead-band quantizer with $\ N = 5 \ quantizing \ levels $  which uses $\ b = 3 \   bits $ for each produced symbol. The quantizer gets as input a uniformly distributed signal $\ X $ at the interval $\ [-10,10] $. How can I calculate the signal to noise ratio $\ SNR $ of the quantizer ?


texts: Detection of sine signals with random amplitudes

texts:
Problem statement:

I am designing a NP detector for the following detection problem:

$\mathcal H_{0}: x[n] = A_0\cos(2\pi f_0n) + w[n]$
$\mathcal H_{1}: x[n] = A_1\cos(2\pi f_0n) + w[n]$

where:

$f_{0}: \text{known parameter}$
$w[n] \sim \mathcal N\left(0, \sigma^2\right), \text {WGN}$
$A_{0}\ \ \sim \mathcal N\left(a_0, \sigma_0^2\right)$
$A_{1}\ \ \sim \mathcal N\left(a_{1},\sigma_1^2\right)$

What didn't work for me:

o I have tried to  directly apply Neyman-Pearson theorem in the way it's done in "Fundamentals of Statistical Signal Processing, Volume II, Detection Theory" by Steven M. Kay. But I have failed to derive test statistics for optimal detector.
o I have tried to use the amplitude of the spectral component for $f_{0}$ as a test statistics and estimate threshold of making decision from given $P_{FA}$ (probability of false alarm). It worked fine for me, but I'm still concerned that this detector is not optimal.

Question:

Is there some way to design optimal NP detector for such problem?
Any ideas, advices, links will be much appreciated.


texts:
The frequency spectrum of cosine function is as below.

While adopting a square window function to the cosine function, it is same as implementing convolution of their frequency spectrum as below. The thick dash line is the sum of them. Is it correct to say the points on this dash line are the FFT result of corresponding frequency? I feel it is not true, but don't know where the problem is.


texts:
I need to record sound using 2 (or 4, eventually) microphones and then calculate the time delay of arrival of sound at the different microphone pairs. But first, I need a microphone array. And since I'm just starting out, I would like to connect just 2 mics first and calculate the time at which sound arrives at them.
Can I do this using a sound card?
And I heard sound cards will have typically one or two input channels. What are these then? http://www.ebay.com/bhp/external-sound-card-for-laptop
The 'channels' mentioned along with the above devices are not exactly input channels?


texts: What are the state of the art methods for person segmentation and person pose identification

texts:
I know that HOG is pretty state of the art for person detection, but in its original form HOG is only a detection method, and I don't believe it is particularly fast either (there is a significant speedup which takes the sliding window and rejects unlikely possibilities early, but it has been patented).

What is the current state of the art approach for determining the pose of a person (including close poses like a head shot, and hand shot, etc...), and for segmenting the person from the background?

The setting here is still images.


texts: Breaking down a higher order filter to simple cascade of averages

texts:
I am having the following filter structure, which is a second-order decimation:

Which actually implements the following equation/code: (ans-0.2727)
$$\frac{2}{N(N+1)}\sum_{i=1}^{N} (N+1-i).d(i)$$

d=[1 0 1 0 1 0 1 0 1 0 1];
N=10;
sum2=0;
for i=1:1:N
   sum2=sum2+((N+1-i)*d(i));
end
second_order_ans=sum2/(N*(N+1));

However, I wanted to break this to two separate for loops (or) a cascade of two first-order decimations.

The first-order decimation corresponds to the following equation and code: (ans-0.5)
$$\frac{1}{N}\sum_{i=1}^{N} d(i)$$

d=[1 0 1 0 1 0 1 0 1 0 1];
N=10;
sum=0;
for i=1:1:N
sum=sum+d(i);
end
first_order_ans=sum/N;

I tried the following for cascade but I get the output as 0.6182 and far away from 0.2727.

%second order - cascade
for i=2:1:(N-2)
sum2_c1(i)=sum2_c1(i-1)+d(i); %First counting
end

for i=2:1:(N-2)
sum2_c2(i)=sum2_c1(i)+sum2_c2(i-1); %Second counting
end
sum2_c=0;

for i=2:1:(N-2)
sum2_c=sum2_c+sum2_c2(i); %adding up
end

sum2_ansc=((sum2_c)*2)/(N*(N+1));


texts: wavelet decomposition for time series signal

texts: Is it possible to use stationary wavelet decomposition as a tool to extract wavelet features for a time series? I can see how it works for image cases, but for a time series prediction problem say $x_i$, $i = 1,2,3,\cdots,n$ and you want to predict $x_{n+1}$, is it possible to use stationary wavelet decomposition to extract features and if yes why it makes sense?


texts:
I have a program which takes in data from an oscilloscope, and due to reflections in the medium, we will get random signals at random intervals. The two grey windows at the right bottom show the original data (below red line), and after FFT forward and backward (above red line).

The data set is originally 10,000 points, so the FFT has 20000 doubles stored. You can see in the command shell that the backwards FFT uses 19000 of those 20000 values. My issue is at the green rectangle I get waves before and after the real signal, which means I can hardly denoise it due to the signal itself being a rather sharp and high frequency.

Would I first find and extract the areas where the signal itself is present, and then run the FFT on that subset, or how would I go with denoising such an area?

The noise itself is not visible, since the 10k data points are scaled to 500 pixel, but its random electric noise coming from the oscilloscope.

Link to the image in big: http://i.imgur.com/3fHo8vP.png

The other two images show the similar data being captured, once being averaged 128 times, one a single set. I need to find a way to denoise the single capture.


texts: Calculate the Inverse DTFT of the DTFT Derivative in Terms of $ x \left[ n \right] $

texts:
Consider the signal $ x \left[ n \right] $ and its DTFT transform $ X \left( {e}^{j \omega} \right) $.
Assume $ X \left( {e}^{j \omega} \right) $ is differentiable.

What is the Inverse DTFT of:

$$ j \frac{d}{d \omega} X \left( {e}^{j \omega} \right) $$

Express the result in terms of $ x \left[ n \right] $.


texts: How can I realize a complex communication system in a real time application?

texts:
So I am working on a project where I have to model a communication system in MATLAB, then implement it on an FBGA.

My transmitter goes something like that:

M        = 4;   N = 1000;
Symb     = randi([0,M-1],N,1);
Mod_Symb = qammod(Symb,M);

Now, as you can see the output is complex consisting of four possibilities [+1+1i, +1-1i, -1+1i, -1-1i].  The question is (which I don't know after 4 years mobile com. school!) how do I go from those complex representations to a real time system? I mean, when I want to feed my antenna amplitude data, how do I do that?


texts: How does window length influence cut-off frequency of filter and group delay?

texts:
I'm attempting to interpret the following function. To my understanding, this aims to implement an FIR low-pass filter with a Nuttall window. Then, it filters the signal by a simple application of the convolution equivalence theorem.
static void GetFilteredSignal(int half_average_length, int fft_size,
    const fft_complex *y_spectrum, int y_length, double *filtered_signal) {
  double *low_pass_filter = new double[fft_size];
  // Nuttall window is used as a low-pass filter.
  // Cutoff frequency depends on the window length.
  NuttallWindow(half_average_length * 4, low_pass_filter);
  for (int i = half_average_length * 4; i < fft_size; ++i)
    low_pass_filter[i] = 0.0;

  fft_complex *low_pass_filter_spectrum = new fft_complex[fft_size];
  fft_plan forwardFFT = fft_plan_dft_r2c_1d(fft_size, low_pass_filter,
      low_pass_filter_spectrum, FFT_ESTIMATE);
  fft_execute(forwardFFT);

  // Convolution
  double tmp = y_spectrum[0][0] * low_pass_filter_spectrum[0][0] -
    y_spectrum[0][1] * low_pass_filter_spectrum[0][1];
  low_pass_filter_spectrum[0][1] =
    y_spectrum[0][0] * low_pass_filter_spectrum[0][1] +
    y_spectrum[0][1] * low_pass_filter_spectrum[0][0];
  low_pass_filter_spectrum[0][0] = tmp;
  for (int i = 1; i <= fft_size / 2; ++i) {
    tmp = y_spectrum[i][0] * low_pass_filter_spectrum[i][0] -
      y_spectrum[i][1] * low_pass_filter_spectrum[i][1];
    low_pass_filter_spectrum[i][1] =
      y_spectrum[i][0] * low_pass_filter_spectrum[i][1] +
      y_spectrum[i][1] * low_pass_filter_spectrum[i][0];
    low_pass_filter_spectrum[i][0] = tmp;
    low_pass_filter_spectrum[fft_size - i - 1][0] =
      low_pass_filter_spectrum[i][0];
    low_pass_filter_spectrum[fft_size - i - 1][1] =
      low_pass_filter_spectrum[i][1];
  }

  fft_plan inverseFFT = fft_plan_dft_c2r_1d(fft_size,
    low_pass_filter_spectrum, filtered_signal, FFT_ESTIMATE);
  fft_execute(inverseFFT);

  // Compensation of the delay.
  int index_bias = half_average_length * 2;
  for (int i = 0; i < y_length; ++i)
    filtered_signal[i] = filtered_signal[i + index_bias];

  fft_destroy_plan(inverseFFT);
  fft_destroy_plan(forwardFFT);
  delete[] low_pass_filter_spectrum;
  delete[] low_pass_filter;
}

There are two major points that I don't understand in this algorithm:

Half average length. What is it? It is used to calculate the filter order and the comment says it influences the cutoff frequency. I realise that the window length limits resolution when applying to a signal. My confusion is that all DSP guides I look up usually assume some filter order N as a design parameter, rather than selecting it based on a desired cut-off frequency. Is there a relationship between N and then cutoff frequency then in this case/generally?
Compensation of the delay. I'm not familiar with the details of filter design, but after reading up, I found that it is essential to shift a time-domain signal after FIR filtering. Filtering usually has frequency-dependent group delay, but FIR filters have the advantage of a constant group delay being $ \frac{N-1}{2} $, where N is the number of taps or order of the filter. Assuming $ \text{half_average_length} \cdot 4 = N $, it follows that $ \text{index_bias} = (\text{half_average_length} \cdot 4 - 1) /2 $ . I realise this wouldn't make sense as it's not a whole number, but then there is something wrong with my interpretation of $ \frac{N-1}{2} $. What is wrong with my interpretation of this formula?


texts:
I want a highly damped highpass filter (damping of at least $2$), with a cutoff somewhere around $1\textrm{ Hz}$ ($51.2\textrm{ Hz}\quad f_s$)

How do i go about designing the filter with adequate control over the damping?

My best guess was to use a standard second order response and frequency transform it, giving

$$
\frac{s^2 \omega_0}{\omega _0^2 s^2 + 2\zeta\omega_0 ^2 + 1}
$$
This gives me the required damping (after using MATLAB to transform it to the discrete domain using Tustin's method with tustin), but the frequency response is terrible (the $\rm dB$ point is much, much higher than it should be).

Is there any way to solve this?
Can I design a FIR filter with control over the damping?
Is there a higher order $s$-domain function that still gives me control over the damping?

Notes: filter order and phase delay are not overly important to my design.


texts: Autocorrelation of a telegraph process/constant signal

texts:
I am trying to calculate the autocorrelation function for the telegraph process, but I somehow don't get the right results.

I am calculating the correlation function of a signal X which can take values 0 and 1 by calculating the spectral density $S(\omega)=|FFT(X)|^2$. The autocorrelation function is then given by the inverse fft of the spectral density $iFFT[S(\omega)](\tau)$.

If I do this the autocorrelation functions never drops down to zero for $\tau \rightarrow\infty$.

Indeed if I take a constant signal $X = 1$, the autocorrelation should be $0$. But if I use the steps mentioned above the autocorrelation function will give $1$ instead.

I really don't understand this, could you explain why I get these results?


texts: Suppose the input signal is represented as $$y=x+w$$ where $x$ is the signal to be estimated and $w$ is AWGN. I design a high pass filter $H$. What does $$H^{T}Hw$$ represent? Will it be just AWGN but attenuated or will it AWGN completely distorted. We can also assume that $\sigma_w = 0.2$


texts:
This is my first question on dsp. I am looking for something very simple here and I am sure there are libraries that will do pitch recognition and FFT but my requirement is very simple.

I have a PC and I am using the audio input for analog signal feedback and output for controlling the device. I am using windows as the OS.

I want to ensure that I have full control on the audio output and no other program or the OS will stream sounds into the same channel.
I want to output precisely at a known time interval
I want to check if the feedback matches with what is expected of the device.

Any snippet of code or an article to do the same will be useful. I can provide more details if there is something missing in my question.


texts:
I am using impulse responses from http://sound.media.mit.edu/resources/KEMAR.html

These are available as either raw or diffuse-field equalized HRTFs.

My question is: Which do I want? Raw, or diffuse-field equalized?

I am using these IRs to project single notes (generated synthetically using a ring resonator) to 12 locations around the listeners head.

π


texts: Using clustering to detect anomalies in vibration data

texts:
I'm looking into using clustering-based method to detect anomalies in vibration signals. The idea is to extract features for every single sliding time-window of a time series of normal vibration data, and then find clusters of the feature vectors. Given a new period of vibration data, the model will classify the period as anomaly if its corresponding feature vector does not belong to any of the normal clusters.

Is there any literature on under what condition such methods will fail? My intuition is that clustering does not capture time-dependence of data. For instance, if I use Gaussian Mixture Model for clustering and find two clusters $N_1$ and $N_2$, the sliding windows of my normal vibration data may well form an alternating pattern, such as $N_1 N_2 N_1 N_2...$, while an anomalous vibrational time series has pattern $N_1 N_1 N_1...$, even though the data points in the anomalous vibrational data will belong to the normal clusters.

I'm not sure if this intuition is valid, though, for natural vibration. Say, if I plan to use the clustering-based approach to detect anomalies in machine vibration, is it even physically possible to have different vibrations that have same clusters as described above?


texts: FFT question as it relates to Neural Networks / Supervised Learning Models

texts:
So I've been researching wavelet transforms and FFT.  I want to feed wavelet transforms of a 1D signal in time into a neural net and train against a target variable at each time step.  The idea being the WT/FFT can help the net pick up on data structures.

Only issue being I have no idea how to feed the coefficients of the FFT into the model because my approximation signal is length 87 while the original series had 44,000 samples.

I looked into zero-padding to restore the number of samples but that is apparently a technique used on the time-series itself before decomposing with FFT.

Any help appreciated, hard to find material on exactly this subject.

paper im recreating: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0180944