texts: Difference Between Convolution and Multiplication of Frequency Response and Frequency Spectrum

texts:
Suppose I have  impulse response like [1/25,2/25,3/25,4/25,5/25,4/25,3/25,2/25,1/25].
I did convolution with 600 samples of test signal (it seems that I did some filtering).
Then I calculated the same thing in frequency domain:

Found frequency spectrum of test signal and frequency response.
Multiplied them and via IDFT convert back to time domain.

Then I plotted difference between these two signals and I got that difference between
these to signals is about |0.4|.

Is it ok? It seems a little bit huge for me -- the difference is about one-tenth.

function generateConvolution2(x, impResp)
{
    var xVal = x();
    var irVal = impResp();

    var convolution = [];
    var data = [];

    for (var i = 0; i < xVal.length + irVal.length - 1; i++)
    {
        convolution[i] = 0;
        for (var j = 0; j < irVal.length; j++)
        {
            if (i - j < 0) convolution[i] += 0;
            else convolution[i] += irVal[j] * xVal[i - j];
        }

        data.push(convolution[i]);
    }
    return data;
}

//calculate dft magnitude and spectrum of a function
function magnitudeAndPhase(length, dftLength, func)
{
    var data = [];
    for (var j = 0; j <= dftLength; j++)
    {
        var f = 0;
        var re = 0;
        var im = 0;

        for (var i = 0; i < length; i++)
        {
            f = func(i);
            re += f * Math.cos(2 * Math.PI * j * i / length);
            im -= f * Math.sin(2 * Math.PI * j * i / length);
        }

        //calculate magnitude
        var mag = Math.sqrt(re * re + im * im);

        //calculate phase
        var phase;
        if (re === 0)
        {
            re = 1e-20;
        }
        phase = Math.atan(im / re);
        if (re < 0 && im < 0) phase -= Math.PI;
        if (re < 0 && im >= 0) phase += Math.PI;
        data.push({mag: mag, phase: phase});
    }
    return data;
}

function multDFt(length, dftLength, func1, func2)
{
    var data1 = magnitudeAndPhase(length, dftLength, func1);
    var data2 = magnitudeAndPhase(length, dftLength, func2);
    var result = [];
    //1)multiply magnitudes
    //2)add phase
    for (var i = 0; i <= dftLength; i++)
    {
        result.push({mag: data1[i].mag * data2[i].mag, phase: data1[i].phase + data2[i].phase });
    }

    return result;
}

//using magnitude and phase calculate IDFT
function IDFT2(mainLength, length, dftLength, testSignal, imprResp)
{
    var rea = [];
    var ima = [];

    var dftMagResult = multDFt(length, dftLength, testSignal, imprResp);
    //create re and im part of coefficients
    for (var i = 0; i <= dftLength; i++)
    {
        rea.push(dftMagResult[i].mag * Math.cos(dftMagResult[i].phase));
        ima.push(dftMagResult[i].mag * Math.sin(dftMagResult[i].phase));
    }

    //prepare re and im coefficients
    rea[0] /= length;
    rea[dftLength] /= length;
    for (var j = 1; j < dftLength; j++)
    {
        rea[j] /= (length / 2  );
        ima[j] /= -1 * (length / 2  );
    }
    var vala = [];
    for (var i = 0; i < mainLength; i++)
    {
        vala.push(0);
    }
    //calculate idft
    for (var j = 0; j <= dftLength; j++)
    {
        for (var i = 0; i < length; i++)
        {
            vala[i] += rea[j] * Math.cos(2 * Math.PI * j * i / length);
            vala[i] += ima[j] * Math.sin(2 * Math.PI * j * i / length);
        }
    }

    return vala;
}

UPDATE: Added my js code. Don't bother with my js style (I'm beginner) and some inconsistent naming. Every func. return array because it is needed for plotting libs.

Thanks in advance.


texts:
I've read the thread: Cut off frequency of moving average

And I use the second answer in my algorithm to calculate the 3dB cut off frequency of my filter, which works great, as my filter length is usually above 300. I verified it with the step response.

But I would like to have a source or derivation for this formula.

I tried by hand with taylor series stopping after the second and third term. I come close but not exactly to the formula and mapple gives me a valid but extremly complicated result.

Hope you guys can help.

kind regards

Slev1n


texts: Understanding the constraint on distance between microphones in microphone array for Direction of Arrival estimation

texts:
The distance between microphones in a two-microphone array for direction of arrival estimation, should be at least less than half of the wavelength of the signal that is being detected by the microphone, when the signal is periodic.

I am having trouble understanding this constraint on the distance between microphones. As far as I remember, the path difference between the sound arriving at one microphone and the other shouldn't be a multiple of half the wavelength, because if $x = \frac{\lambda}{2}$, the phase difference between the two signals for a continuous periodic signal would be zero, and thus we wouldn't be able to measure the time difference between the sound arriving at the two micropones.

However, I am unable to understand how a less than constraint appeared, and on the distance between the microphones and not on the path difference of the sound arriving at the two microphones.


texts: Avoiding transients during filtering: Padding or Gustafsson's method?

texts: I'm pretty uneducated when it comes to filtering, more of a stats guy, but I have to bandpass some signals and after a bunch of reading today I understand that the scipy.signal.sosfiltfilt function seems to be a gold standard. One thing I notice is that it doesn't provide the option to use the Gustafsson method for reducing transients as does regular ol' scipy.signal.filtfilt does. Which led me to wonder, is there much to argue for using the Gustafsson method over simple padding?


texts: so I've been attempting to do a FFT on some data and I'm seeing a peak at 0Hz, which I can't really comprehend why since I'm quite new to signal processing. I read in my file using pandas, saved them to arrays and used scipy's FFT function on the y data. But when plotting it against a frequency linspace the results seem odd to me. I'm not quite sure how to get the time domain into the frequency domain. Any help would be much appreciated, thank you. (Note: The data isn't perfectly horizontal).


texts:
I have a question regarding the scrambling code, What's the difference between spreading and scrambling? As I know, following the below figure, the scrambling doesn't affect the signal bandwidth:

As shown above, the chip rate after the scrambling code is still the same,I mean it's not spread. My question, does scrambling is performed only by multiplying the sequence resulted from spreading, element by element, with scrambling code? or how is it processed ?

Thank you


texts: Possible ways to determine the direction of moving object?

texts:
I'm trying to determine the moving direction of an object in x-axis only. Figure below shows part of my effort.

However, I got stuck. As you can see from the middle image, the green represent the dot moved to right (positive x), the black represent overlapped area and the purple represent the reference image relative to the moving image.
In bottom image, I got two directions for that one object. Left one is going through positive transition, and right negative transition. I mainly use the intensity differences to determine the direction., i.e. Moving image intensity - reference image intensity.
Is there other possible ways to determine the direction of a moving object (to and fro in x-direction only)? Or is there some ways of correction I can apply to?
Thanks in advance guys!


texts: Understanding Digital Audio - "From Audio Device" in Matlab

texts:
I understand that digital audio has both a sample rate (time "discretization") and a bit depth ("amplitude" discretization). Also, it can have multiple channels (stereo would be 2 channels).

But things that I can't find readily in the internet:

What is the range that is discretized in the bit depth? My experience is in motions, so I know if I have an accelerometer with a range of +-10 g's, and I discretize it in 8bit, well I would have a resolution of 20g/256. What would be the equivalence in audio?
Matlab's Simulink has a block for reading audio. It's output has the dimensions of MxN(xTime), where M is the number of "consecutive samples" and N the number of channels. What is "consecutive samples" (also called in the same help file as "Frame Size (samples)"). If I set this value to one, I would get what I would expect from a Mono sound, just a time history of amplitudes. But if I set it to two, it's almost as if they are both the same with minor differences in fast transitions. Is this meant for to get a average of the signals to be somewhat more "true to sound"?

Here is the link for Simulink's block: http://www.mathworks.com/help/dsp/ref/fromaudiodevice.html


texts: What is the meaning of imaginary values in the time domain?

texts:
Say I take a real valued function in the time domain, take the Fourier transform to get its frequency domain representation, do some work on it in the frequency domain, which leads to imaginary values, and then take the inverse Fourier transform to get a new time domain signal.  What does it mean if the values of the time domain signal are imaginary?  I understand that imaginary values in the frequency domain allow for the representation of phase of each frequency (https://stackoverflow.com/questions/10304532/why-does-fft-produce-complex-numbers-instead-of-real-numbers).  But what do imaginary numbers mean in the time domain?

I can tell you what my particular application is, but I believe that this is a general question that others might come across in the future.


texts:
in a mock exam of dsp was asked this question, and even in knowing how to compute DFT's of the plot, I wondered if there was a fast way of pairing these plots to their respective impulse response?

Thanks for reading!


texts: PSD Periodogram smaller segments compared to larger segment

texts:
I am calculating the power spectrum of a heart rate variability via FFT and a Welch window.

I have a sample of 500 points as the input signal.
The signal was broken up into two 128 point segments. For each segment the multiplication with a window, FFT and scaling was performed.
Results from segment 1 and segment 2 were analyzed.

These data were compared to the same analysis with one larger segment of 256 point.

As a result the input data from values 1-128 and 128-256 are compared the values from 1-256.

I would have expected that the averaged result for two samples 1-128 points and 128-256 points would be the same than the result for one sample with 1-256 points.  However, this is not the case.

The value of the spectral frequency bands with high frequencies is much larger for each of the two samples than for the one sample (factor 1.7 larger). Lower frequencies are better fitting, but not exactly.
(Smaller Segment 1: power 3500 ms^2, Smaller segment 2: 2700 ms^2, Larger Segment: power 1700 ms^2).

I compared my self-made calculation with a commercial calculation and the one sample variant produces correct values.

Does anyone has an idea what could be the reason for this phenomenon?
What could be the error? How should it normally look like?

It is also not so clear how the add or average the two smaller spectra.

The scaling and binning seems to be correct in each case.


texts: Math.Net and alglib returning different FFT outputs by default

texts:
I am developing an application in C# with spectrogram drawing functionality.

For my fist try, I used MathNet.Numerics, and now I am continuing to develop with alglib. When I changed from one to the other, I noticed that the output differs for them. Mathnet uses some kind of correction by default, which alglib seems to omit. I am not really into signal processing, also a newbie to programming, and I could not figure out what the difference exactly comes from.

MathNet default output (raw magnitude) values are ranging from ~0.1 to ~274 in my case. And with alglib I get values ranging from ~0.2 to ~6220.

I found that MathNet Fourier.Forward uses a default scaling option. Here is says, the FourierOptions.Default is "Universal; Symmetric scaling and common exponent (used in Maple)." https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/FourierOptions.htm If I use FourierOptions.NoScaling, the output is the same as what alglib produces.

In MathNet, I used Fourier.Forward function: https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/Fourier.htm#Forward In case of alglib, I used fftr1d function: https://www.alglib.net/translator/man/manual.csharp.html#sub_fftr1d

What is the function that I could maybe use to convert alglib output magnitude to that of MathNet, or vice versa?
In what cases should I use these different "scalings"? What are they for exactly?

Please share your knowledge. Thanks in advance!


texts: Bounds of higher order derivative of bandlimited signals

texts:
Here is a bandlimited function f(t) with bandwidth Ω:

The function f(t) is bounded in [-A,A].
Then the bound of the derivative of f(t) is bounded as:
|f'(t)|≤2πΩA.
So, what is the bound of its n-th order derivatives?
Is it as follows?

where n is the order of derivative.
If yes, then are the bounds of higher order derivatives much larger than the bound of f(t)?
I tried MATLAB simulations but the results were contrary.


texts:
I have two wave files each containing a sample of a note. I wish to compare, using a program, each sample and present how loud the first sample is with respect to the second if each sample is played back to the listener with the same physical conditions (volume, distance from speaker etc.). For example, if the second sample is twice as loud as the first, I wish to present "2" to the user. I'm not interested in being very accurate. I simply want to give a reasonable estimation.

Lets say the first file contains only a 1kHz tone with an amplitude of 0.1 (0.1*sin(2*pi*1000*t)), while the second contains a 1kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*1000*t)). So the difference is 10db correct? Does this mean that the second sample will be, roughly speaking, perceived as being twice as loud, with respect to the first?

I understand that phons are with respect to spl and that therefore if I have only the two wave files to work with, I can't talk in such terms since I've no idea what spl existed at the microphone when the recordings were made - or more importantly, the spl that exists at the listeners ears when they hear the recording. However, lets say I have a 3kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*3000*t)) in my first wave file, while the second file contains a 10kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*10000*t)). Isn't it fair to say that the first file will be perceived as being louder? Should it be fair, based on the graph to say that the 3kHz tone will be perceived as being roughly 4 times louder? Since in the graph you have to, roughly speaking, increase the 10kHz signal by 20db to get the same loudness as that of a 3kHz signal.

Again lets say that the first file contains a 3kHz tone with an amplitude of 0.1 (0.1*sin(2*pi*3000*t)), while the second contains a 10kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*1000*t)). Is it fair to estimate that the difference in loudness is a factor of 0.5? That is, the first file will be perceived as being twice as loud as the second.


texts: Fourier Transform of Kernel Density Estimation - Convolution Theorem?

texts:
I am reading this paper about density estimation (Appendix A), where the authors apply a Fourier transform to the estimated probability density (the $X_j$ are a sample of $N$ data points drawn independently from an unknown probability density function $f(x)$):

$$
\hat{f}(x) = \frac1N \sum^N_{j=1} K(x-X_j)
$$

The authors now note that they apply the convolution theorem and the Fourier transform (where the Fourier transform of the function $g(x)$ is defined by $\Phi_g(t) = \int \exp\lbrace\mathrm{i}xt\rbrace g(x) \mathrm{d}x$), and obtain:

$$
\Phi_{\hat{f}}(t) = \kappa(t) \frac1N \sum^N_{t=1} \exp\lbrace\mathrm{i}tX_j\rbrace
$$

where $\kappa(t)$ is the Fourier transform of the Kernel $K$.

Now, I am a bit stuck on how they obtained this form and how exactly they used the convolution theorem here. In particular, what happened to the argument $X_j$ in the kernel?

EDIT: Ok, I think I got it. Is this the right approach?

We apply the Fourier transform and obtain:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \int K(x-X_j) \exp\lbrace\mathrm{i}xt\rbrace\mathrm{d}x
$$

We can easily rewrite this equation as an integral over a Dirac delta distribution:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \iint K(y) \delta(y-(x-X_j)) \exp\lbrace\mathrm{i}xt\rbrace\mathrm{d}x\mathrm{d}y
$$

Integrating over $x$ and using the symmetry of $\delta$ we obtain:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \int K(y) \exp\lbrace\mathrm{i}(y+X_j)t\rbrace\mathrm{d}y
$$

We see that we can move factor involving $X_j$ outside the integral, and are left with the result:

$$
\Phi_{\hat{f}}(t) = \frac1N \sum^N_{j=1} \exp\lbrace\mathrm{i}X_jt\rbrace \int K(y) \exp\lbrace\mathrm{i}yt\rbrace\mathrm{d}y
$$

where the first factor is $\Delta(t) = \frac1N \sum^N_{j=1} \exp\lbrace\mathrm{i}X_jt\rbrace$, and the second factor is $\kappa(t) = \int K(y) \exp\lbrace\mathrm{i}yt\rbrace\mathrm{d}y$, exactly as in the paper.

I am still not sure in how far they applied the convolution theorem here?


texts: Scale FFT frequency range for a bars graph

texts:
I'm trying to show frequency spectrum for audio data as a bars graph using FFT.

If I have data for a sine wave at 440 Hz and with sampling rate 44100 Hz and analyze its first 1024 samples I get an array with 1024 FFT values.
I know that the first half of those values (512) represents the frequency range equal to half the sampling rate (22050 Hz).

So to draw a bars graph I just take the first 15-20 values from the FFT array, calculate the frequency range for each element and plot it.

In the case of FFT size 1024 I frequencies each ~43 Hz wide.
But if I change FFT size to 2048 each frequency range is ~21.5 Hz. So the peak of the spectrum of the sine wave will be outside my 15-20 bars.

So I wanted to ask if it's possible to scale/normalize the X axis for a frequency range and a constant graph bars count?

The example limiting data would be:

frequency range 0-500 Hz
16 bars on the graph

I guess it's some kind of bins/magnitudes grouping. But should I sum them or get max? And what could be done with magnitudes when grouped bin number isn't integer? (I hope my explanation is not too confusing, it's hard to choose the correct terminology)

Thanks.

Here is my test code in Python:
(I used a static int array imported from another file that represents a sine wave)

import scipy.io.wavfile as wavfile
import scipy
import scipy.fftpack
import numpy as np
from scipy.signal.windows import *
from matplotlib import pyplot as plt

from imp.sine_440_hex import *

def run():
  # fp = 'e:/sine_440_hz.wav'
  # FS, data = wavfile.read(fp)

  # data - array of int for a sine wave at 440Hz and amplitude 0.5
  DS = len(data)
  FS = 44100

  FN = 1024
  # FN = 2048
  MAXS = 32768

  MN = 15   # number of first frequencies in the FFT result

  ndata = np.array(data)
  ndata = ndata / MAXS

  w = np.hamming(FN)
  y = ndata[:FN] * w

  freqs = scipy.fftpack.fftfreq(FN, 1/FS)
  mags = abs(scipy.fft(y))

  dbfs = 20 * np.log10(mags * 2 / np.sum(w))
  dbs = dbfs + 120

  # -------- plot --------
  fr_res = freqs[:MN]
  db_res = dbs[:MN]

  ax1 = plt.subplot(211)
  plt.plot(fr_res, db_res)
  plt.grid(True)
  plt.xlabel('Frequency [Hz]')
  plt.ylabel('Amplitude [dB]')

  ax2 = plt.subplot(212)
  plt.bar(fr_res, db_res, fr_res[1]-10)
  plt.grid(True)
  ax2.set_ylim([80, 120])
  ax2.set_xticks(fr_res)
  ax2.set_xticklabels(fr_res.astype(int))

  plt.show()

# -------
run()

And example graphs:

For FN = 1024

For FN = 2048


texts: How to remove historical values from exponential filter

texts:
Suppose I have two time series $A:=[0,0,0,4,5,6]$ and $B:=[1,2,3,4,5,6]$. I implement the following filter (with initial value at the first element):

$$ f(t) = f(t-1) + \alpha*(f(t)-f(t-1))$$

with $\alpha \in [0,1]$.

Is there any way I can use filtered values of $B$ to recover the filtered values of $A$? For example suppose I can store all historical filtered values of $B$ for any choice of $\alpha$ in memory, is it possible to reconstruct what the filtered values for $A$ would have been for some choice of $\alpha$?

Intuitively this would involve subtracting the some filtered value(s) of $B$ at time step $3$ from the filtered value(s) of $B$ at time step $6$, to recover an approximation of $A$.


texts: I have run an FFT now want to get the mean value in the 8hz - 13hz range?

texts:
I have a  matrix [66 x 100,000] (it's some EEG data, across 66 electrodes) and want to find the mean power in the 8-13hz range. How can I do this?

I have so far: x = fft(data) What do I have to do next?

Furthermore, when I plot this I get the following::

I was hoping for something more like this:


texts: Multirate filtering attenuates frequencies at overlap

texts:
I am trying to build a filter with different frequency bands in Matlab. I would like the filter to be as invisible as possible, but after the filtering, a great amount of frequencies are being attenuated and I believe it is where the filters overlap, but also one can see that the low frequencies are also attenuated though there is not any overlap there.

The following images show the results:

Original signal

Filtered signal

Could you help me understand what I am doing wrong?

This is the code I am using:

[x,Fs] = audioread('drum-loop.wav');

% define constants  - - - - - - - - - - - - - - - - - - - - - - - - - - - -

ntaps = 255;                                                               % number of filter-taps

% sampling frequencies for each branch
fs_A = 48000;                   %
fs_B = 24000;                   %
fs_C = 6000;                    %
fs_D = 1500;                    %

% decimation for each branch
dec_A = 1;                      %
dec_B = 2;                      %
dec_C = 4;                      %
dec_D = 4;                      %

% decimation filters  - - - - - - - - - - - - - - - - - - - - - - - - - - -
x_B = decimate(x, dec_B);
x_C = decimate(x_B, dec_C);
x_D = decimate(x_C, dec_D);

% interpolation filters   - - - - - - - - - - - - - - - - - - - - - - - - -
% ref: https://es.mathworks.com/help/dsp/ref/mfilt.html
intfilter_B = mfilt.firinterp(dec_B);                                      % decimate signal for branch B
intfilter_C = mfilt.firinterp(dec_C);                                      % decimate signal for branch C
intfilter_D = mfilt.firinterp(dec_D);                                      % decimate signal for branch D

% define filters  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ovlp1 = 2000;
ovlp2 = 500;
ovlp3 = 250;

% branch A: band pass filter.
bp_A = designfilt('bandpassfir','FilterOrder',ntaps-1, ...
         'CutoffFrequency1',fs_B/2-ovlp1,'CutoffFrequency2',(fs_A/2-1), ...
         'SampleRate',fs_A);
% branch B: band pass filter.
bp_B = designfilt('bandpassfir','FilterOrder',ntaps-1, ...
         'CutoffFrequency1',fs_C/2-ovlp2,'CutoffFrequency2',(fs_B/2-1), ...
         'SampleRate',fs_B);
% branch C: band pass filter.
bp_C = designfilt('bandpassfir','FilterOrder',ntaps-1, ...
         'CutoffFrequency1',fs_D/2-ovlp3,'CutoffFrequency2',(fs_C/2-1), ...
         'SampleRate',fs_C);
% branch D: band pass filter.
bp_D = designfilt('lowpassfir','FilterOrder',ntaps-1, ...
         'CutoffFrequency',(fs_D/2-1), 'SampleRate',fs_D);

% apply filters
x_A_f = filtfilt(bp_A,x);
x_B_f = filtfilt(bp_B,x_B);
x_C_f = filtfilt(bp_C,x_C);
x_D_f = filtfilt(bp_D,x_D);

% summation filter
y = x_A_f + interp(x_B_f,dec_B) + interp(interp(x_C_f,dec_C),dec_B) + ...
    interp(interp(interp(x_D_f,dec_D),dec_C),dec_B);
y = y/max(y);                                                              % normalize

% write audio
audiowrite('output.wav',y,48000);

% plot figures - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% compute fft
figure;
subplot(3,1,1);plot(x);title('original signal')
NFFT = 2^nextpow2(length(x));                                              % Next power of 2 from length of y (in samples)
Y = fft(x,NFFT)/Fs;
f = Fs/2*linspace(0,1,NFFT/2+1);
subplot(3,1,2);spectrogram(x,blackman(128),60,128,1e3,'yaxis');title('spectrogram')
% subplot(3,1,3);plot(f,2*angle(Y(1:NFFT/2+1)));title('phase')
subplot(3,1,3);plot(f,abs(Y(1:NFFT/2+1)));title('fft')

% compute fft
figure;
subplot(3,1,1);plot(y);title('filtered signal')
NFFT = 2^nextpow2(length(y));                                           % Next power of 2 from length of y (in samples)
Y = fft(y,NFFT)/Fs;
f = Fs/2*linspace(0,1,NFFT/2+1);
subplot(3,1,2);spectrogram(y,blackman(128),60,128,1e3,'yaxis');title('spectrogram')
% subplot(3,1,3);plot(f,2*angle(Y(1:NFFT/2+1)));title('phase')
subplot(3,1,3);plot(f,abs(Y(1:NFFT/2+1)));title('fft')


texts: How to prove that Dilution of a time domain signal results in duplication in Frequncy domain?

texts:
While reading a description on FFT from this book, I found the following statement:

  ...diluting the time domain with zeros corresponds to a duplication of
  the frequency spectrum.

How can I prove this intuitively or otherwise? I have read some properties of Fourier Transform, such as a  circular shift in time domain and frequency domain etc, but I am unable to deduce how duplication in frequency domain will happen. Following illustrate this:

Time domain samples

x(0)=a
x(1)=b
x(2)=c
x(3)=d

Frequency domain samples

X(0)=A
X(1)=B
X(2)=C
X(3)=D

Dilution with 0 in time domain

x(0)=a
x(1)=0
x(2)=b
x(3)=0
x(4)=c
x(5)=0
x(6)=d
x(7)=0

Frequency domain samples having duplicate values

X(0)=A
X(1)=B
X(2)=C
X(3)=D
X(4)=A
X(5)=B
X(6)=C
X(7)=D

Frequency domain samples have duplicate values. Why?


texts:
A radio signal recording of a wireless communication system (e.g: Wi-Fi traffic) is beaconized, channelized and subject to noise.

When working with such an RF signal, numerically transformed to a stochastic timeseries, to which extent can I consider the signal to be stationary?
What properties are relevant from a communications point of view?


texts: Recognizing object of interest in object tracking using Lucas-Kanade method

texts: I am able to find out the velocities of pixels using the Lucas-Kanade method. Now I want to use it for tracking. I am using a Kalman filter taking $x$, $y$, $v_x$ and $v_y$ ($(x,y)$ is position of a particular object pixel and $(v_x,v_y)$ is velocity of a pixel) as state variables of the Kalman filter. Now the problem I am facing is how will I get to select the pixels that are representing the object to be tracked?


texts: Relationship between sampling and number of pixels?

texts: In case of digital images, what is the relationship between sampling and number of pixels? For upsampling, the number of pixels increass, and for down sampling, the number of pixels decrease??


texts:
my book about RADAR says that:

If a Continuous wave radar sends a sine wave at frequency $f_T$ to a moving object (at speed V), a frequency $f_R$ is received. Their difference is called Doppler frequency and it is equal to $f_D = f_R - f_T = - 2 \frac{V}{\lambda}$.

It is quite clear. Then, it says that the situation is slightly different if a radar which sends rectangular pulses modulated with a sine wave at frequency $f_T$ is used. In such a case, the continuous wave is sampled by a train of strict rectangular pulses spaced by period PRT (Pulse Repetition Time):

Then, it says that:

If $PRF = \frac{1}{PRT}$ (Pulse Repetition Frequency) is less than $2f_D$, spectral aliasing will occur (Nyquist Theorem).

First Question: why is the Nyquist theorem applied to the doppler frequency $f_D$? It should be applied to the sampled sine wave, which is at frequency $f_T$.
Then, it shows the aliased spectrum:

and says that:

The spectral aliasing causes spectrum folding at $f_D > \frac{PRF}{2}$ in the main interval ranging from $-\frac{PRF}{2}$ to $+\frac{PRF}{2}$

Second Question: why does aliasing cause folding? As far as I'm concerned, spectral aliasing means that the original (correct) spectrum is replicated and those spectra overlap (and hence it is not possible to separate them correctly).

I don't understand the link between my last two pictures.


texts: Signals and systems book with linear algebra approach

texts:
I am currently taking a course on discrete-time signals and systems. I am using B. P. Lathi's Signal Processing and Linear Systems, which I don't like at all, as it doesn't draw any paralells to linear algebra. The concept of eigenvalues and eigenvectors in linear transformation is never discussed, for example.

Can you please recommend a book on signals and systems with a linear algebra approach?


texts: Get the discrete-time poles and zeros from continuous-time poles and zeros

texts:
How do you implement the following function:

$$[Z^d, P^d, K^d] = \text{fcn} \,(Z^c, P^c, K^c),$$

where $Z^c = [z^c_m, ..., z^c_1]$, $P^c = [p^c_m, ..., p^c_1]$, and $K^c$ are zeros, poles, and gain of the continuous-time TF, and similarly $Z^d$, $P^d$, and $K^d$ are zeros, poles, and gain of the equivalent discrete-time TF with sampling time of $T_s$.

I have tried the following but it did not give me the correct answer:

$Z^d = e^{Z^c \times T_s}$

$P^d = e^{P^c \times T_s}$

$K^d = 1-e^{-K^c \times T_s}$    ---> (is the minus sign right?)

I am looking for some conversion formula that can handle all cases (including single/multiple poles or zeros).

Note: I do not want to use any Matlab built-in functions like c2d(), zp2tf(), or tf().


texts: Approximation using a Fourier transform with low pass filter

texts:
I need to approximate a function f, but I cannot do so with frequencies that exceed 1kHz

What is the best approximation I can get? Is taking the Fourier transform then zeroing any term above 1kHz the best approximation? Or can I fiddle with the lower order terms and get a better fit?

This is a made up scenario, but I have to prove the same concept with Walsh transforms. I am fairly certain that the lower order terms form the best approximation from random twiddling and hill climbing searches, but I need proof.

I believe the proof is something very similar to a least squares regression proof, but I can't get it. Has this problem been solved before? At least in the Fourier domain?


texts: I have a recording of a single note played several times on a violin. I'm working with Matlab to estimate the amplitude of the first N harmonics. I already know which note is played but there can be fluctuations in frequency and noise (very little, actually). Can you suggest me the more robust method in this kind of situation? At the moment I am windowing the signal and using a peak detection algorithm on the FFT of each frame in order to find the amplitude of the harmonics. Is there anything better?


texts: What's the meaning of the continuity in spectrum analysis?

texts:
We all know that for any, suitable, kind of signal $f(t)$,  there corresponds a Fourier transform function $F(j\omega)$ such that
$$
f(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} F(j\omega) e^{i\omega t} d\omega
$$

So, $F(j\omega)$ can be recognized as the component of $f(t)$ along $e^{i \omega t}$. This is similar to orthogonal analysis. Like a cartesian vector $\vec{a} = (1,2,3)$ can be represented as
$$
\vec{a}=  1\hat{x}+2\hat{y}+3\hat{z}
$$
where  $\hat{x}$, $\hat{y}$, and $\hat{z}$ are the orthonormal basis vectors (unit vectors along x-y-z directions).

However, these two things are different in one perspective.Among $x$,$y$,$z$, there is no concept as 'distance'. For example, in Fourier transform. we can say that the distance between $\omega_1 = 1$ and $\omega_2=2$ is smaller than the distance between $\omega_2$ and $\omega_3 =4$.

This algebra structure in Fourier transform definitely has some extra meaning, but what is it? Does it mean that the Fourier transform makes "more sense" compared with other orthonormal basis analyses where there is no easily-definable distance between basis functions?


texts: How to shift a signal on the $y$-axis to match another signal in MATLAB

texts:
I have two signals recorded in the time-domain then converted to the frequency domain. These signals are $z$-axis accelerometer measurements recorded from an Arduino 101.

Now visually they are fairly similar (with the peaks and falls sort of lining up), however the difference is in their amplitude, e.g.:

I will be comparing the power of the two signals, so I need them to be in-line with each other, ie. shift either one of the signals to fall around the same area.

I'm not too sure what the terminology for this is, but any help would seriously be appreciated.


texts: What is a clear, definitive reference for learning about WOLA?

texts:
I understand that WOLA stands for Weighted OverLap Add, and that it is used as a channelizer for STFT (short time Fourier transform). I am comfortable with DFT and FFT, and the concept of performing an FFT in a window which slides along in time, in order to get the frequencies of the signal at different points in time.
However I am looking for a good, clear, definite reference for learning about WOLA in more depth, which contains the relevant equations, some diagrams/further intuition, an explanation of why WOLA is better than just taking blocks of samples, applying a window, and taking the FFT, and ideally also some (pseudo) code. Information about the inverse transform (synthesis) would also be very useful. Ultimately I want to be able to implement this for myself, while understanding how it works, and why it is a sensible thing to do.
I have tried looking for this already, but was surprised not to find what I was looking for. Maybe WOLA also goes by a different name?
Reference to a chapter of a DSP book, a website, or a paper would be great. Alternatively if someone is willing to provide all of this information directly as an answer to this SE question, that would also be great.


texts: How to Clean Black and White Images Which Are Results of Edge Detection

texts:
I want to extract text from some images with messy backgrounds.

First step was switching to gray scale, applying an edge detection algorithm, and then cleaning up a little bit.

Here's an example of the result:

I want to use the tesseract OCR to extract the text from it, but it's not good enough to find it, as the background is quite messed up still. (I tried with a cleaner version and it's doing ok).

So, could anyone recommend an algorithm, or a type of algo, that would help me clean up this image and just get the text? (it will always have this font).

One idea I had was to try to remove all the clumps that are "larger" than the font, which is pretty thin.

ps: hopefully it's the right stackexchange site.


texts: How do I implement a nonorthogonal quadratic spline wavelet into discrete wavelet transform?

texts:
I am trying to use the discrete wavelet transform for signal processing (time series data from a plant). I would like to use a mother wavelet that is not in Matlab. Mallat and Zhong 1992 described a nonorthogonal quadratic spline wavelet I would like to use. They name some coefficients and impulse responses, but I'm not sure how to use these (reference below). To detect and treat abnormal sudden changes, the wavelet must satisfy a very particular property. We call a smoothing function any function θ(x) whose integral is equal to 1 and that converges to 0 at infinity. This function must be twice differentiable and define ψ_1 the first order derivative and ψ_2 the second order derivative.

ψ_1=(dθ(x))/dx  and ψ_2=(d^2 θ(x))/(dx^2 )

The wavelet transform is computed by convolving the signal function f with a wavelet function (dilated at a scale s).

W_s^1 f(x)=fψ_(1,s) (x) and W_s^2 f(x)=fψ_(2,s) (x)

The extrema points of the first wavelet transform (WT) gives us the zero crossings points of the second WT. The extrema of W_s^1 f(x) and the zero-crossings of W_s^2 f(x) are the inflection points or the sharp variations points of f(x) smoothed by θ(x).Those edges in the signal can only be found if the wavelet is the second derivative of a smooth function. Furthermore, the first- and second-order WT are used in the detection of steady-state periods.

So I would like to code a quadratic spline wavelet function and a cubic spline scaling function like those in Mallat and Zhong (1992) article. Could anyone help me ?

Mallat, S., & Zhong, S. (1992). Characterization of Signals From Multiscale Edges. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(7), 710-732. doi:10.1109/34.142909


texts:
The Burg algorithm implemented in PRAAT for formant estimation also works via LPC. First the sound is resampled to a sampling frequency of twice the value of maximum formant. Then a pre-emphasis is applied and finally, PRAAT applies a Gaussian like window to compute linear prediction coefficients through the Burg algorithm.

I understand that pre-emphasis is helpful for formant extraction by LPC since the human voice is not spectrally flat, but am confused about why Praat is resampling.
https://github.com/praat/praat/blob/bac77ad76481c65703ade8d7ddb190e9d67048d6/fon/Sound_to_Formant.cpp#L381


texts: How to achieve a smoothed frequency phase response with no outliers

texts:
I have two sampled signals, which I then perform wavelet analysis on, by synthesising the Filterbank, perform the dft on each kernel, and the dft on the sampled signals, and multiply together.

Using this I can get successfully extract the magnitude of the signal, and either take an average or identify the time response of a frequency.

Now I would like to extract the phase.  Taking either an instance in time, or averaging the mean signal I can get the phase and plot it, using atan2(im / re) for each signal, and then subtracting one from the other. However it appears to be very noisy.  Although I can identify the phase trend, there are a lot of sporadic results that lie far out of the trend.

Researching into this a little, I have found there can be issues with rounding small floats, and it is good to set a threshold for the phase, so everything below the threshold is rounded to 0.  This returns similar results.

What I am witnessing looks like some kind of phase noise, reading this it is my understanding this occurs when the frequency is not an exact integer number of cycles for the dft size. This is where I start to get confused as I have my Filterbank frequencies and I have the dft size. How do these two relate?   Would I need to remove all results where the centre frequency of each filter does not meet an integer number of cycles for the sample length?

Am I on the right path here, or is there something else that would be causing my phase noise.  Failing this, is there any decent way to filter out outliers in a signal.

If I go about not performing Wavelet analysis, and do a straight forward dft on an incoming signal and extract the phase, I get similar results, so I think phase noise is along the right path,  please point me in the right direction!


texts: Why is the Hamming and Hanning window the most widely used for speech spectograms?

texts:
When taking the Short-Time Fourier Transform of the speech signal, windowing is applied to taper away the signal frame at both ends and mitigate introduced artifact at expense of the maximum obtainable resolution. Window design is then motivated controlling the trade-offs or spectral-time resolution, and various figures of merits, like the mainlobe width, sidelode width, sidelobe roll-off per octave is introduced in the literature.

In terms of a speech signal, usually a Hamming or Hanning window is employed, and the only justification I have found for it is that speech signals are narrow band signals. How could I demonstrate that a Hamming window is indeed better than, say, a rectangular window for a speech signal?

In a typical natural speech application, a speaker will exhibits frequencies between 0 and 200 Hz, a little higher perhaps for females. Utterances are usually a lot longer than that to produce required frequency, so we can afford to take more frames. I don't know a required resolution for speech. Thus, is it always better to use a Hamming/Hanning window on speech? Is there a clear reason I should tell someone using a rectangular window on producing spectograms, why he/she should not do it (apart from reproducibility of certain results in academic literature)?


texts: Confusion calculating length of Kaiser window

texts:
I am reading Digital Signal Processing Using Matlab by Ingle and Proakis (3rd ed.)
In Chapter 7, Section 7.3,
I am confused why he is adding 1 to main lobe width for calculating length of window.
I have attached snapshots below.
First/above snap is of matlab code of example 7.8 and 2nd snap is of code of example, both snaps contain the confusing line as highlighted.


texts:
Possible Duplicate:
  Free alternatives to SIFT/SURF for commercial use

There are same very good algorithms for interest point detection/description, like:

SIFT (Scale-Invariant Feature Transform)
SURF (Speeded-Up Robust Features)
MOPs (Multi-Scale Oriented Patches)

Regrettably, at least the first two are patented and cannot be used commercially without paying huge fees to the authors. I am currently negotiating about the third one, which is patented as well.

In a search for patent-free algorithm, I found there is a multitude of patents like this one. Huh?! Is there merely anything usable then?

Maybe at least Harris-Laplace or some basic detector can be used. But how to be sure?


texts: Spectral Leakage, Phase and other questions about FFT

texts:
I'm currently implementing a voice recognition algorithm using FFT and now some questions came up to me:

Spectral Leakage: I know what it is, why it appears and how to solve it using a Windowing Function, but is there another work around for this issue?
Phase: I know what it is, but where can I use it and how?
Getting the real amplitude (I know it has something to do with spectral leakage, but I will ask this separately): Let's say I have a frequency spectrum with a frequency resolution of 250 (sampling rate = 32000 Hz, frame size = 128) and the input I sent to FFT is a sinusoidal function whose frequency is 567.5 Hz (thus, the number of cycles are non-integer, causing the spectral leakage). My spectrum would be something like this:

So my question is: How can I retrieve the original amplitude of the frequency 567.5 Hz?


texts: Fourier transform and impulse function $\delta(\omega)$

texts: Why does impulse function $\delta(\omega)$ keep occurring in the Fourier transform expression of standard functions like $\sin(t)$, $\cos(t)$, constant function, unit step $u(t)$ etc? (can someone intuitively explain what it means)


texts: Where does the error come from in PID control?

texts:
Let's say my thermostat is 50 degrees Fahrenheit. I know that 1000 units of kinetic energy are needed to raise the temperature from 80 degrees Fahrenheit. So I calculate how much I need to increase my furnace to produce 1000 units of kinetic energy, and then the internal temperature should be 80 degrees.

What's wrong with what I just described? There's no PID control involved, no error term.


texts: How do KLT and SVD differ from each other?

texts: What are the differences between KLT and SVD? It looks more or less similar to me. Did I misunderstand the concepts? Explain me more about SVD and if possible provide me the link for SVD (Image Processing). Thank you.


texts: Why DFT is used for approximating CTFT when you can approximate CTFT-integral itself?

texts: I was using MATLAB for approximating FTs. Why DFT is used if we can approximate the transform-integration using summation.


texts: How to get around the circular shift property of Discrete Fourier Transform?

texts:
I understand that when we introduce a linear time shift using DFT on a finite sequence, the algorithm assumes that the signal repeats itself outside of the given range. Here is an example explaining this circular shift property of DFT (Oppenheim, 1998):

Now say if I sampled a finite length (X samples) of a pure sinusoid, but the sinusoid itself is not perfectly repeating when joining the beginning and end together. Hence, when I want to introduce a time shift in the frequency domain, the  DFT algorithm will assume that my signal repeats itself every X number of samples. So I get this:

The output signal (red) is not a smooth sinusoid because my original signal does not repeat itself perfectly.

My question is, is there a way for me to apply the time shift to the pure sinusoid without having this discontinuity in the output? That is, can I have another algorithm that continues my signal smoothly for me after the time shift has been applied?

Thank you all for your time!


texts: How do I implement IFFT to filter low frequencies from a group,after Freq detection by FFT?

texts:
I use a combination of Visual Basic Programming interfaced with FAMOS which is a DSP software like MATLAB but much easier. I am trying to separate signals from a group, which requires filtering each signal. I used FFT and then IFFT on a 600s signal but the IFFT was longer than the original.
I don't understand if this changes the originality of the signal in terms of shape,which might affect the further processing of the signal.
After re-sampling the original 10min signal according to Shannon-Niquist theorem and using rectangular windowing, I detect the frequency using FFT.

After I apply IFFT to the signal, does this change the phase or amplitude even slightly that might vary from the original signal?
How do I filter one particular freq before IFFT after detecting it from FFT (Maybe I use peak detection)? Is there a general method of doing it?
DSP is totally new to me,and I am still attempting to learn.


texts: We are  usually asked about the full system bit rate. However, we use channel encoding in practical systems, for example conventional coding with rate $(1/3)$. My question is about the full system bit rate in such system which employs channel coding. Should it be calculated including the coding bits or after extracting the code bits? I mean if we have system whose bit rate of $12bps$ and employs channel encoding of rate $(1/3)$. Should the bit rate be $12bps$ or $4bps$ since the other $8$ bits are used for coding?


texts: Time-reversal and delay commutativity and resulting signals

texts:
I'm reading Schaum's DSP book, and in Fig 1-3 they demonstrate why shifting and reversal are order-dependent, showing a couple of simple systems, 1) delay followed by reversal, and 2) reversal followed by delay.

The equations they present as follows:

1) (Time delay -> Reversal)

$x(n) \rightarrow x(n-n_0) \rightarrow x(-n-n_0)$

2) (Reversal -> Time delay)

$x(n) \rightarrow x(-n) \rightarrow x(-n+n_0)$

I'm fine with shifting and reversal not being commutative due to the different results above, but shouldn't the results be opposite? I.e.:

1) (Time delay -> Reversal)

$x(n) \rightarrow x(n-n_0) \rightarrow x(-(n-n_0)) = x(-n + n_0)$

2) (Reversal -> Time delay)

$x(n) \rightarrow x(-n) \rightarrow x(-n - n_0)$

I didn't find this in the errata for the book which is why I'm thoroughly confused. Which is correct?


texts: 'xcorr' function in matlab not working for zero signal

texts:
I came across an interesting problem when using 'xcorr' in MATLAB. The expected signal and generated signal both are at zero DC. So, I expect correlation coefficient to be 1 or perfectly correlated. To my surprise, It is showing the highest error and time lag !!! Pls help !

Thanks,


texts: Finding the $\mathcal Z$-transform of $((-\frac{1}{3})^n + \frac{1}{2})^n \mu[n-2]$?

texts:
What is the $\mathcal Z$-transform of $\left(\left(-\frac{1}{3}\right)^n + \frac{1}{2}\right)^n \mu[n-2]$?

Doing raw computations with large $n$ gives a sum of $0.6030$ which doesn't seem right, and I'm not sure how to manipulate the expression to get some clean result.


texts: Symbol timing synchronization using a high sampling rate

texts:
Symbol timing synchronization usually relies on a digital PLL procedure, an iterative process that relies on detecting a timing error signal that is used to control an interpolation filter that in its turn tracks this error until lock is achieved.
 Suppose this feedback algorithm is totally ignored and brute force oversampling is used so that the maximum amplitude sample at the output of the matched filter would be the one to process further. so no tracking, no lock. Ignoring any small quantization error, are there any other issues using this method?


texts: Separating/recovering base signal from two mixed signals, given phase information

texts:
I have collected two signals, $B_1(x)$ and $B_2(x)$. The signals start and end at the same $x$-values. Each signal $B_i(x)$ contains:

a base signal $b(x)$, which is the same for both, and
a signal, either $w(x)$ or $w(x+\phi)$, which are identical except for the shift $\phi$ along $x$.

Signals $b(x)$, $w(x)$, and $w(x+\phi)$ are non-periodic.
Essentially, I have:
$B_1(x) = b(x) + w(x) \\
B_2(x) = b(x) + w(x + \phi)$,
with only access to $B_1(x)$, $B_2(x)$, and $\phi$.
My question: Is is possible to use the phase-shift knowledge from my two signals to recover the base signal $b(x)$? And if so, can you point me in the direction of a method that would work?
I've read about or tried to implement methods such as blind source separation, independent component analysis, and stationary subspace analysis, but none seem to take advantage of the phase information, or work at all for my data.
Note: If it helps, I am able to collect more data with different $\phi$-values, to obtain more $B_i(x)$ signals.
I've attached a figure with toy data (my real data is messier and longer). The top plot shows $b(x)$ (what I'd like to recover), the middle plot shows $w(x)$ and $w(x+\phi)$, and the bottom plot shows the measured signals $B_1(x)$ and $B_2(x)$, representing the combination.


texts: Derivation of Nyquist Frequency and Sampling Theorem

texts: I have been looking through different sites and questions over the internet about Sampling theory, but couldn’t find the clear definition of how nyquist frequency condition is derived? It would be great if someone could direct me to the derivation of the condition that sampling frequency should >= to maximum frequency in the signal


texts: FFT function isn't working on imported signals greater than a second

texts:
I have an fft function in matlab/octave that works great on an imported mono wav/audio signal that is 1 second long, but when the signal imported is more than 1 second the frequency values become incorrect.  To check this I created a 2 second signal using audacity the first half of the 2 second signal (1 second of it) was 500hz at 0.3 amplitude and the second half of the 2 second signal (1 second of it) was 200hz at 0.8 amplitude.  I also checked it with audacity's analyse spectrum option.

When I use the function on this 2 second signal it shows that the max frequency is at 999.95hz when it should show the max frequency is at 200hz.  Can anyone help me get this function to work with signals greater than just 1 second.

I've included the function along with the code that calls the function.

The FUNCTION:

function [freq,amp,ampinv,phase,phase_radtodeg,phase_degtorad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(vp_sig_orig,Fs,phase_goal)

    vp_sig_orig = vp_sig_orig - mean(vp_sig_orig); %remove the mean of your signal due to dc offset

    vp_sig_orig=vp_sig_orig';
    vp_sig_len=length(vp_sig_orig); %get sample rate from vp fs_rate needs to be an even number?

    % Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
    nfft= 2^(nextpow2(length(vp_sig_orig)));

    % Take fft, padding with zeros so that length(fftx) is equal to nfft
    fftx = fft(vp_sig_orig,nfft);
    sigfft= fft(vp_sig_orig);
    sigifft=ifft(sigfft);
    sigphase = unwrap(angle(sigfft')); %get phase of orginal signal

    % Calculate the number of unique points
    NumUniquePts = ceil((nfft+1)/2);

    % FFT is symmetric, throw away second half
    fftx = fftx(1:NumUniquePts);

    % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
    mx = abs(fftx)/length(vp_sig_orig); %replaced for testing from stackexchage

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
    % The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
    if rem(nfft, 2) % odd nfft excludes Nyquist point
    mx(2:end) = mx(2:end)*2;
    else
    mx(2:end -1) = mx(2:end -1)*2;
    end

    %yamp=(mx(1,:)/max(abs(mx(1,:)))*1); % keep at 1, dont use in function use in script file directly
    amp=mx;
    ampinv=abs(amp-max(amp));

    % This is an evenly spaced frequency vector with NumUniquePts points.
    freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft;
    freq=(freq_vect'); %take of round if you want exact numbers

    %get phase of new signal
    phase = unwrap(angle(fftx)); %get phase of orginal signal

    %phase stuff phase adj stuff not fully done yet I don't know the correct way to use it leave at 0 for now
    phase_radtodeg=rtrad2deg(phase); %also check sigphase to see if it's the same.
    phase_degtorad=rtdeg2rad(phase_radtodeg); %for checking
    phase_adjindeg=phase_goal-phase_radtodeg; %used to create new phase to get to goal phase
    phase_adjinrad=rtdeg2rad(phase_adjindeg);

Example: that calls the function

imported_sig_L=imported_sig(:,1)';
[xfreq_orig,yamp_orig,yamp_inv,phase_orig,phase_in_deg,phase_in_rad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(imported_sig_L,fs_rate,0);

array1=[xfreq_orig,yamp_orig,yamp_inv,phase_orig,mod(phase_in_deg,360),phase_in_rad,phase_adjindeg,phase_adjinrad]; %have degs only go up to 360

%Sort Array Section

array1_sort_amp_norm=sortrows(array1,-2); %sort by amplitude
array1_sort_amp_inv=sortrows(array1,-3); %sort by amplitude inverse

%get max freq checking to make sure freq isn't zero if so go to next row
if (array1_sort_amp_norm(1,1)>0); % To get max amplitude
    maxfreq_normA=array1_sort_amp_norm(1,1)
elseif   (array1_sort_amp_norm(1,1)<=0); %
    maxfreq_normA=array1_sort_amp_norm(2,1)
end

if (array1_sort_amp_inv(1,1)>0); % To get inverted max amplitude
    maxfreq_invA=array1_sort_amp_inv(1,1)
elseif   (array1_sort_amp_inv(1,1)<=0); %
    maxfreq_invA=array1_sort_amp_inv(2,1)
end

Thanks


texts: What does Cooley–Tukey FFT algorithm calculates in order to compute DFT

texts:
0
I'm relative new to this subject, I've watched many videos explaining FFT and DFT and read some articles.
I wanted to see how I could implement FFT in C++ and then I found this code, it works but I don't fully understand it, for example, I'm not sure what is he using the nested-for loops for, I thought it's for Matrix multiplication but maybe it's just using the DFT formula for each element of the series, in addition why does he needs to reverse the bits?
if someone could explain in "simple words" the idea behind this code I'll appriciate it because I just started learning it and I've come crossed super complicated stuff which I don't understand.
Here is the FFT Function I'm using:
void fft(Iter_T a, Iter_T b, int log2n)
{
    typedef typename iterator_traits<Iter_T>::value_type complex;
    const complex J(0, 1);
    int n = 1 << log2n;
    for (unsigned int i = 0; i < n; ++i) {
        b[ReverseBits(i, log2n)] = a[i];
    }
    for (int s = 1; s <= log2n; ++s) {
        int m = 1 << s;
        int m2 = m >> 1;
        complex w(1, 0);
        complex wm = exp(-J * (PI / m2));
        for (int j = 0; j < m2; ++j) {
            for (int k = j; k < n; k += m) {
                complex t = w * b[k + m2];
                complex u = b[k];
                b[k] = u + t;
                b[k + m2] = u - t;
            }
            w *= wm;
        }
    }
}


texts: How determine the delay in my signal practically

texts:
I am trying to practically determine filter group delay so that I can try to reduce or adjust for it. I have tried all kinds of methods but still can't get the right result. I know that theoretically to get the group delay you should differentiate the phase response with respect to frequency. How can I use this theoretical result to determine what delay will actually occur to my whole signal envelope? I ask because the theory does not seem to apply in practice in the way I have understood it.

I have tried to align filtered and unfiltered signal to see an actual shift but I do not seem to see any delay.
I have tried adding a spike (Ok which is zero frequency) just so that I can compare the position of the spike in the original and the filtered signal, but the spike does not move in the filtered signal.

Additionally, some filters have both negative and positive theoretical group delay. In those filters should I take the average across all frequencies to determine the effective group delay?

I have added an example below in MATLAB. I used a linear phase FIR filter (a comb filter) which should have a delay of 20 samples according to MATLAB across all frequencies. But I have looked at the plots with the filtered and original signal on the same axis but could not get 20 samples at several frequencies.

%% Filter
b=[1,zeros(1,39),-1];%y(n)=x(n)-x(n-40)
a=1;

subplot(3,1,1)
grpdelay(b,1)

% Simulate
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';

spike=zeros(1,length(t));
spike(300)=0.02;%Feature
spike(150)=0.02;%feature

x1=sin(2*pi*49*t).*wi+spike;
x2=sin(2*pi*25*t).*wi;
x3=sin(2*pi*2*t).*wi;

x=[x1,x2,x3];

y=filter(cb,1,x);
subplot(3,1,2)
plot([x',y'])
title('x vs filtered x');
legend({'x','x-filtered'})

% show spike
subplot(3,1,3)
plot([x',y'])
title('x vs filtered x, zoomed to see spikes');
legend({'x','x-filtered'})
xlim([130,350])

My overall question is how can I practically measure the effective signal envelop delay due to filtering? I would like to match this measured delay with the theoretically calculated delay.


texts:
I know this question has been asked before. It is, however, so confusing, that I'd like to give this another try:
I have come across the following 2 definitions of the DTFT:

So, the first line is the definition given in my lecture. The second line was given in a tutorial on the difference between DTFT and DFT. The first, obvious difference is the limit of the sum. While the first goes from -infinity to +infinity, the second treats a finite signal. I think DTFT are applied to infinite signals. Now, ,the second difference appears in the exponential function.  In the second line, we are using small omega (2*pi*f) multiplied with n, which should be the sample index. In the first line we are using 2*pi*f divided by fs. I think the signal's highest frequency component divided by fs gives the number of samples. So... no, I still don't see the link to n. Can somebody please help me understand this? Although this question appears a lot in this forum?


texts: What is the best method for measuring of room reverberation time?

texts: I want to measure a reverberation time in the room. What is a best way to do that? What kind of signal should be used? I know about balloon and pink noise, but are there any better methods?


texts:
i would a confirm of this: with polyphase structure is possible only design a filter bank with a INTEGER oversampling ratio? For non-integer i've seen the weighted overlap-add metod, is right?

thanks
S.


texts: How to correctly use sinc interpolation in Matlab?

texts:
What is the right way to use sinc interpolation for a given discrete signal $x[n]$? Following is the sinc interpolation formula:
$$x(t) = \sum_{n=-\infty}^\infty x[n] \mathrm{sinc}\left(\frac{t-nT}{T}\right)$$
A simple explanation with a real world example of how and when to use a sinc interpolator will be very helpful.


texts: number of subcarriers for an OFDM system

texts:
For an OFDM system, it is necessary to know the number of subcarriers. For example for the subcarrier $k^e$, we must have: $B_k <B_c$ with $B_c$ the channel coherence bandwith. And $B_k=\frac{1} {NT+I}$ with $I$ is the guard interval, $T$ is the symbol time and $N$ the number of subcarriers.

Here is my questions:

For example, for $N=3$, we have $B_k<B_c$. Is the number of subcarriers small?  Because I always see, N = 64, 1024, ...
Is there a minimum number of subcarriers?


texts: Downsampling in MATLAB - why scaling filter coefficients?

texts:
Trying to undertsand how MATLAB does resampling, looking at toolbox/signal/resample.m.

This is the part where (downsampling) filter coefficients are calculated:

  fc = 1/2/pqmax;
  L = 2*N*pqmax + 1;
  h = firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ;
  h = p*h/sum(h);

(In this case pqmax and p both represent downsampling ratio.)

Not sure what is the purpose of last line where all filter coefficients are scaled by downsampling ratio over coeffs sum, p/sum(h)? Does anyone know the theory behind?

MATLAB firls
MATLAB kaiser
firls theory: Least Squared Error Design of FIR Filters


texts: How to Combine 8 N/8 FFT's into one N FFT

texts:
I need to make in FPGA (using Verilog) an FFT. Input data is N=8192 points at 1 GSPS. However, the FPGA operates at 125 MHz, therefore the data is split into 8 channels (each one at 125 MHz). This splitting of data and computing the N/8 FFT is not problematic for me, it is already done.

What I don't understand is how to combine the outputs of the N/8 FFTs to create one N FFT. I have created a schematic showing the flow of data, where blue blocks represent what is done and green blocks represent what I don't understand.

I understand generally the Cooley-Tukey FFT algorithm and Butterfly diagrams as they relate to 8 point or 16 point data, but I don't understand how these can be expanded to a 8192 input sequence.

Any help regarding the theory, math, or FPGA implementation behind how the green blocks are implemented is greatly appreciated!


texts:
in Digital Image Processing restoration ,from my understanding after reading articles , it's seems that the noise must be dominant in the histogram of the image , which is make it possible , to get the mean and variance of the noise from the histogram it self.
My question is  if I have a little Gaussian noise in the image , which i means that the noise is not dominant in the image , it's seems that it's wrong to extract the mean and variance from the image , because if we do this ,  we finding the mean and variance of the distorted  image it self which is included the noise , yields that we  affect on the original image if we use the model , which is included the model of the noise and the original image together, So how could we avoid this problem ?


texts: Generating a colored waveform of an audio file?

texts:
I would like to ask on what approaches could one attempt to generate a colored waveform out of an audio file. Below are two examples of waveform coloring by commercial products along an audio sample of that music at this position:

Traktor Pro

Pioneer Rekordbox

Audio sample at this position:

https://youtu.be/42t0Sdy-NUw?t=96

By looking at the output, one can infer the following:

low frequency content is colored in red
mid frequency content is colored in green
high frequency content is colored in blue

Obviously, this is a very simplistic analysis of the situation and doesn't really tell how it's been done.

Question:

What kind of algorithm is at work under the hood to produce such colored waveforms?

Disclaimer:

I have asked a similar question a while ago, now I feel it's time to ask it again for these reasons:

question has been asked 6 years ago!
the EchoNest web service used for detecting audio features has closed shop
couldn't really make good use of @pichenettes answer by lack of knowledge
still little public information on the web about this subject
there might be fresh minds that can suggest interesting approaches

(let me know how this question can be improved to be complementary to the one previously asked)


texts: Help or suggestions with Line detection in microscopy images

texts:
I'm doing a personal project in which I want to select some very specific lines in a microscopy image.
The first picture is the raw image, the second picture is the image for the lines (in yellow) I would like to detect. Not straight lines yes, but lines indeed.

Now, I made a very basic code for line detection that kind of "works" detecting lines, but it takes a lot of undesired detected lines in the process. See an example with another image (and the code below it, in Python, using OpenCV).

import cv2 as cv
import numpy as np

#Load and normalize .tif image:
img = 'imagenes.tif'
img = cv.imread(img, 0)
img_scaled = cv.normalize(img,  np.zeros((800, 800)), 0, 255, cv.NORM_MINMAX)

#Resizing images:
pct = 30
w = int(img.shape[1] * pct / 100)
h = int(img.shape[0] * pct / 100)
dim = (w, h)
resized = cv.resize(img_scaled, dim, interpolation = cv.INTER_AREA)

# First stage: Gaussian blur.
gaussian = cv.GaussianBlur(resized,(5, 5),0)
# Second Stage: edge detection with Canny.
low_threshold = 90
high_threshold = 150
edges = cv.Canny(gaussian, low_threshold, high_threshold)
# Third stage: Hough Line transform.
rho = 1
theta = np.pi / 180
threshold = 15
minLineLenght = 15
maxLineGap = 15
myBlankLines = resized.copy()
lines = cv.HoughLinesP(edges, rho, theta, threshold, np.array([]),
                    minLineLenght, maxLineGap)
for line in lines:
    for x1,y1,x2,y2 in line:
        cv.line(myBlankLines,(x1,y1),(x2,y2),(255,0,0),5)

cv.imshow('image_resized', resized)
cv.imshow('Stage 1: Gaussian', gaussian)
cv.imshow('Stage 2: Edges', edges)
cv.imshow('Stage 3: Lines', myBlankLines)
cv.waitKey(0)
cv.destroyAllWindows()

Now. I don't mind building a filter from zero and avoid using OpenCV or any other image processing/computer vision library. But I feel a bit lost on what direction I should go. I'm not sure if:

Doing more pre-filtering.
Adjust parameters for Hough Lines Transform.
Something else related with Machine Learning perhaps.

What would you suggest? Thanks for any heads up you might have on this interesting task.


texts: Output of system defined by differential equation

texts:
I don't fully understand how the output of a system can be derived from the system's differential equation and a given input.

For example:

$$y(0-) = 1 $$
$$y'(0-) = -2$$
$$ u(t) :   \text{Heaviside function} $$
$$ S:  y''(t) +2y'(t) +y(t) = u''(t) -2u'(t) + u(t)$$

Can someone show me a good way to find $y(t)$, given $u(t)$ and $S$?


texts:
I am trying to write a face aligner for as a preprocessing step before feeding my face into a neural network. I have found (using Python) the landmarks of a face via dlib's implementation of Vahid Kazemi and Josephine Sullivan's ensemble of regression trees to predict face landmark estimation. I have a tried a number of different kinds of transforms to no avail. This is my basic pipeline:

Find faces in image and crop image to only having the face :

image = <some_image_with_one_face>
face_detector = dlib.get_frontal_face_detector()
detected_face = face_detector(image)
w,h = 160 # so that facenet can use it
for face_rect in detected_face:

    # First crop the face
    left = face_rect.left()
    top = face_rect.top()
    right = face_rect.right()
    bottom = face_rect.bottom()

    new_face_rect = dlib.rectangle(0, 0, right-left, bottom-top)
    cropped_image = image[top:bottom, left:right, :].copy()

    # Get the the face's pose
    pose_landmarks = face_pose_predictor(cropped_image, new_face_rect)`

Find the landmarks (get_landmark_points returns a 68x2 array

landmarks = get_landmark_points(pose_landmarks, N_LANDMARKS, dlib_point=False)

Define some src and dst points:

top_of_nose = landmarks[27]
left_eye = landmarks[36]
right_eye = landmarks[45]
bottom_lip = landmarks[57]

src = [top_of_nose, left_eye, right_eye, bottom_lip]
dst = [[np.int(0.5 * w), np.int(h/3)],\
       [np.int(0.3 * w), np.int(h / 3)],\
       [np.int(0.7 * w), np.int(h / 3)],\
       [np.int(0.5 * w), np.int(h * (2.0/3))]]

Estimate the transform and also transform the rest of the landmarks

transformed_crop = cv2.warpAffine(cropped_image, transform, (w, h))

# Get the transformed landmarks
transformed_landmarks = np.reshape(cv2.transform(np.expand_dims(landmarks, 1),\
                        transform), (N_LANDMARKS, 2)).astype(np.float32)

# Add the boundary points
transformed_landmarks = np.append(transformed_landmarks, boundary_points, axis = 0)

Here is input and outputs of 2 inputs:

Although it kind of works, it isn't perfect as although the eyes and lip are of the same height in both images, their x position is different. I was wondering if there is some way to improve it, or am I using the wrong technique? It would be helpful if someone pointed me in the right direction, thanks!


texts: I have been working on Digital signal processing since Feb-2013. I have good command over it. Now I want to expand work in field of IMAGE and Video processing. Please suggest me best practical guide books, blogs, online links and other supporting matters related to IMAGE and Video processing .


texts: Choosing sampling frequency and base interval to eliminate leakage

texts:
I'm studying for an exam and have some trouble with an exercise.
So I have a continuous-time signal that is

$$s(t)=\sin(2\pi f_{s1} t) + \sin(2\pi f_{s2}t)$$
with $f_{s1} = 1.4\text{ kHz}$ and $f_{s2} = 1.45\text{ kHz}$.

After sampling the time continuous signal I'll have to get the DFT of the discrete signal. The task is to chose the sampling frequency and the base interval of the discrete signal so that there will be no leakage-effect after calculating the DFT.

I also know that $\frac{N}{\text{period}}$ has to be an integer.

How would you approach something like this? I don't expect to be spoonfed the answer, hints would be nice though.


texts: DFT calculations returning different results

texts:
My coworker ran a DFT query using a software in C downloaded to his PC. His results was the following chart:

I am using a JS tool from corbanbrook/dsp.js

My results were very different.

any ideas why, or where to look?


texts: Detect FSK high/low bits in presence of noise with only simple processing

texts:
I'm currently in the process of writing DSP code, based on the amazing RTLSDR project, and using this as a vehicle to teach myself the basics of DSP (I have an extensive background in software development, but none whatsoever in physics, EE or DSP.)

Specifically, I am trying to implement FSK from scratch using only the most basic building blocks of DSP, and writing all the code by myself. Up until now I have managed to do the following:

Extract a stream of I/Q samples via the RTLSDR interface
Calculate the angle of each I/Q pair on the unit circle by calculating atan2(q, i)
Calculate the phase by subtracting each pair of two consecutive I/Q samples atan2(q_2, i_2) - atan2(q_1, i_1)

So now I have a stream of phases, and I can basically "see" the signal going round and round on the unit circle with values between [0-2*PI). If we look at FSK, this stream will be moving slightly faster whenever a high bit is transmitted, and slightly slower on the low bits. This is where I'm stuck.

The implementation I'm testing right now basically counts how many samples have passed since the last time a full circle has completed, but this seems very prone to errors due to noise.

How do I calculate, from this stream, whether I am currently on a high bit, or a low bit? Since I am just teaching myself here, I'm looking for the most basic method for this step, no complicated algorithms, just yet. Thanks!


texts: On the simplification using trigonometric functions

texts:
Assume I have a matrix $D$ whose its entries are as below :

Where $A$ and $B$ can be written using using the trigonometric functions for (1) as:

My question, Is it possible to simplify (1)
 more? how can we do that ?


texts:
I am studying channel modeling, at the moment, in the book: "MIMO-OFDM wireless communications with MATLAB", the explanation about the model is as follow:

The code for generation of Doppler spectrum is
% fd = Maximum Doppler frequency
% Nfft= Number of frequency domain points
df = 2*fd/Nfft; % frequency spacing
% DC component first
f(1) = 0;
y(1) = 1.5/(pi*fd);
% The other components for one side of the spectrum
for i = 2:Nfft/2,
f(i)=(i-1)*df; % frequency indices for polynomial fitting
y([i Nfft-i+2]) = 1.5/(pi*fd*sqrt(1-(f(i)/fd)^2));
end
% Nyquist frequency applied polynomial fitting using last 3 samples
nFitPoints=3 ;
kk=[Nfft/2-nFitPoints:Nfft/2];
polyFreq = polyfit(f(kk),y(kk),nFitPoints);
y((Nfft/2)+1) = polyval(polyFreq,f(Nfft/2)+df);

I cannot understand the code y([i Nfft-i+2]) = 1.5/(pi*fd*sqrt(1-(f(i)/fd)^2));, it seems like generating the Dopple spectrum, but I cannot find any similar formula for this line of code
What is the purpose of the last two line of code?

polyFreq = polyfit(f(kk),y(kk),nFitPoints);
y((Nfft/2)+1) = polyval(polyFreq,f(Nfft/2)+df);


texts: Thresholding algorithm for separating skin from clothes in python

texts:
I have an several images that look like the following:

I would like to separate the skin from the trousers here without specifically training a segmentation model. Can anyone suggest any such algorithm that can potentially achieve this task? Also a dummy example that can demonstrate this would be great. Thank you


texts:
I have a rotating shaft where I need to move from the time domain to the angular domain to look for vibration faults. I have the data from a tachometer.

From another post I found these steps:

  The steps involved for time domain --> Angle domain are:

  Use the tachometer signal to obtain a speed signal. This is the
  inverse of the difference in time for each trigger point.
  Depending
  upon the variability of the tachometer signal least squares cubic
  spline fitting is often used. Matlab spline library will work as will
  matlab central libraries. I use the fastBspline library.
  Integrate the
  speed signal to obtain shaft angle position.
  Determine the equal angle
  increments that you want to use. Keep the Nyquist theorum in mind, in
  other words there should be roughly the same number of points in the
  angle domain as in the time domain.
  Switch the (time, angle) data to
  (angle, time). Resample the (angle, time) data to obtain uniform angle
  increments. Once this is done, do all the same processing you would
  have done, but now you are in the angle order domain.

Questions

In regards to step 2: What is returned if I integrate the speed signal (which is an array)? Surely that would not be a function or an array?

In regards to step 5: Why is it necessary to change from (time, angle) to (angle,time)? How is this plot related to the original vibration signal (time, acceleration)?


texts: Does MATLAB code retain its accuracy when is executed in dsp kit?

texts:
Ladies, Gantlemen,
Please let me ask your help on my problem. MATLAB code (at https://metrw-pitch.blogspot.com/ you can see code in .mlx (for run it you want delete comments) and .m format) executing in MATLAB online R2020b, outputs good accuracy even for low fundamental frequencies (fundFreq) but code (originally in my answer to a question of mine) of similar algorithm in C++ (acctually I have translated this code into Matlab's language), executing in Visual Studio 2019, outputs bad accuracy in low (lower than 150 Hz) frequencies.
Question is whether MATLAB's code executing in dsp kit retains its accuracy. I'm homeless (in France) and get internet access in Social Services and Post Offices like this I'm just now. So it's impossible for me to run code in dsp kit. You are kindly requested test MATLAB's code in your dsp kits and say me its output for fundFreq 50, 60, 70 Hz. Be aware when change FundFreq, change, if necessary, FFTfundFreq and/or GridDemiSpan so that grid's span covers FundFreq.
I have asked similar question on 8 Feb 2021 at 13:25, at MATLAB but nobody answered yet.
Copy of Matlab's code without comments is following.
format long;
SampFreq = 16000;
Segm = 1:1600;
FundFreq = 50;
FFTfundFreq = 41;
GridDemiSpan = 10;

FirstHarmAngles = FundFreq*2*pi/SampFreq*Segm+1.9*pi;
SinFirstHarmAngles = sin(FirstHarmAngles);
SecondHarmAngles = FundFreq*2*2*pi/SampFreq*Segm+0.9*pi;
SinSecondHarmAngles = sin(SecondHarmAngles);
ThirdHarmAngles = FundFreq*3*2*pi/SampFreq*Segm+0.3*pi;
SinThirdHarmAngles = sin(ThirdHarmAngles);
Xn = 170000*SinFirstHarmAngles+220000*SinSecondHarmAngles+...
150000*SinThirdHarmAngles;

Freqs = FFTfundFreq-GridDemiSpan:0.1:FFTfundFreq+GridDemiSpan;
MagnSqrd = ones(1,201);

for f = 1:201
   Angles = Freqs(f)*2*pi/SampFreq*Segm;
   XnCos = sum(Xn.*cos(Angles));
   XnSin = sum(Xn.*-sin(Angles));
   MagnSqrd(f) = XnCos.^2+XnSin.^2;
end

[maxMagnSqrd, maxMagnSqrdIndex] = max(MagnSqrd);
GRIDfundFreq = Freqs(maxMagnSqrdIndex);
disp(GRIDfundFreq);

With regards and friendship, Georges Theodosiou


texts:
I am working on a matlab program which needs me to detect the pattern of surges in a signal and then remove them. I tried out various filters like median, hampel, butterworth etc. but couldn't get an effective output.
Here's a snap of my signal :

Now, I need to filter out the surges. Could anyone please suggest me a method on how to approach this problem? Sampling frequency is 4khz and No. of data points is 40k.


texts: the DFT of a periodic signal represented by a fourier series

texts:
If I have a signal represented by a Fourier series(like in the photo), which is sampled with $T_s$:

$$x[n]=x(t=nT_s) = \sum_{m=-\infty}^{\infty}a[m]e^{j2\pi m(nT_s)/T_0} $$

How do I find its DFT? I can't understand that.
Sorry for the English mistakes.


texts: Design one of 4 type of FIR filter by window function?

texts:
There're 4 type of FIR to chose from

symetry ,even
symetry ,odd
antisymetry, even
antysymetry, odd

Also, there are many type of window like rectangular, Hanning, Kaiser,etc So how to use a window function to design the N-length filter and what's the point of doing it? That's what my teacher ask me to do but I don't understand the reason

What I confuse here is that for example I want the ideal filter $H(\omega)$ that have passband in $\omega_a <\omega<\omega_b$ then what I need to do is just use the inverse Fourier transform to the frequency response of the ideal filter to find $h[n]$ that is infinite then use the window function to make it finite. Then what's the point of make it look like the form of one of those 4 type? Moreover, what I need to do after doing inverse Fourier Transform and window function to make it in the form of FIR digital filter with N-length?


texts: Least Squares Re Synthesis of Short Time Fourier Transform (STFT)

texts:
Working with my own MATLAB implementation of the short-time Fourier transform (STFT), I've managed to write code for the analysis step where a 1D time-domain signal $s[t]$ is progressively windowed, taken into the Fourier domain and arranged in a 2D matrix $S[t,\omega]$.  Each element in the 2D matrix is a function of time $t$ and frequency $\omega$.

This presentation is a rather nice overview of the STFT, and gives a number of equations detailing the analysis and re-synthesis steps.

However, I would like to be able to arbitrarily modify $S[t,\omega]$, and then use the re-synthesis to get back $s[t]$.

I believe that I should be able to change $S[t,\omega]$ in whatever way that I want, and then obtain $s[t]$ by the re-synthesis procedure.  This seems to be very similar to the idea of the phase vocoder.

As noted in section 3.1 of the presentation, the time-domain signal $s[t]$ can be recomposed using a least-squares procedure.  This is given as Equation (6) in the 1984 paper by Griffin and Lim.

The least-squares procedure is required to be applied when $S[t,\omega]$ is modified in some way.

Question:

What does Equation (6) of the Griffin and Lim paper mean?
What steps do I follow to numerically implement Equation (6)?

In the presentation, the equation is written in a slightly different way:

$$x(n)=\frac{\sum_{l=-\infty}^\infty w\left(n-lI\right)y\left(lI, n\right)}{\sum_{l=-\infty}^\infty w\left(n-lI\right)^2}$$

Note that $x(n)$ is the re-synthesized time-domain sequence, $w(n)$ is the window function, and $y(n)$ is the time-domain version of a column of the 2D matrix.

Steps:

From the presentation, here are the steps that I think are required to do the re-synthesis:

Let w_n be the discrete window vector and y_n(:,k) be the time domain vector computed using the IFFT on a column k of the 2D matrix.  Both w_n and y_n(:,k) are the same length.
Then, using Matlab syntax, we compute the point-by-point multiplication:
w_n .* y_n(:,k)

Is this the numerator of the expression above?
What happens during steps 3 and 4?
What do the infinite summations signify?


texts:
Is there a way to measure each pixel's content in an image somehow?

Like can I interpret each color as some kind of a single value, rather than a triplet?

Perhaps I can think of them as $\mathbb{R}^3$ vectors and then e.g. take their norms?


texts: How to determine ALPHA, smoothing constant, of a LPF?

texts:
I am implementing a simple LPF in Java.

But I have problems in choosing the value of ALPHA.

public class LowPassFilter {

    /*
     * Time smoothing constant for low-pass filter 0 ≤ α ≤ 1 ; a smaller value
     * basically means more smoothing See:
     * http://en.wikipedia.org/wiki/Low-pass_filter#Discrete-time_realization
     */
    private static final float ALPHA = 0.2f;

    private LowPassFilter() {
    }

    /**
     * Filter the given input against the previous values and return a low-pass
     * filtered result.
     *
     * @param input
     *            float array to smooth.
     * @param prev
     *            float array representing the previous values.
     * @return float array smoothed with a low-pass filter.
     */
    public static float[] filter(float[] input, float[] prev) {

        if (input == null || prev == null)
            throw new NullPointerException("input and prev float arrays must be non-NULL");
        if (input.length != prev.length)
            throw new IllegalArgumentException("input and prev must be the same length");

        for (int i = 0; i < input.length; i++) {
            prev[i] = prev[i] + ALPHA * (input[i] - prev[i]);
        }

        return prev;
    }
}

ALPHA is defined to be dT/(dT+RC), where dT is the event delivering rate and RC is the time-constant of LPF.

To find the value of ALPHA, I need to know the values of dT and RC.

My LPF should have a cut-off frequency of 4Hz and the sampling rate is 350Hz.
Anybody can help determine the proper ALPHA value?

Thanks in advance!


texts:
Disclaimer: I posted this question on StackOverFlow but learned that it is better suited here.
I am a bit new to image processing so I'd like to ask you about finding the optimal solution for my problem, not help for code.
I couldn't think of a good idea yet so wanted to ask your advice.
I'm working on a project using OpenCV on counting the vehicles from a video file or a live camera. Other people working on such a project generally track the moving objects then count them but instead of that, I wanted to work with a different viewpoint; asking user to set a ROI (Region of interest) on the video window and work only for this region (for some reasons, like to not deal with the whole frame and some performance increase), as seen below. (The user can set more than one ROI and the user is asked to set the height of the ROI about two times of a normal car by sense of proportion.)

I've done some basic progress so far, like background updating, morphological filters, thresholding and getting the moving object as a binary image something like below.

After doing them, I tried to count the white pixels of the final thresholded foreground frame and estimate whether it was a car or not by checking the total white pixels number (I set a lower bound by a static calculation by knowing the height of ROI). To illustrate, I drew a sample graphic:

As you can see, it was easy to calculate the white pixels and check if it draws a curve by the time and determining whether a car or something like noise.
Problem
I was quite successful until two cars passed through my ROI together at the same time. My algorithm incorrectly counted them as one car. I tried different approaches for this problem and similar to this like long vehicles but I haven't found a satisfactory solution.
Question
Is it possible to handle this task by this approach of pixel value counting? If not, what do you suggest?


texts: What are the units of the product of two signals?

texts:
In a standard mixer, we multiply two signals $x_1(t)$ and $x_2(t)$, each with units of volts $(V)$, and get at the output:

$$y(t)=x_1(t)\,\cdot\,x_2(t)$$

I have always thought that the output of mixing operation was another voltage signal. But how can this be? Looking at the units of the factors, we multiply a quantity of volts, $x_1(t)$, by another quantity of volts, $x_2(t)$. So shouldn't the units of $y(t)$ be $V\cdot V=V^2=W$?


texts: Despeckling Ultrasound images using Diffusion Filtering

texts:
I am working on Despeckling Ultrasound images. I have been trying out some basic filters and their hybrid combination.. Now I would like to try Diffusion Filtering techniques to despeckle the ultrasound images. I don't know anything about the diffusion filtering techniques (Isotropic diffusion or Perona Malik Anisotropic Diffusion Filter). I would like to learn from the basic. Please suggest some books to learn about Diffusion techniques.. Please keep in mind that I am not an advanced learner.

I did google and found this book (Despeckle Filtering Algorithms and Software for Ultrasound Imaging By Christos P. Loizou, Constantinos S. Pattichis, Costas Pattichis) very useful. But I couldn't able to buy it. If anyone have the PDF of this book (even 1 or 2 chapters), please provide me the link. Also provide some useful links to learn these techniques


texts: Nyquist Plot for transfer functions with poles at the origin

texts:
I'm learning Nyquist plots and something has been seriously bugging me when treating poles or zeros in the origin. Nyquist plots obtains information based on the argument principle which states

"If f(z) is a meromorphic function inside and on some closed contour C, and f has no zeros or poles on C, then

$${\displaystyle \oint _{C}{f'(z) \over f(z)}\,dz=2\pi i(N-P)} \oint_{C} {f'(z) \over f(z)}\, dz=2\pi i (N-P)$$
where N and P denote respectively the number of zeros and poles of f(z) inside the contour C, with each zero and pole counted as many times as its multiplicity and order, respectively, indicate."

So we just ignore the fact the transfer function contains a pole over the contour $C$? How can we treat them normally when it clearly violates the argument principle?


texts: What does it mean by DN in salt-pepper noise?

texts: What is DN on page 245 of these notes shown in the slide below?


texts:
In the above image I have created a time series (in cyan) that is a simple addition of a major cycle (in green), a minor cycle (in cyan) and a linear trend (not shown separately). What I would like to do is filter this time series to remove the minor cycle so what is left would be the trend plus major cycle (in red). In practice no information would be available about the minor cycle, although the period of the major cycle would be known from applying code for the Hilbert transform.

One way in which I might do this is:

apply a bandpass filter set to the major cycle period, on the original data
subtract this bandpass from the original series
apply the Hilbert transform to the remainder to get period of minor cycle
apply another bandpass filter set to this minor cycle period, on the original data
subtract this second, minor cycle bandpass from the original data
and hopefully end up with the desired red output

First off, would this be a valid approach, and secondly is there another way in which I could filter out the minor trend?


texts: Does this make sense for cross-correlation in frequency domain?

texts:
If I have a signal A at 2 kHz and signal B mixed at 2 kHz + 3kHz, and then i take the FFT of them both, and then cross-correlate them many times, the 3kHz peak should disappear? Since this peak is not present in signal A.

Thanks


texts: Tikhonov Regularization Alternate Formulation

texts:
2.2 Tikhonov Reqularization

  Tikhonov regularization, named for Russian mathematician Andrey Tikhonov, attempts to fix the issue that arises when the least squares method is used with an ill-posed inverse problem by adding an additional constraint on the minimization:
  $$\min_x\{\Vert Ax-b \Vert ^2_2 + \lambda^2\Vert x\Vert^2_2\}$$
  In the new minimization, the additional constraint is the squared norm of $x.$ $\lambda$ is a non-negative constant decided in advance, acting as a weight on the strength of the additional constraint. Some alternate ways to write Tikhonov regularization include
  $$\bbox[lightgreen]{(A^T A + \lambda^2 I)x=A^{T}b}$$
  and
  $$\min \left\Vert\genfrac{[}{]}{0pt}{}{A}{\lambda I}x - \genfrac{[}{]}{0pt}{}{b}{0}\right\Vert$$

Can someone explain the alternate formulation in green. I am not able to convert this formulation to standard one. I get one Atransposedinverse extra. The second formulation can be easily converted to original formulation, though. Thanks.


texts: Finding a better threshold in a more accurate manner in OpenCV

texts:
I have a grayscale image,
the dominant colors in the background are shades of white and light gray, I have tried using otsu binarization but it gave no good results, so I decided according to the histogram below to write manually:
img_gray[img_gray > 130] = 255
so all white/gray shades have been converted to white. it worked impressively well, I was wondering is there an efficient way getting this threshold without re-inventing the wheel, I know I can apply: np.bincount(img_gray.ravel()) and try to find the beginning of the ramp before the peak.

Important fact to say is:
The images' characteristics are there will be always two peaks, the first will reside between 0 to 50 (it's the small one) and the other somewhere in the light white shades.
Is there an approach to find this threshold?
Thank you all


texts: Elimination of image blurring using noise model restoration techniques

texts:
How can the blurring be eliminated from image using the noise model restoration process techniques?

To my knowledge the answer must be filter.If yes,then which filter is used?


texts:
I have position $\theta(s) = \frac{s^2+1}{s^4+2s^3+3s^2+4s+5}$ and velocity $\dot\theta(s) = \frac{s(s^2+1)}{s^4+2s^3+3s^2+4s+5}$.

How do I construct 2-state-observer with sampling time $T_s$?

$\left[
\begin{array}{c}
  \theta_{k+1}\\
  \dot\theta_{k+1}\\
\end{array}
\right]
=A\left[
\begin{array}{c}
  \theta_{k}\\
  \dot\theta_{k}\\
\end{array}
\right]$


texts: How to interpret results of Fourier transform by using Python

texts:
Connected with this topic , I would like to have some clarifications about the Fourier Transform and the Python tools that make it.
Having my data set whose plot is showed in the above link, my array of data is stored in rate. Then, if I check the result from the Fourier transform:

yfft = np.fft.rfft(rate)
y_smooth = np.fft.irfft(yfft)
print y_smooth

I get:

[  3.20000028e+01+0.j          -7.92633533e+00-7.75032655j
-5.34245688e+00-3.82390808j  -2.39880935e+00-0.60037314j
-8.71249166e-01-0.04362804j  -8.52507194e-02+1.09973189j
-1.58332047e-01+1.15256295j  -3.38909974e-01+0.96327143j
 5.51323015e-02+0.43936749j   3.81205886e-01-0.02489853j
 2.99110153e-01-0.26319101j   3.28283652e-01+0.11668814j
 1.08849205e-02-0.11317937j   3.79952011e-01+0.22278952j
 7.37720998e-02-0.26810083j  -2.14247586e-01-0.25950865j
-1.32242152e-01+0.j        ]

So, I am quite confused about this result, and I would need some tips how to interpret these numbers.

1) The array contains $17$ complex number : does it mean that my signal can be approximated by a superimposition of $17$ harmonics?

2) Are these numbers the power of each harmonic, or what?

3) are they ordered by the frequency order, or by power? Let's say I want to put to zero the higher orders, which ones should I put to zero?


texts: Pairwise Error Probability simulation versus SNR for Space-time code

texts:
I have a MIMO system using Space-time coding. I am required to produce simulation results in the form $\text{PEP}=f(\text{SNR})$, meaning a curve representing the Pairwise Error Probability (PEP) for each SNR value (say from 0 to 25 dB).

The PEP is the probability of mistaking a matrix $M_a$ with a matrix $M_b$ which is $P(M_a \rightarrow M_b)$.

Given that I have a matrix constellation in my system of cardinality 16 (16 matrices in my code book) how am I to represent $\text{PEP}=f(\text{SNR})$?

Should represent the result of just one pair of matrices ? (exemple $P(M_1 \rightarrow M_2)$ would this be fair given that I have other pairs of matrices?

Or is there a method to group the PEPs of all the pairs into one curve ?

Attached is an example of the curve I am looking to produce.
thanks


texts: Fast way to calculate a few points in convolution?

texts:
Say I have two 1000-point sequences. I can compute the entire convolution fairly efficiently using FFTs (keeping in mind that zeroes have to be padded on to avoid doing a circular convolution).

However, I'm only interested in about 500 points in this convolution. To be more specific, say the indices of the convolution are 0 to 1998; I'm only interested in 1000 - 1499 inclusive.

Is there some way to calculate only these points, without calcualting the entire convolution and throwing most of it away?


texts:
Possible Duplicate:
  Real Discrete Fourier Transform
  What is the most lucid, intuitive explanation for the various FTs - CFT, DFT, DTFT and the Fourier Series?
  Discrete-time Fourier transform

I read about the Discrete Fourier Transform. I understand how to compute it given a sequence of numbers(digital signal), but I still didn't get what's the meaning of the new sequence I get as the output.

I mean, I know that the Fourier Transform takes an analog signal representation in the time domain and returns its representation in the frequency domain, its spectrum.
If I look at the analog signal in a time-window $T$, sampling it $N$ times , then the sampling interval will be $t_0 = T/N$ and the sampling frequency will be $f_s =1/t_0$.
The sampled values of the analog signal give me a new digital signal, lets say $s_n$

So what is the meaning of the DFT of this digital signal $s_n$? How the DFT sequence related to $s_n$? How it is related to the original analog signal? What is the meaning of these frequency components I get? Should they be taught of an approximation the previous analog signal frequency components or as sampled values of its spectrum?

Thanks a lot!


texts: Plotted ECG signals are not around Amplitude 0 line

texts:
For my project I have to have few datasets of ECG signals and I did get them from Physionet website.

What I Got

So I downloaded few ECG signals and plotted them in Matlab. Below shows one plot (I lead) I got.

As you can see, this signal is not in the Amplitude 0 line. But all the ECG recordings I have seen and saw in research papers are horizontal. Like in below image.

My Question

I'm new to ECG signal processing. So I'm curious to know,

1)Whether datasets having this kind of abnormal behavior of ECG are wrong?

2)Is there any Matlab function which can convert such data to view like in 2nd image I attached?

UPDATE

To download the dataset,

https://www.physionet.org/cgi-bin/atm/ATM

Go to the above link and select "PTB Diagnosed ECG database" as the Database and select a patient record (E.g: s0010_re).


texts:
What would happen if I am using the maximum frequency as the sample rate for sampling a pure tone sine wave?
For example, a $10\ \rm kHz$ sampling frequency for a $10\ \rm kHz$ monotone sine wave. What effect would aliasing produce at the output?


texts:
can we say there is no sound produced when the amplitude is 0?
if yes what can we say when the amplitude is negative?


texts: Relating Fourier and power-law descriptions of noise

texts:
I'm trying to understand this PDF, which is in English but has Polish authors, so there may be translation issues confusing me. The subject is the frequencies of sentence lengths in novels.

Eq. (1), $S(f)=\left|\sum_{j=1}^{j_\max}l(j)e^{-2\pi ifj}\right|^2$, claims to present a power spectrum "of the series $l(j)$ representing lengths of the consecutive sentences $j$". It is then claimed empirical $S$-against-$f$ plots exhibit Eq. (7), $S(f)=1/f^\beta$.

I expect that might be meant to be a proportionality relation, but let's put that aside for the moment. To my mind, these equations are incompatible, because (1) implies $S(f+1)=S(f)$, with a finite peak in $S$ at integer $f$, including $f=0$, by the triangle inequality in $\Bbb C$. By contrast, a $1/f^\beta$ law is monotonic decreasing on $f\ge0$, with a divergence at $f=0$.

Figure 1 suggests values of $f$ from about $10^{-3}$ to $10^{-1}$ are considered, which on log-log axes allows one to make a power-law approximation, but I don't see the merit in that. But supposing this is a standard procedure, how would I estimate $\beta$ for $l$ from my own time series? Is there, for example, a computationally cheap library function to estimate $\beta$?


texts: Why could window function help to get more accurate specific frequency amplitude?

texts:
In order to calculate amplitude of sinusoidal signal at a frequency around 60 Hz (but may be not 60 Hz), window function is adopted while using DFT.

hk is the coefficient of window function and Kdc is the sum of these coefficients.
Np is the number of sample in one period when frequency is 60 Hz.
(XL1, XL2, XL3 are time series sample of three phase electrical quantities.
Sqrt(2) is due to some electrical reason.)
As far as I know, window function is helpful to reduce discontinuity when frequency deviates a little bit. However, I am really confused why this could help to get a more accurate sinusoidal signal amplitude.
I tried to derive, but after getting the steps below, I don't know how to prove V1 is closer to A than V2.
h1 and h2 are two types of window function. h1 corresponds to rectangular window and h2 represents triangle window.


texts: I want to calculate the bit period in my RF packets that I transmit them with 50Kbaud (Symbol rate=bit rate in my case) , how can I calculate the bit period?