texts: Kalman Filter - Implementation, Parameters and Tuning

texts:
First of all, this is the first time I try to make a Kalman filter.

I earlier posted the follwoing question Filter out noise and variations from speed values on StackOverflow which describes the background for this post. This is a typical sample of values I'm trying to filter. They don't necessarily have to decrease which is the case here. But the rate of change is typically like this

X------- Y
16 ---233.75
24 ---234.01
26 ---234.33
32 ---234.12
36 ---233.85
39 ---233.42
47 ---233.69
52 ---233.68
55 ---233.76
60 ---232.97
66 ---233.31
72 ---233.99

I have implemented my Kalman Filter according to this tutorial: Kalman Filter for Dummies.

My implementation looks like this (pseudocode).

//Standard deviation is 0.05. Used in calculation of Kalman gain

void updateAngle(double lastAngle){
  if(firsTimeRunning==true)
     priorEstimate = 0;               //estimate is the old one here
     priorErrorVariance = 1.2;        //errorCovariance is the old one
  else
     priorEstimate = estimate;              //estimate is the old one here
     priorErrorVariance = errorCovariance;  //errorCovariance is the old one
  rawValue = lastAngle;          //lastAngle is the newest Y-value recieved
  kalmanGain = priorErrorVariance / (priorErrVariance + 0.05);
  estimate = priorEstimate + (kalmanGain * (rawValue - priorEstimate));
  errorCovariance = (1 - kalmanGain) * priorErrVariance;
  angle = estimate;              //angle is the variable I want to update
}                                //which will be lastAngle next time

I start out with a prior estimate of 0. This seems to work fine. But what i notice is that the kalmanGain will decrease each time this update is run, which means that I trust my new values less the longer my filter is running (?).
I don't want that.

I went from just using a moving average (simple and exponential weighted) to using this. Right now I can't even get as good results as that did.

My question is if this is the right implementation and if my prior error variance and standard deviation looks good according to the sample values I have posted? My parameters are actually just picked randomly to see if I could get some good results. I have tried several different ranges but with poor results. If you have any suggestions to changes I can do, it would be really appreciated. I'm sorry if there is some obvious things missing. First time posting here too.


texts: Most efficient phase linear x4 resampling algorithm

texts:
I'm working on a program in which I need to do the following about 20 times for each block of audio:
- Upsample from 192 to 768 kHz (x4)
- Clip
- Downsample back to 192 kHz (/4)

What I'm using now:

FFT block size 4096
Add 0's at the end of the FFT block
IFFT block size 16384
Clip audio
FFT block size 16384
IFFT block size 4096 (this ignores the higher freqs so it automatically lowpass filters)

For the FFT's I'm using Intel's IPP library, and all the other code (adding 0's, clipping) is written using AVX intrinsics.

Still, doing this 20 times is heavy, and I know that it is possible to do it more efficiently, because in the 1st FFT step I'm adding 0's which are then used in the calculations (while they could be discarded) and at the other end, 75% of the FFT values that are calculated are discarded as well. So, even with FFT's, this should be possible more efficiently, which leads me to think that there must be more efficient resampling algorithms out there. But I haven't been able to find one yet.


texts: DFT (FFT) of a Real Even Function Doesn't Yield Real Only DFT Signal

texts:
There are a lot of similar questions already, but I can't find one answering this problem.

I used the following MATLAB code to perform the FFT on a Lorentzian input signal

x = linspace(-10,10,2^10);
y = 1./(x.^2 + 1);
F = fft(y);
IF = imag(F).^2;
power = trapz(IF); %integrating the imaginary part squared

Now, from everything I've read, the FFT of a real, even function is supposed to give $0$ imaginary part. But in this case, I get

power = 3.86

which doesn't seem to be just a small numerical inaccuracy. What's going on here?

Attached is the imaginary spectrum after fftshift.


texts: temporal high-pass filtering to approximate phase derivative

texts:
I am currently following the algorithm as based on a previously published research article "A Motion Detection Algorithm Using Local Phase Information" by Lazar et al. The basic overview of the algorithm is as follows:

I am able to do everything upto taking the 2D FFT. Once I take the FFT of a block of size 32 x 32 pixels, (Gaussian window has already been applied to this block), I get an array with the contents as follows:

For every entry in the 32 x 32 array, I understand that you can extract phase using the angle() function of MATLAB. However, I am not quite understanding

what the frequency(wx, wy) is referring to
how to temporally high pass filter phase information to approximate the change in local phase information.

Thank you.


texts: Plotting Hilbert Spectrum from Hilbert Huang Transform

texts:
Before anyone calls it a repeat post from past posts. I'd like to point out I've read all of those; however, I have questions on making a hard coded version of the hilbert spectrum.

I do indeed have all the ingredients aka the hilbert transform and empirical mode decompositions.

Here is my code below:

def hilbert_huang(self, imf, res):
    # Reconstruct the signal
    sig = np.sum(imf, axis=0) + res
    # Hilbert Transform of IMFs
    hil_imf = np.array([sp.fftpack.hilbert(imf[i]) for i in range(len(imf))])
    # Instantaneous Amplitude
    inst_amp = np.array([np.sqrt(imf[i]**2 + hil_imf[i]**2) for i in range(len(imf))])
    # Instaneous Phase
    inst_phase = np.array([np.arctan(hil_imf[i]/imf[i]) for i in range(len(imf))])
    # Instaneous Temporal Frequency
    inst_freq = np.array([(1/(2*np.pi))* np.gradient(inst_phase[i]) for i in range(len(imf))])

    return inst_amp, inst_phase, inst_freq, sig

However, I'm not sure how to construct the spectrum from all this, I know that the formula is given by
$$ H(\omega, t) = \Re\left(\sum_{i=1}^{n} a_{i}(t) e^{j \int \omega_{i} (t) dt}\right) $$

But I'm not sure how to fix $\omega$ , given that $\omega$ is an array with values , in the wikipedia article it says to let $\omega = \omega_{j}$ where $j$ indicates which IMF it's from but again $\omega$ is a vector, so do I just find the max? Let me know if I'm not clear or need to provide any more mathematical expressions or code.


texts: Contour and area, raw (spatial) and central image moments

texts:
I recently started using image moments for image processing of binary images. I read that the $0^{th}$ order contour moment is the perimeter and the $0^{th}$ order area moment is the area. These raw moments are both given by:

$M_{ij} = \sum_{x}\sum_{y}x^iy^j$.

This means If I have an image like this (but binary, foreground pixels shown in blue), the $0^{th}$ moment will correspond to the the perimeter, since it's an image of a contour :

On the other hand, if I have an image like this (foreground shown as while), I will get the area of the object as the $0^{th}$ moment:

Since I want to use the contours to get more properties, I also calculate the higher order ($1^{st}$, $2^{nd}$, $3^{rd}$ order) raw contour moment. I want to use these to get the central moments.

Formulas I am using to get the central moments are:

$\mu_{00} = M_{00}$

$\mu_{01} = 0$

$\mu_{10} = 0$

$\mu_{11} = \frac{M_{11}}{M_{00}} - x_c*y_c = \frac{M_{11}}{M_{00}} - (\frac{M_{10}}{M_{00}})*(\frac{M_{01}}{M_{00}})$

The formulas for calculating central moments are using raw moments. My question is: Which raw moments are used to calculate central moments, area or contour?. My guess is area moments, since the $0^{th}$ order central moment is also equal to the area, which is in fact the $0^{th}$ order area moment.

Additionally, can I calculate the central moments based on contour raw moments?


texts: Are reconstruction filter always needed?

texts:
in mixed (digital and analogue application) where you have to convert a signal from continuous to digital time domain and then back again to continuous time an anti-aliasing LPF is needed and in general also a LP reconstruction filter.

Are there some applications (if so what are) in which the reconstruction filter is not needed or not strictly necessary?


texts:
I have two signals, one whisper data and one normal voice data of the same sentence.

I'm trying to build a voice conversion, but just want to start off with some simple signal manipulations.

I wanted to get the excitation of several vowels first to apply the excitation to the whispered data; basically to get a Electrolarynx kinda robotic sound. Therefore, I splitted the signal into smaller signals with only vowels. So now I have, e.g. an 'i' signal in whispered data and normal voice data.

How can I get the excitation now? These are some buzz words I heard concerning this topic:
- calculating the fft of both signals
- inverse fft

I'm pretty clueless here, so any help/tip is appreciated.


texts: What is the Compression Algorithm used by Facebook?

texts:
I would like to know what is the Compression Algorithm used by Facebook for compressing the images? I uploaded a image of size 242 Kb and I downloaded the same image from my Facebook album. The size of downloaded image is 46 Kb and the subjective quality is similar, I wasn't able to differentiate between the original (uploaded) and downloaded (from Facebook) image. The downloaded image was downsampled though, from 1226 x 1600 to 552 x 720. So Facebook compresses the image at a good ratio while we are uploading it to our Facebook account without affecting the subjective quality much.

I wonder what compression algorithm they are using for this purpose. What would be your best guess?

                           Original image: tamil_girl_painting (size 1226 x 1600)

                  Downloaded (from Facebook) image: tamil_girl_painting (size 552 x 720)


texts: Calculating amplitude, time-shift and phase of third sinusoid as a sum of two others

texts:
My apologies in advance, I am new to this digital signal processing. Question is at the bottom

I have two 4000hz sinusoids given by the following formulas:  unknown frequency such that:
$$x_1(t) = 24\cos(2\pi4000(t-t_{m1})) \\ x_2(t) = 28.8 \cos(2\pi4000(t-t_{m2})) $$
and a third given by
$$ x_3(t) = x_1(t) + x_2(t)$$
where
$$t_{m1} = 0.000775 \\ t_{m2} = -0.00258 $$

I have calculated the phase of the first two sinusoids (in radians) as $$ \phi_1 = 0.329 \\ \phi_2 = -1.132$$
with the formula $$ \phi_i = \omega_if_i $$ where
$$\omega_i = 2\pi f $$

Finally, I need to calculate the amplitude $A_3$ and phase $\phi_3$ of $x_3(t)$.
I am unsure which formula to use to do so. I believe $A_3$ is calculated as: $$a_3 = \sqrt{A_1^2 + A_2^2} = 37.49$$

How would I calculate the time shift and phase angle of $x_3(t)$?

Thanks in advance.


texts: Should I use CV_HAAR_SCALE_IMAGE while using LBP CascadeClassifier?

texts:
I trained a cascade-LBP to detect Lollipops with 1000 images, now I'm trying to "adapt" the openCV HAAR-Cascade example to use my LBP trained .xml but I'm not sure about the "HAAR_SCALE_IMAGE" flag, I read that optimizes the downscale to increase the HAAR features is it right to use it with LBP?

//-- Detect Lollipops
  lollipops_cascade.detectMultiScale( frame_gray, lollipops, 1.1, 2, 0|CV_HAAR_SCALE_IMAGE, Size(30, 30) );

An example of imput image:

For more info:
http://docs.opencv.org/2.4/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html


texts: FFT strong signal in middle frequency bin

texts:
I've developed a radar that I'm using to determine the distance to remote objects. It uses a custom PCB with an onboard FPGA that performs the DSP algorithms. The data from it is then plotted on a host PC. This appears as a 2D histogram where the y-axis denotes the FFT frequency bins (due to the nature of the radar this is proportional to distance) and the x-axis denotes time. The plot (shown below) gives a very strong signal at a distance halfway between the antenna and max range, which I'm unable to explain.

The actual algorithms performed are: an FIR polyphase decimation filter (downsamples from 40MHz to 2MHz) which produces an output 1024 samples in length. Then I run it through a Kaiser window function with a beta of 6, followed by a 1024-point FFT, the result of which is transmitted to the host PC. For each value of t in the plot (the x-axis) the host PC averages over 30 1024-length sequences (averaged element-wise). Since all inputs to the FFT are real, the output is Hermitian symmetric and so I only plot the first 512 values of each output sequence. The strong signal you see above occurs at bins 257 and 258 (indexed from 0). I've tested the radar in an open space where it shouldn't generate any strong signals. I've simulated all of the FPGA logic, so while I can't be sure it's right (I've only formally verified parts of it), I'd be surprised if it wasn't.

What could be the cause of this? Is there some obvious aspect I'm missing? If any of this is unclear or some part of the information I've omitted is important for answering this question (e.g. the equation relating frequency to distance), please let me know and I'll include it.

Edit: more details on acquired signal

This is an FMCW radar. A frequency synthesizer generates sawtooth ramps from 5.3GHz to 5.9GHz over a duration of 1ms. This signal is simultaneously transmitted and mixed back in with the reflected signal. We then measure the difference frequency to back out the distance.

The FPGA modules are timed such that data is only acquired during the synthesizer's ramp period. First, I enable the ramp and power amplifier and (once enabled) begin acquiring data. The data is processed by the FIR filter and then passed through the kaiser window. Once the last sample passes through the FIR filter the ramp and power amplifier are disabled. The processed data (which were stored in a FIFO) are now run through the FFT and then the resulting output is dispatched in packets to the host PC via USB. I use a header sequence, tail sequence, and duplicated packets to try to avoid data corruption/loss. Once the FFT is finished, the process starts again (ramp and power amp enable, etc.).

The FIR filter should take just longer than 0.5ms to acquire all samples, so it should fall within the frequency ramp period.


texts: Reducing a signal to Acos(omega * t+ phi) with negative omega*t

texts:
I am trying to solve a phasor addition problem, reducing from its original form to
$$X(t) = A \cos(\omega_0 * t + \phi)$$

The original equation is :
 $$X(t) = 2 \sin(\omega_0 * t + 45) + \cos(\omega_0*t)$$

I converted the sine term into cosine. The result that i got was
$$X(t) = 2 \cos(\pi/2 - \omega_0*t - \pi/4) + \cos(\omega_0*t) $$
Which reduces to :
$$X(t) = 2 \cos(\pi/4 - \omega_0*t) + \cos(\omega_0*t) $$
Up to this point I verify in MATLAB that the original function and the function in cosines are equal by graphing them and then calling the $norm()$ function and verify that it is something very close to 0.
http://www.mathworks.com/help/matlab/ref/norm.html

I then perform the phasor addition process.
$X(t) = Re(2e^{j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$

$$X_1= 2e^{j(\pi/4)} = 2(cos(\pi/4) + jsin(\pi/4)) = 2(\frac{\sqrt2}2+j\frac{\sqrt2}2) = \sqrt2 + j\sqrt2$$
$$X_2 = e^{j0} = cos0 + jsin0 = 1$$
$$X_3 = X_1 + X_2 = \sqrt2 + j\sqrt2 + 1 = 2.41421 + j1.41421$$
 To turn it into back into the stated form of $X(t) = A \cos(\omega_0 * t + \phi)$ i find the $R$
$$R = \sqrt{2.41421^2 + 1.41421^2} = 2.79793$$
$$\theta = arctan(1.41421/2.41421) = .529902$$

The final result being $X(t) = 2.79793\cos(\omega_0 t + .529902) $

I then go to MATLAB to plot the original signal and the reduced signal that I found but when I call the $norm(originalSignal - reducedSignal)$ I get a number much greater than 0.

I am trying to figure out why this is returning such a high number. Some of my guesses are the negative $-\omega_0 t$ in the first term of the function when I convert sine to cosine $X(t) = 2 \cos(\pi/2 - \omega_0 t - \pi/4) + \cos(\omega_0 t) $ I don't know if the process changes for a situation like that. Problems I have worked through and have gotten the correct answer always seem to have the $\omega_0 t$ as positive part inside the cosine term.

Below is the Matlab code.

t = [0 : .01 : 1];
originalWave = 2*sin(50 *t + pi/4) + cos(50*t);
reducedWave = 2.79793*cos(50*t+.529902);

subplot(2,1,1)
plot(t, originalWave)
subplot(2,1,2)
plot(t, reducedWave)

norm(originalWave - reducedWave)

Any ideas on what could be wrong with my process or answer?

I appreciate any advice or help. Thanks.


texts: Positive and negative frequencies in DFT due to frequency folding, or due to negatively indexed frequencies?

texts:
When I look for the cause of the mirroring of frequencies in DFT output, I get two types of explanations:

The first one which says the frequencies are mirrored because of the complex exponential which has a positive and negative counterpart. The poritive(real) parts add up and the negative (imaginary) parts cancel out, and that the amplitudes are spilt equally between the positive and negative frequencies.
The second one which says this happens as a consequence of sampling above the Nyquist frequency, i.e., above $\frac{N}2$. But the mirroring would still be present if we take the range $-\frac{N}2$ to $\frac{N}2$.

My question is what is the actual cause for this mirroring, and what actually changes with taking the range as $-\frac{N}2$ to $\frac{N}2$ or $0$ to $N-1$?


texts: What is the theoretical bit error rate for QPSK Matlab simulation when constellation points are transmitted (no symbol duration)?

texts:
I'm simulating a QPSK system in Matlab and I'm trying to calculate the theoretical bit error rate to compare it with the simulated ber. The formula that I'm aware of involves the term Eb/N0, but how do I calculate it, knowing just the SNR and the energy per bit? I seem to need to know the bandwidth of the signal but since I'm just sending constellation points with no carrier (so just complex numbers), there is no symbol duration, so I can't calculate N0 from N0 = Pn/B.

Any help would be greatly appreciated.


texts:
Below is my MATLAB code for OQPSK communication:

    data_in = [ 2 1 3 1 2 3 1 1 0 2 0 2 3 2 3 1 1 2 0 3 1 3  3 0 2 1 ];
    data_modualted=oqpskmod(data_in);
    filtmsg = filter(H,1,data_modualted); % Introduce channel distortion.
    eq1 = lineareq(nweigth,normlms(StepSize)); % Create an equalizer object. eq1.SigConst=oqpskmod(0:3);
    [symbolest,filtmsg_eq] = equalize(eq1,filtmsg,data_modualted(1:trainlen)); % Equalize.
    %%Scatter plot
    h = scatterplot(filtmsg,4,0,'g.');
    hold on ;
    scatterplot(symbolest,4,0,'kx',h);
    scatterplot(eq1.SigConst,1,0,'r*',h);
    legend('Filtered signal','Equalized signal','Ideal signalconstellation')

I am unable to understand which values to use to generate the phase
constellation diagram for my code?  How can i generate the correct
phase constellation diagram?

The first figure is constellation diagram of The filtered signal and the second is the constellation diagram of The Equalized signal


texts: Producing Colored Noise from a given PSD Data

texts:
I need to estimate (extract/produce) the colored noise of a given PSD Data. I have the following procedure to get the desired results.

First of all I calculate the filter coefficiencts using firls command in MATLAB by implementing Helbert Transform firls(filter_order,norm_freq,signal,'hilbert'), the problem I face is to choose the right order for the filter, I am at the moment working with hit and trial method and rougly a value of around 400 is selected which fits quite well with the PSD curve.
Secondly the white noise is passed to this filter and further amplified with a power amp. PSD of the extracted noise is again calculted and matched with the input PSD.

I want to know if this way is alright or not, moreover for passing white noise through the filter I am using Simulink, is there a way not to use time varibales and directly calculate the colored noise. In other word if I have to get rid of Simulink block. I tried to work with lsim and I can get to colored noise but it's amplitude is lesser and I am not sure how can I amplify the output signal so that it matches the input PSD.

just a short flow:

INPUT PSD --> Filter Coefficient Calculated --> White Noise Passed Through Filter --> Colored Noise at a lower amplitude extracted --> Amplified --> PSD of Colored Noise Calculated --> INPUT PSD matched with OUTPUT PSD


texts:
In Matlab we use below function to read an image :

var = imread('IMAGE PATH');

But, how to do same thing in Scilab ?


texts: Uniqueness of Fourier Series Representation and the Fourier Transform of Periodic Signals

texts:
If we are given a signal of the form $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{j k \omega_0 t},$$ can we call it a Fourier Series representation of $x(t)$ right away?

Suppose we are given the signal $x(t) = e^{j2\Omega t} + e^{j4\Omega t}$, this signal can be expressed as $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{jk\Omega t},$$ however, in this case, $2\Omega$ is the fundamental frequency of the signal, not $\Omega$ as the equation above might suggest.

I ask this because in the second edition of Signals and Systems by Alan Oppenheim, he derives the Fourier transform of a periodic signal by considering the impulse train $$X(j\omega) = \sum_{k = -\infty}^{+\infty} 2\pi a_k\delta (\omega - k \omega_0)$$
and applying the inverse Fourier transform to obtain $$x(t) = \sum_{k = -\infty}^{+\infty} a_k e^{j k \omega_0 t},$$ which he says is the Fourier series representation of the signal without further discussion. But how do we know that in this case $\omega_0$ is the fundamental frequency of the signal? Couldn't it be the case that it is a signal of the form given above, namely, $x(t) = e^{j2\omega_0 t} + e^{j4\omega_0 t}$, in which case the fundamental frequenct is $2\omega_0$? Thank you in advance.


texts: What frequencies are present in the Fourier transform of the Dirac impulse?

texts: When I do the Fourier transform of the Dirac impulse I get a pure sinusoid (or complex exponential, however you wanna call it) but I read in several places that all frequencies are present in the dirac impulse and all of them with the same amplitude. How is this possible? Am I wrong when I perform the transform?


texts: Signal in frequency domain with OpenCV dft

texts:
I am experimenting with cv::dft: a 1HZ sinus signal is generated, and displayed in the frequency domain. But for some reason it hasn't got the maximum component at 1Hz. My code is the following:
       const int FRAME_RATE = 20;             //!< sampling rate in [Hz]
       const int WINDOW_SIZE = 256;

       double len = double(WINDOW_SIZE)/double(FRAME_RATE);   // signal length in seconds
       double Fb = 1./len;                                     // frequency bin in Hz

       // Constructing frequency vector
       std::vector<double> f;
       double freq_step = 0;
       for (int i = 0; i < WINDOW_SIZE; ++i)
       {
          f.push_back(freq_step);
          freq_step += Fb;
       }

   // Create time vector
   std::vector<double> t;
   double time_step = 0;
   for(int i = 0; i<WINDOW_SIZE; ++i)
   {
       t.push_back(time_step);
       time_step += 1./double(FRAME_RATE);
   }

   // Creating sin signal with 1Hz period
   std::vector<double> y;
   for(auto val : t)
   {
       y.push_back(sin(1*FRAME_RATE*val));
   }

   // Compute DFT
   cv::Mat fd;
   cv::dft(y, fd, cv::DFT_REAL_OUTPUT);
   fd = cv::abs(fd);

If I plot the signal in time and frequency domain: plot(t, y); plot(f, fd) the result is the following:

The time signal is good, but the frequency signal has maximum around 6HZ instead of 1HZ.
Where did I take the mistake?


texts: How do I filter audio with modified spectrogram?

texts:
Say I computed Spectrogram of some audio input using scipy.signal.spectrogram.

Is there an easy way in python to use it as a filter for my audio input to produce an other audio which would basically be just silence or very near to that?

In practice I would modify that spectrogram before applying it back to the audio to filter out only some audio information.


texts: Is it necessary to apply some window method, to obtain the FFT? Java

texts:
I read that the window is needed to apply filter (Low Pass, High Pass, Band Pass or band Stop),  but not for the FFT.

1.
What advantage does apply a window or disadvantage does not apply?

2.
What is the criteria to choose the size (My LengthW:Length of Han Window)?

3. Where I can find literature that specifically?

My code:

Suppose an array of samples stored in

integer SIZE = 16384; //SomeValue
short[] sSamples = new short[SIZE];
sSamples = obtainSamples();  //Some method to fill my array

Now I need to obtain my FFT.
I arbitrarily chose 4096 for my FFT, then are 8 times (or 8 fft's)

FFTLength = 4096;
SIZE = 16384;
double Overlap = 0.5;
SIZE/(FFTLength*Overlap) = 8;

DIRECT METHOD:

complex cSamples = new complex[SIZE];
cSamples = short2Complex(sSamples);

Now I need to obtain part of array of sample 8 arrays of length:SIZE/4

//getPartArray(from, until) maybe like System.ArrayCopy
complex[] cSamples00 = getPartArray(0,4095);
complex[] cSamples01 = getPartArray(4096,8191);
//similarly for all, from 00 until 07
complex[] cSamples01 = getPartArray(12287,16383);

Finally, Obatin the FFT

//fft(complex[], length)
complex[] cFFT00 = fft(cSamples00,4096);
//similarly for all, from 00 until 07
complex[] cFFT07 = fft(cSamples00,4096);

WINDOWED METHOD:

//int LengthW = ??  // (besides positive and odd)
double[] Wn = new double[LengthW];
//Suppose Han Window,can be Hamming, Kaiser, Blackman, etc
Wn = obtainHanWindow(LengthW);
double[] sWindowedSamples = applyingWindow(Wn,sSamples);

complex cSamples = new complex[SIZE];
cSamples = short2Complex(sWindowedSamples);

Now I need to obtain part of array of sample 8 arrays of length:SIZE/4

//getPartArray(from, until) maybe like System.ArrayCopy
complex[] cSamples00 = getPartArray(0,4095);
complex[] cSamples01 = getPartArray(4096,8191);
//similarly for all, from 00 until 07
complex[] cSamples01 = getPartArray(12287,16383);

Finally, Obatin the FFT

//fft(complex[], length)
complex[] cFFT00 = fft(cSamples00,4096);
//similarly for all, from 00 until 07
complex[] cFFT07 = fft(cSamples00,4096);

The only difference is calculate a window and apply it to the input.

Definition's class and methods

The complex class is defined like:

class complex {
  public double real = 0.0;
  public double imag = 0.0;
}

short2Complex:

short2Complex(short[] intSSample) {
  complex intCSample = new complex[intSSample.length];
  for (int i = 0; i < intSSample.length; i++) {
    complex tcomplex = new complex();
    tcomplex.real = intSSample[i];
    intCSample[i] = tcomplex;
  }
  return intCSample;
}

obtainHanWindow:

double[] obtainHanWindow(int Size) {
  //Size must be Odd and Positive integer
  double[] dHan = new double[Size];
  int HalfSize = (int)Math.floor((double)Size/2);
  for (int i = 0; i < HalfSize; i++) {
    int j = 2*HalfSize - i;
    double Pi2DivSize = 2.0*Math.PI/Size;
    dHan[i] = 0.5+0.5*Math.cos(Pi2DivSize*(double)(i - HalfSize));
    dHan[j] = 0.5+0.5*Math.cos(Pi2DivSize*(double)(HalfSize - i));
  }
  dHan[HalfSize] = 1.0;
  return dHan;
}

applyingWindow

applyingWindow(double[] dWn, double[] dSamples) {
  double dOutput = new double[dSamples.length]
  for (int i=0; i < dWn.length; i++) {
    dOutput[i] = 0;
    for (int j=0; j < dSamples.length; j++) {
      if(i-j>=0) {
        dOutput[i] += dSamples[j]*dWn[i-j];
      }
    }
  }
  return dOutput;
}


texts: Do we need Channel estimation prior to (time domain DFE or linear) adaptative equalization in single tone transmission?

texts: If we are to transmit a single tone signal over a channel impaired by inter symbol interference and awgn, and if the equalization is adaptative (LMS/RLS track) and is done in time domain by LE (linear equalizer) or DFE, do we need to compute the channel impulse response first or the equalizer simply track channel changes when it computes it coefficients ?


texts:
I'm trying to write my own high quality audio sample rate converter. I barely know anything about signal processing though so I need help.
From what I understand I need to sum together normalized sinc functions that touch every sample to find the value at any arbitrary point. I'm guessing this is hard because that would mean more than a billion calculations for a few second long audio file.
So do I only use the sinc functions for the samples that are closest to my x value? Is this what is meant by "windowed sinc"? How many samples should I go in each direction away from my x value? Additionally, my DAW has something called "32 point sinc" resampling. What is the "32 point" supposed to mean in this case?


texts: How to predict next image in a sequence, or predict the deblurring of an image?

texts:
I have 3 versions of an image at lower resolutions - 480p, 720p, and 1080p. Is it possible to use those 3 images to predict the next highest resolution image (the 4k image)?

My first guess was that I could calculate the difference (or momentum or acceleration?) of each pixel between images and use that to predict the new 4th image. I assume this would require me to give all pictures the same number of pixels. I also assume that the differences in each pixel between images would need to be smooth. Would this work?

My other guess was that maybe there is some kind of statistics or machine learning technique that I could use, but I don't know anything about that stuff.

Thanks!

Ps - this man is not my real data. My real images are an intensity map (temperatures over an area) at different resolutions. I suspect that the introduction of colors might complicate things, so should we pretend that the images below are black-and-white images?

Here are the pictures:


texts:
A question about using fir1/fir2/firpm to design FIR filters in MATLAB.  The first parameter is the 'filter order'. For my filter design I pass in 250 as this argument.  However, when looking at the impulse response of the filter, it is a type 1 FIR filter.  So my question is, is the filter order not the same as the number of taps? Is the number of taps one less, so 249 in this case, which gives the type 1 symmetry?

Cheers


texts: How do Anti-Aliasing Filters filter the signal?

texts:
I have programmed a anti-aliasing filter that uses the fourier-transformation and makes a perfect cutoff on my desired cutoff frequency for downsampling without aliasing.
I have read about other types of filters for anti-aliasing. Here i have a nice picture of different filters from this site:

The dotted line seems for me the perfect or ideal filter, but in practical it is not possible to filter a signal perfectly right? And therefore we can use different filters with different orders.
Does this graph mean, the different filters do not filter completely on the perfect cutoff frequency and always contains a little aliasing?

The programmed FFT-Filter is used as followed in python:

    def lowpass(x, binmax):
        N = len(x)
        return np.fft.irfft(np.fft.rfft(x, axis=0)[:binmax], N, axis=0)

x48hr is a Signal sampled with 48kHz and do a cutoff at 4kHz
lr_signal = lowpass(x48hr, binmax=len(x48hr)*4//48) #hr

I use this filter for generating low-resolution samples to train an convolutional neuronal network to learn the mapping between lr and hr to make a hr signal reconstruction.

In general I am actually trying to understand how much aliasing my measurements have and how I can use my measurements to study different types of filters and the aliasing in them. I want to increase the resolution of my signal, therefore I need aliasing in my low-resolution signals. I have asked another question regarding my problem here if you want to know more about it.


texts: Extract Impulse Response of a System From its Input/Output Waveform

texts:
I am a new engineer for signal processing field.

I want to extract the system impulse response from input and output of the system, i.e. so-called system identification, more specfically, for impulse response. E.g. I could input waveform to system as following:

And I also got the output waveform from the system as following:

Please note both waveform are four periods, but there is unknown delay existing in between.

How could I extract the impulse response of system from those two?

Update in 2020/04/27

I followed the comments from Marcus, and used following code, but the result is obviously not correct, not even close.

x = load("C:\Users\Desktop\bit_wf.txt")
x = x(:,2)
plot(x)

% convolution between bit_waveform and impluse response
% wf = conv(bit_wf, imp_rsp);
y = load("C:\Users\Desktop\wf.txt")
y = sbr(:,2)
plot(y)

ryx = xcorr(y,x)
plot(ryx)

rxx = xcorr(x,x)
plot(rxx)

imp = deconv(ryx, rxx)

plot(imp)

If I deconv(ryx, rxx), I only get a constant number. I understood the theory is there, but my problem is how to implement it correctly. Thanks, hope someone could help me there, really appreciate.


texts: Difference between video mosaicing and video stitching

texts: What is the difference between video mosaicing and video stitching ,is this both are similar or is there any methodology difference?


texts: Hardware implementation of multi rate FIR filter

texts:
I have designed a multi rate fir filter in MATLAB and I want to implement it in verilog.This is the block diagram.
Total number of taps are 391. Is there any way to automatically generate verilog code? MATLAB hdl coder is giving a lot of errors.
Also, My input sampling rate is 44.1 kHz and output sampling rate is 48 kHz. Signal passband and stopband are 1kHz and 1.1kHz.
So what should be the clock frequency of G(z),I(z), K(z) and P(z) (I,K,P(z) are the anti imaging filter). please help!


texts: What is an example of using aliasing to your advantage when recovering an input signal?

texts:
Suppose you have an arbitrary analog input signal $x_a(t)$ guaranteed to have frequencies within a bandwidth $[f_1,f_2]$ Hz.

Suppose your sampling frequency $F_s$Hz, and sample $x_a(t)$ to produce $x(n)=x_a(\frac{n}{F_s})$.

Then by the sampling theorem, you can successfully recover (or reconstruct) the signal so long as $f_2 < \frac{F_s}{2}$. Otherwise you will experience aliasing and cannot recover the original signal $x_a(t)$.

It appears (from this article) that there are "tricks" that can be employed that take advantage of aliasing scenario and can still recover the signal $x_a(t)$ completely.

The specific part of the article of concern is pasted below:

  Using Nyquist aliasing as benefit

  The trick is to use the aliasing (or frequency folding) to your
  advantage. By undersampling the data converter, higher-frequency
  content will be aliased into all of the lower Nyquist zones (see
  Figure 2). You will need to make absolutely sure that nothing ends up
  in the lower bands – any noise or frequency components in the lower
  zones will also be aliased into the first Nyquist. The good news is
  that the data rate from the data converter is only a fraction of the
  required RF input sample rate if this were a first Nyquist system.
  Under sampling greatly reduces the data rate of the samples supplied
  to the digital signal processor (DSP) or FPGA.

Could somebody explain this to me? Why is the Texas Instruments link not violating Nyquist's Theorem?


texts: Solving the PDE of a signal in an optical fiber

texts:
Disclaimer Please let me know if this is not the right place to ask this question.
Disclaimer Part 2 I have no background in signals. My focus is in mathematics.

In my programming class, we are studying the pulse propagation in an optical fiber. The fiber can be modeled as a line marked by a one-dimensional coordinate $z$. The pulse at the fiber input ($z=0$) is given as a time-dependent field envelope $A(0,t)$. I need to derive the magnitude of the field envelope at any location and at any time $A(z,t)$. We are told to use numerical methods (i.e. Fourier transform) to solve the partial differential equation

$$\dfrac{\partial A}{\partial z}
= \dfrac{\alpha}{2} A
- \beta_1\dfrac{\partial A}{\partial t}
-j\dfrac{\beta_2}{2}\dfrac{\partial^2 A}{\partial t^2}$$

Here, the first term on the right-hand side models the signal attenuation due to scattering, the second term models the group velocity, and the third term models chromatic dispersion (i.e. different colors travel at different speeds in the fiber, resulting in pulse broadening). My teacher has said to use a Gaussian pulse and to discretize the pulse over the range $\left[-\frac{t}{2}, \frac{t}{2}\right]$ such that $t >> \sigma$ $\left(t=10\sigma\right)$.

I have been able to show mathematically that the solution should be of the form
$$A(z,t) = \text{ifft}\left[\tilde{A}(0,\omega)\exp\left(\left(\frac{\alpha}{2} - j\omega\beta_1 + \frac{j\omega^2}{2}\beta_2\right)z\right)\right] $$
(where $\tilde{A}(0,\omega)$ is the Fourier transform of $A(0,t)$) but am having trouble actually implementing these results in MATLAB.

What I have so far:

sigma = 0.1; % standard deviation of the pulse
N = 256;  % number of points sampled
t = linspace(-10*sigma, 10*sigma, N); % time vector
A0t = (0.25/sqrt(2*pi*sigma^2))*exp(-(t.^2)/(2*sigma^2)); % descritized gaussian pulse

Fs = 1000; %sampling frequency

(I don't know what this number should be)

A0w = fft(A0t);
A0w = abs(A0w(1:round(N/2)));

w = (0:(N-1))*(Fs/N);  %frequency vector
w = w(w < Fs/2);

(I don't know why the last three lines are necessary. I saw it here.)

alpha = input(...);
beta1 = input(...);
beta2 = input(...);

lambda = ((alpha/2) - (1j*w*beta1) + (1j*(w.^2)*beta2/2)); %eigenvalue of ODE in frequency domain

end = input(...);

z = linspace(0,end,length(A0w))'; % spatial vector

Azw = A0w.*exp(lambda.*z); % solution in the frequency domain

Azt = abs(ifft2(Azw)); % solution in the time domain

Using alpha = 10, beta1 = 50, beta2 = 18, and end = 10 results in the surface plot below, which I do not believe is correct. ($t$ and $z$ are plotted against their indices rather than the actual values because $t$ is not the same size as $A(z,t)$.)

If somebody could please explain to me why I need to do the things that I am doing in my code (the italic comments) as well as what I should be doing otherwise to get a correct output, or if what I am getting is correct, it would be much appreciated.


texts: Conjugate reciprocal pairs of zeros and poles in FIR design

texts:
Assuming the impulse response $h[n]$ of an FIR filter is real for all $n$,

Why are zeros and poles in FIR design found in reciprocal and conjugate pairs?
Is the assumption necessary for this phenomenon to take place?


texts: Calculate signal to noise ratio from the Hillbert envelope of a stacked signal

texts:
I am working with seismic waveforms and more specifically with a stack of years of daily seismograms, that I would call the reference waveform. I am trying to calculate the SNR for each daily waveforms by using the envelope of the reference, but I have doubts on how I should calculate the noise.
Concretely, I am trying to reproduce what was done in this paper (page 8).

I believe It should be simple, but I have poor basics in signal processing and struggle with notations.


texts: How to design very narrow band, very sharp low pass filters - Only DC needed?

texts:
I am using a DSP processor to sample signals at 108kHz and I only want to get the DC part out of it. Such a high sampling rate is used because I want to use oversampling to reduce quantization noises.

Two filters are to be designed:

One has 10Hz stop-band frequency; new data are read every 0.1 second;
One has 0.1Hz stop-band frequency; new data are read every 10 second;

What I did

I used several stages of decimation and downsample the signals to 300Hz. Then I use FIRs to filter the signal:

For the 10Hz bandwidth: 128 stage filters are used
For the 0.1Hz bandwidth: 768 stage filters are used.

My question

Performance is the most important to me: The sharper low-pass filters are, the better.

Could I ask for suggestions in implementing such a filter?
Is it possible to design a 0.01Hz low pass filter? If so, the new data can be read every 100 seconds?
I heard that using IIRs could be a solution, is it?

Thank you very much!

Is averaging the best DC filter?

I tried to compare with matlab:

Setting 1:
%Use Matlab fir1 function

    size = 1080000;
    h = fir1(size, 0.00000000001);
    fvtool(h)

Setting 2:
%Build an averaging FIR

    size = 1080000;
    h = ones(1,size);
    h = h *1/size;
    fvtool(h)

The theory should be both two settings have similar performance. However, it seems that fir1 has narrower stopband than the simple average filter. It seems that the average DC filter will include more noises.


texts: Determine Frequency Content of Variable Rate Signal (Non Uniform Sampling of Signal)

texts:
What is the proper way to calculate the frequency content or spectrum of a signal with a variable (but measured) sample rate?  I.e. get similar results as you would from an FFT, but for data not evenly spaced in time.

I know I can resample the data by interpolating on an evenly-spaced time vector then apply a FFT, but this seems less than ideal.


texts: FM modulation audio 50 Hz – 16 kHz MATLAB

texts:
I'm trying to modulate in FM an audio signal $x(t)$ coming from a microphone which works from 50 Hz to 16 kHz in MATLAB.

But the problem is that I want a FM Bandwidth of 25 kHz at the output of modulator.

I have some problems to do this:

The FM bandwidth for any signal using Carson's rule is:

$$BT = \left(\frac{f_d \max\vert x(t)\vert}{W}+ 1\right)W$$

where $W$ is the maximum bandwidth (15.950 kHz in this case) of the modulating signal and $\max\vert x(t) \vert $ is the max value of the modulating signal in $Vp$.

The thing is I have to adjust the deviation sensibility $f_d$ (slope) to get the desired FM bandwidth, but how to know $\max\vert x(t)\vert$ in the microphone datasheet? or even using an amplifier between microphone and modulator?


texts: Why the RMS of a PSD curve is the root of the area below

texts:
I will try to explain what is my level of understanding of this problem, please correct me if I'm wrong:

RMS is the Root Mean Square, it represent the mean value of the input signal.
PSD is the measurement of the responses that shows me at which frequencies most of the energy is concentrated.
The area below a curve is the integration of that function.

My situation is that several random vibration tests are performed. These tests are called random tests because of the input signal. In contrast to a sine test where the structure is excited with a sinusoidal input, only one frequency is excited at a time, here 'all' frequencies are excited at the same time.

In this case PSD is measured in ${{g^2}/{Hz}}$ and RMS in ${g_{RMS}}$. Armed with that it easy to see that if you multiply PSD per the frequency range and you take the root of the result you will get something in ${g}$'s, but I don't know how exactly derive the famous relation:

\begin{equation}
{g_{RMS}=\sqrt{\int_{f_1}^{f_2}PSD(f)df}}
\end{equation}

The understanding that I have is very basic and it would be great if someone give me a clear idea of the relations among the RMS, PSD and the real signal. Thank you very much.

In the figure I have plotted a standard random signal.


texts: Matching two signals with small time drift

texts:
I have two time-domain signals sampled at 100 Hz that were measured using two different oscillators and therefore have a time drift between them. I have two synchronization points, one at the start and one at the end of the measurement (usually about 24h long).

The drift is very small, about 1 second for a 24h measurement, i.e. about 0.001% or 1/86400. I assume a linear model for the drift.

The length of the measurement and the relative scale of the drift prohibit usual resampling techniques: I'd have to up- and down-sample the longer signal by large factors.

My approach would be to drop every 86400th sample of the longer measurement, the questions I have are:

Is there a better approach?
What's the frequency-domain equivalent of dropping a sample? I imagine it would shift the spectrum by a tiny amount and introduce some high-frequency noise.

To be honest, the second question is more academic than practically relevant; most of the analysis is done in small chunks and I can just exclude the chunks with the dropped samples from the analysis.


texts: Discrete Fourier transform of complex time series

texts:
I am working through previous exam questions for a class I'm currently taking in signal theory, and one problem, although it looks easy, stumps me.

PROBLEM

Two real $N$-periodic time seires $f_n$ and $g_n$ are given, where $n = 0, 1, 2, ..., N-1$.

In addition, the complex time series $h_n$ is given as

$$h_n = f_n + i g_n$$

and the discrete Fourier transform

$$H_m = F_m + i G_m$$,

where $m = 0,1,. . ., N-1$.

Find the discrete Fourier transofrm $F_m$ and $G_m$ to the time series $f_n$ and $g_n$ expressed through $H_m$

ATTEMPTED SOLUTION

As I said, I'm really stumped here.  Other than writing:

$$F_m = H_m - i G_m$$

and

$$G_m = i (F_m - H_m)$$

I don't see what else to do.  Is this all there is to it?  It seems way too easy, so I'm pretty sure I'm not doing this correctly.

If anyone can give me any help, I would greatly appreciate it!


texts: Can LMS filter converge to sparse solution?

texts: Consider received signal $r=h*s$, where $h$ is the channel and $s$ is transmitted signal. Now from $r$ ( $s$ is known) can we predict $h$ using LMS algorithm under all conditions?  I was under impression that both $s$ and $h$ needs to be stationary to estimate $h$ properly. In case, $h$ is sparse (few non zero only) LMS cannot converge to sparse $h$. Is it true or false? Can someone explain my doubt pleas?


texts: Removing white noise by taking the mean of many samples

texts:
I'm learning about signal processing, and I am attempting to remove white noise from an output of a blackbox system.

As you can see from the above image (left is the input sin(t), right is the output) there is quite a bit of noise.
I am attempting to verify that it is white noise. To do, I have run the system 100 times with the same sin(t) input and have taking the mean of the outputs, which resulted in the following plot:

If it is white noise, I believe I should expect that the mean of many outputs will approximate the noise-free output for larger and larger sample sizes.
However, as you can see, there is still quite a bit of noise and I'm not too sure why. At first I thought the data between samples wasn't lining up when taking the mean, but I don't believe that should actually matter.


texts: How to calculate spectrogram from nonstationary, nonuniformly sampled signal

texts:
I have a high resolution 3D spatial signal (xy-elevation data collected via Terrestrial Laser Scan) with very nonuniform sampling (holes/shadows due to vegetation; point spread varies as a function of distance from scanner). The original point density varies from ~0-100 pt/cm2 and has been gridded/downsampled to 1 pt/cm2. Holes can be up to a few meters wide.

I now need to calculate a 1D spatial spectrogram in Matlab (in either x or y vs elevation) without interpolating data over the holes or downsampling to the lowest sampling rate since these would introduce high frequency artifacts or decrease resolution (I'm interested in the spectral characteristics at the cm-to-m-scale). The signal is also nonstationary, so I don't think I can just use a resampling or bootstrapping method.

The built-in Matlab functions I am familiar with (FFT, spect, pwelch,) result in NaN output if there are any NaNs in the signal, so I'm considering using Matlab's Lomb-Scargle periodogram algorith (plomb) in overlapping windows.

Is there another existing spectrogram method that uses least squares fitting, frequency-dependent windowing, or some other method to ignore NaNs/holes? Or is there any reason why using overlapping plomb windows would be a poor solution here?


texts:
In my previous post I asked for help for a Fourier transform of $$ t \text{rect} ( t- \frac{1}{2} ) $$ and I think I’ve understand the process. Now I’ve 2 another similar Fourier transform to do , I already solved both , but I don’t have the correct result. Can someone tell my if the 2 results I obtained are wrong ? Thank you

If $$ x(t) = t \text{rect}(t) $$ I obtained $$ X(f) = \frac {\text{sinc}(f) - \frac{1}{2} [e^{i \pi f } + e^{-i \pi f } ]}{i 2 \pi f } + \frac{1}{4} \delta(f) $$

And for $$ x(t) = \Big(t + \frac{1}{2}\Big) \text{rect}(t) $$ I obtained $$ X(f) = \frac{ \text{sinc}(f) - e^{- i \pi f } }{ i 2 \pi f } + \frac{1}{4} \delta (f) $$


texts:
Generalized/power means may be used to construct moving-average filters with different properties than regular one which is based on arithmetic mean. This observation seems to be trivial - even Wikipedia mentions that.
The problem is that power means are usually defined for sets of positive numbers - while in DSP we often deal with negative numbers too. Therefore, many of well-known inequalities and properties of power means may not be applied to sets with negative numbers.

This is of interest to me since I have found that in my partical area of application power means of higher order give better results than arithmetic mean. I have been able to prove some basic properties and I more or less understand where does the improvement come from. However, I would be interested in learning. However, everything I have found seems to focus on positive numbers. So my question is:
Are there any papers on using generalized means for moving average filtering or some other signal processing? Any papers on purely mathematical properties of generalized means for sets with negative values may be relevant too. Thanks!


texts:
The signal I'm working on includes some noise which is I assume Gaussian white noise with zero mean and an arbitrary variance. So, it is like;

y[t] = a[t]*x[t] + n[t]

I'm searching for the unknown variance of the noise component of the signal to use in Kalman filter. My idea is using Whitening filter over observation y[k] and finding variance of noise component n[k] afterwards.

Is it a common practice to use whitening filter for noise estimation and is there any counter-argument not to use it?


texts: What is averaging and how can it be done?

texts:
I'm studying about the analysis of (mainly) fMRI and EEG data. Multiple times it's mentioned that to reduce noise, you can use averaging, but nothing more detailed than that.

There never is literally stated that you can use this or that technique to average.
So, what exactly is averaging? What techniques can be used to achieve it?


texts: Given a system with Transfer Function and its desired output. Is it possible to find the required Input?

texts:
I'm currently working on an DC servomotor. I managed to find it transfer function using System Identification theory and Matlab.

The transfer function is given as
$$H(Z) = \frac{1 - 0.4952z^{-1} + 0.01826z^{-2}}{1.612z^{-1} - 0.02492z^{-2}}$$

Now, I have a pre-defined velocity profile. That mean at each sampling time $nT_s$, I have to reached a certain speed. In other words, I have a matrix of required output at time $nT_s$.

Now my question is, is it possible to calculate the input to the transfer function $u(t)$ so that the output approximates the desired velocity? If it is possible, could you kindly give me some broad instruction or some reference so that I can study on it? If no, is there any possible theory that can solve my problem, that is find required Input to get desired Output?

My process so far: Well, I simply inverse that transfer function and ask Matlab to find the required input for me. However, the inversed is not proper and thus it is not possible to do that. So I add some fast-orders (first time hearing this, not sure what it means) to the denominator to make it proper. The function runs, but the results is quite... not right.

This is my desired Output vRequired, stored in a $1$x$714$ matrix and plot for clarity.

And this is the result with the function lsim(reservedTf, vRequired):

where reservedTf is my transfer function reserved with fast-orders:

$$  H^{-1}(z) = \frac{z^{-1} - 0.4952 z^{-2} + 0.01826 z^{-3}}{1.612*10^{-10} + 3.224*10^{-5} z^{-1} + 1.612 z^{-2} - 0.02492 z^{-3}} $$

I very much appreciate all your help. Please tell me if there is anything wrong with my post here.


texts: How do I calculate the autocorrelation of a train pulse having only the pdf?

texts: I have a pulse train of a fixed period, but a random amplitude following a Gaussian distribution $(\mu=0,\sigma=1)$. I know how to work with a random period or delay, but how do I proceed to calculate the autocorrelation with a random amplitude?


texts: Most generally-accepted audio term, for "pre" and "post" lengths of a symmetric filter's impulse response?

texts:
Consider a FIR audio filter with symmetric impulse response.
In my case, the filter kernel width is 15 frames on each side (so including the center frame, 31 frames total).  When examining the audio impulse response of this filter, it includes (A) the initial interval where it approaches the center, then (B) the centered impulse output frame itself, then (C) the subsequent interval following the center.
I have seen numerous terms referring to intervals (A) and (C), including

(A) lookahead, ramp-in, processing delay, pre-echo, positive filter width, ring-in, pre-response;
(C) decay time, ramp-out, post-delay, negative filter width, ring-out, causal response.

... and many others.
In the audio world, is there a term that is considered most accurate and/or most clearly descriptive for these two intervals?


texts: I understand the four different type os FIR filters however i don´t understand why do linear-phase FIR filters offer opportunities for efficient realization.


texts: Can you present the convolution of sinusoidal with itself?

texts:
Ladies, Gentlemen,
Because I am homeless (in France), and get internet access only in public libraries with many restrictions in timing etc, I can not write down even a simple convolution. So I ask you post here the convolution of  some sinusoidal function with itshelf.

For example: 0.951056516295, 0.587785252292, -0.587785252292, -0.951056516295, 0.

Regards


texts: Filter long audio data using cascading biquad filter

texts:
I would like to lowpass filter a long pcm wav file. I designed a 20 order butterworth lowpass filter in Octave and transform the zero pole gain data to a series of second order sections. Then I am applying these SOS filters to each processing block on the long audio. The script is like this:

clear all; close all;
[y, fs, bit] = wavread('sweep.wav');

# design butter
[z, p, k] = butter(20, 1200 / (fs / 2), 'low');
[sos, k] = zp2sos(z, p, k);

# frame index
framesPerBuffer = 512;
currentFrame = 1;

# final condition buffer
zbuf(2) = 0;

while currentFrame + framesPerBuffer < size(y, 1)
    # read one frame
    data = y(currentFrame:currentFrame + framesPerBuffer - 1, 1);
    for i = 1:size(sos, 1)
        data = filter(sos(i, 1:3), sos(i, 4:6), data);
    end
    y(currentFrame:currentFrame + framesPerBuffer - 1, 1) = data;
    currentFrame += framesPerBuffer;
end
# apply gain to filtered signal
y = y * k;

The "sweep.wav" input is a mono pcm file with 48kHz sample rate and 16 bit depth with duration of 40 seconds. It is a step sweep signal with each sine tone duration about 1 second, frequency range from 200-12000 Hz. The overall wave shape is like this:

After the filtering process above the output signal has the wavform of this:

Since butterworth filter has the flat passband/stopband, the waveform is clearly not correct. When I listen to the filtered wave it is like a filtered wav plus some low frequency modulation.

I have also tried to add final condition z buffer for the filter function but it is still perform bad. Is there anything wrong for me filter process? How can I get the correct filtered data out of this long PCM? Thank you very much!


texts:
A DT system is described as follows:
$$y[n] = T(x[n]) = x[n] + \sum_{l=1}^\infty (0.9^l(x[n-l]+x[n+l])) $$

I have found that the Discrete-time necessary and sufficient condition for stability is.
$$ \sum_{n=-\infty}^\infty |h[n]| = ||h||_{1} < \infty $$

and in my case:

$ y(n) \le  B_{x} + \sum0.9^{L}(B_{x}+B_{x}) = B_{x} (1+2 \frac{2}{1-0.9} - 2) = 19B_{x} < \infty $

Is it enough to prove system stability?


texts: sampling a sinusoidal signal at the Nyquist rate

texts: To sample a sinusoidal signal at the Nyquist rate, we have to take two samples per cycle. Is there any condition on the position of the two samples? is it necessary to take them at the peaks?


texts: Current fluctuations with reduced capacity

texts:
I simulate a neuron by compartments, where each compartment is represented by an electrical circuit. Myelinated nerves have reduced capacitance - so for example, if $n$ denotes the membrane layers, then the capacitance is C/n measured in [uF/cm^2]. Shouldn't random current fluctuations become smaller when capacitance is reduced? However, they do become larger - i.e. the random SCh fluctuations for a capacitance of
$C = 1 [uF/cm^2]$ are smaller than for $C = 1/5 [uF/cm^2]$ for example.
But does that make sense?

I model a neuron with a multi compartment model of the Hodgkin Huxley type - this means that the membrane voltage of each compartment can be determined with the HH equations. this figure shows the soma and with two adjacent compartments. If I reduce the capacitance of the soma, then (measured at the end of the axon) the fluctuations in the membrane voltage increase instead of decreasing.


texts: Modelling transmitter and receiver impairments in communication system

texts:
A paper I read studied the effect of transmitter impairments and receiver impairments by modeling both as additive white Gaussian noise. The noise for transmitter impairmenets are added before the antenna before transmission. The noise modelling receiver impairments is added as an additional white Gaussian noise to the thermal white Gaussian noise. A conclusion was made as follows

  The main reason for this is that the noise added at the transmitter (modeling Tx EVM) experiences the same frequency- and time-selective fading as the desired signal. In contrast, the noise added at the receiver (modeling Rx EVM) does not experience frequency-selective fading (only time-selective fading).

Does the quoted conclusion above make sense?


texts: How to change a continuous signal to form a discrete signal?

texts:
I have a continuous signal $$ x(t)= 2\cos(200\pi t)+8\cos(250\pi t)$$ and it is sampled at 400Hz.

How is this changed into the discrete time domain in order to determine the two discrete frequency components of $x[n]$ ?


texts:
For STFT, we impose window of certain size onto the original signal, then we perform fft on each window. The uncertanty about frequency and time is determined by the width of the window, however, I can't understand what is the point of having overlap windows...

If we have a signal, for instance, why can't we just divide the signal into 6 trunks (non-overlapping window), and then we perform fft on each of those trunks?

Maybe let me make it more clearly in my application. I am going to mostly dealing with 60Hz power line wave, and occationally, we want to monitor the 180Hz transient effect at the power line. Since the signal will be mostly periodic, should I use window then?


texts:
When I read papers about channels, I usually see something Quasi-static Channel,

Does "Quasi-static Channel" mean "time-variant channel" or "time-invariant channel"?


texts: Doppler Shift and the Pulse Repetition Interval

texts:
Suppose that I have a signal $x(t)$ consisting of $N$ pulses at a given frequency $\omega_{c}$. I space them $T$ seconds apart, ie the Pulse Repetition Interval (PRI) for the waveform is $T$.

Now this signal bounces off a target and experiences a large Doppler shift, such that the compression of the pulses are evident. My question is -- does the PRI change as well?

In other words, are each of the waveforms still $T$ seconds apart, or does the PRI scale with the Doppler shift?

I guess the question really is, should the model of the second pulse be:

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T - \tau}{v}) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T}{v} - \tau) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t}{v} - T - \tau) $$

Where $\tau$ is the target delay and $\frac{1}{v}$ is the Doppler shift.

Note: A similar question is asked in the DSP Stackexchange question
Doppler Shift of Phase Coded Pulse Compression Waveform


texts:
I'm reading 'The Scientist and Engineer's Guide to Digital Signal Processing'. It says for a noisy signal, the important parameter is not the deviation from the mean, but the power represented by the deviation from the mean.

So, if I were to consider Amplitude it would just be voltage deviations from the mean, but for power it becomes Voltage^2 deviations. Then an average is taken and a square root to give the standard deviation [It would just be an average deviation in the case of the Amplitude].

I'm trying to understand why considering just the Amplitude is not good enough?
Why is it that when I'm dealing with noisy signals, I need to consider power calculations to try and reduce the noise?


texts: How to implement bandpass filter on complex valued signal?

texts:
I am using Scipy to implement bandpass filter but it assumes that positive normalized frequency is passed but I & Q samples range from [Fc-Fs/2,Fc+Fs/2] where Fs is sampling frequency & Fc is centre frequency but desired bandpass filter should filter from [f1,f2] such that they can negative too. How to pass this band information ?

My approach is

Generated a window using Scipy firwin
Convoluted it with signal using Scipy fftconvolve


texts: (FFT) Is it good practice to ignore each bin's phase when characterizing a signal

texts:
I have been recently working with the FFT and realized that each frequency within a bin corresponds to a phase. I understand why a phase that corresponds to each frequency is needed because it allows one to construct the original signal from those sinusoidal waves.

Currently I am trying see if a specific signal corresponds to a specific sound.

My Question

When characterizing a signal with FFT, is good practice to just use the frequency of the bins and ignore the corresponding phases?


texts: Where 2-D DFT is necessary in image processing applications instead of 2-D DCT?

texts:
As we know the 2-D DCT is a real-valued kernel and less computational complexity and it used in several image processing applications like image compression etc. 2-D DFT is a complex kernel and high computational complexity, So in which image processing application of 2-D DFT is more advantageous than 2-D DCT?
Why DFT is instead of DCT in image processing applications?
Any suggestions would be greatly appreciated!
Thanks
Sukumar.


texts: Determining effect of powers of DFT matrix conjugate on input without multiplying

texts:
Suppose I have a DFT matrix $F$ of dimensions $N \times N$ where $N$ is a power of $2$. I convert it into a unitary matrix $\displaystyle U = \frac{F}{\sqrt{N}}$ and compute $S = U^* U^* $.

$S$ seems to have a pattern. $S$ has at most one $1$ per row. Applying $S$ on my $N$-sample signal $x$ essentially provides me either a permutation of $x$(if $N$ is an even power of $2$) or a subset of $x$(if $N$ is odd power of $2$). By subset, the length of $x$ remains the same, but some of the elements get truncated to $0$

Translated to MATLAB code for $N=8$,

N=8;
F = dftmtx(N);
U = F/sqrt(N);
S=ctranspose(U)*ctranspose(U);

$S$ above looks like the following:

My problem is the following: Given $N$, can I determine how to pick elements of $x$ without having to resort to multiplying by $S$ ? (For the curious, let's say the matrix multiplication is expensive for my situation).


texts: How does resizing an image affect the intrinsic camera matrix?

texts:
I have a camera matrix (I know both intrinsic and extrinsic parameters) known for image of size HxW. (I use this matrix for some calculations I need).

I want to use a smaller image, say: $\frac{H}{2}\times \frac{W}{2}$ (half the original).
What changes do I need to make to the matrix, in order to keep the same relation ?

I have, $K$ as the intrinsic parameters, ($R$,$T$ rotation and translation)

$$\text{cam} = K \cdot [R T]$$

$$K = \left( \begin{array}&a_x &0 &u_0\\0 &a_y &v_0 \\ 0 &0 &1\end{array} \right)$$

$K$ is 3*3, I thought on multiplying $a_x$, $a_y$, $u_0$, and $v_0$ by 0.5 (the factor the image was resized) , but I'm not sure.


texts: Pyramid Match Kernel: How to fit grid around data points?

texts:
I am planning to use Pyramid Match Kernel method for object recognition from depth images: I am going to extract a set of feature vectors $x \in \mathbb{R}^d$ for each object instance and then I want to use the Pyramid Match Kernel as described in http://jmlr.csail.mit.edu/papers/volume8/grauman07a/grauman07a.pdf in Grauman's paper, in order to use the set of feature vectors for each object with a SVM classifier.

I have a problem to understand how exactly we are supposed to build bins around the data points to begin with. The paper gives the following definition:

What does "inter-vector" distance mean here? I think of the following: We consider all $x_i$ values for each $X$ in $S$. Then we find two vectors with smallest distance $d$ ,(unique ones) and then scale all vectors with $1/d$. Then $D$ is the maximum absolute value of a vector element in $S$.

Is this what the authors meant here? What about if we are going to test a new vector $y$, how should we scale it?


texts: I have two sets of inputs, A (10 values) and B (20 values), and for each point (A,B) I have a measurement to make a (10x20) table of measurements C. Is there a way to use a Kalman Filter to improve these measurements if I have no physical rule to go from the inputs A and B to the value of C?


texts: How to find the best Q format representation?

texts: I'm currently working in a project in which I use the coefficients of a IIR Elliptical digital filter (9th order) and the voice signal from a recording of my cellphone but I'm having issues finding the best Q format to work its implemenation, any suggestions? Till now I've been analysing which are the maximum and minimum values of those two sources in orden to determine which would be the biggest or smallest result from their product, is this a good approach to solve this problem? How can anyone determine the best Q format to work in a system?


texts: How can I redefine median filter function instead of using the available function in Matlab? I want to smoothen my sound chart


texts: How to remove text region in a document image? or extract image from a document image

texts:
Given a document image (i.e. newspaper), how to extract photos in it or remove text region?

I think traditional OCR methods may not be suitable here, as I don't need to recognize the text, and OCR is not accurate and slow. I believe text region (i.e. text blocks) and image region should be distinguishable by some threshold based methods in image processing. Any suggestions or example codes in OpenCV will be appreciated. Thanks！

BTW, what if the background color is not white, or the background color of certain blocks are not white?

Example image:


texts: How to generate such sounds as in example?

texts:
I found these files on Internet:

1.wav
2.wav
3.wav
4.wav

How to generate such sounds?

Geno Chen wrote "Acoustic Grand Piano (Instrument #0)"

What to do next?


texts: Texture-like measures for quantifying density of data in binary images

texts:
Consider the following black-and-white image. It depicts a freehand sketch.

I wish to characterize the "density" of sketch strokes. For e.g. the hair strokes are densely grouped together. So are strokes near the wrists and the necklace stone. Other strokes are somewhat scattered and "far", e.g. the nearly vertical strokes depicting the dress.

Is there a good measure (texture-like ?) which can be used to quantify the above notion ?


texts: Why is doing fft on a non-stationary signal a problem?

texts:
Why is it a problem to do frequency analysis on a non-stationary signal?

what makes the frequency interpretation incorrect?


texts: QAM constellation slicer with unknown carrier phase

texts:
As shown in the image, I have a 16QAM constellation that is misaligned due to a phase rotation. In this case you can see that rotation amount is approximately $\theta = \pi/4$, but this won't be the case in general. For real-world data the phase could be a slowly varying function of time, $\theta(t)$, so that it is not enough to apply some fixed correction factor.

I am aware of differential mapping schemes that solve the phase ambiguity problem due to the constellation having $\pi/2$ symmetry, but it seems $\theta$ must still be known to perform the slicing.

One suggested solution was to try to map the received constellation point to the nearest QAM constellation point and feed a phase-locked-loop with the result, but it is not clear how this would perform when $\theta$ varies over time.

What techniques exist to recover the symbols? I have already tried various carrier recovery schemes based on feedback loops, with no success, and am interested in decision directed approaches that may avoid having to find the phase.


texts: Mysterious notation used in "BLIT" CCRMA paper

texts:
I'm working on a project involving synthesizing audio, which led me to a promising Stanford CCRMA paper on band-limited synthesis. However, this double-arrow notation keeps popping up throughout some pretty critical portions of the document, and I've been unable to find any mention of it in the paper itself (you can imagine running a Google search for this might be tricky). It seems as though it has something to do with sampling/resampling (as in sampling the continuous variable $t$ at a rate of $nT_s$. But this doesn't seem to explain it's use in $(\leftrightarrow 1)^n$ (for lack of a better symbol). I'm hoping someone here might be able to shed some insight.

Here's a link to the paper for further examples:
https://ccrma.stanford.edu/~stilti/papers/blit.pdf


texts: Is the system $y(n)=x(n^2)$ time invariant?

texts:
According to my solution,

$$y(n) = x(n^2) \tag{1}$$

for $y(n,k)$ i.e output for delayed input $x(n-k)$

$$x_2(n) = x(n-k)$$

so $$y_2(n) = x_2(n) = x\big((n-k)^2\big)$$

and for delayed output signal $y_1(n)$, replace $n$ by $n-k$ in equation (1), so we get,

$$y_1(n) = x\big((n-k)^2\big)$$ and therefore system is time invariant

But in the answers to the book in which this question it says the system is time variant. Can anyone point out the mistake in my steps, and give correct way of solving this sum?


texts: Understanding Eigenface Face Recognition Algorithm Notations

texts:
I have been reading the paper from Matthew A. Tusk and Alex Pentland

https://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf

I got stuck in this part. I have a couple of questions.

Assume my images are 256 by 256.

1- What are G1,G2...Gm ? Row vectors? Column vectors? (1*65536) or (65536 * 1) ? Or are they 2D vectors (256*256)?

2- Then it says A = [ F1 F2 ... Fm] . Now I do not get this notation. Again A is what ? A matrix of matrices? or just a huge column/row vector?

Thanks


texts: The output of a Digital to Analog converter has a magnitude roll-off in frequency as a Sinc function with the first null at the sampling rate. Why is this the case?


texts: Creating highpass from lowpass in Python

texts:
I want to create high pass filter from low pass filter in Python. My code:

h_lowpass = lp_design_window(fc, N, window)
dirac_delta = np.concatenate((np.ones(1), np.zeros(h_lowpass.shape[0] - 1)))
h_highpass = dirac_delta - h_lowpass

where lp_design_window is the function with arguments: cutoff frequency, number of coefficients and window function - and returns coefficients of low pass filter using window method. It attenuates desired frequencies, but not enough. Is it something wrong with dirac_delta? Or maybe I should search for the problem somewhere else?


texts: Understanding of Random Process/Random Variable

texts:
At a simpler level to my previous question, I wanted to confirm my understanding on Random Process based on Random Variables using an example.

So, I took this example:
If we consider a dice, which has six sides denoted by fi (i can be 1,2,3,4,5,6). So, assigning outcomes f1=10,f2=20,f3=30,f4=40,f5=50,f6=60 thereby constructs a Random Variable X.

Now, if we assign the outcomes: f1=15,f2=25,f3=35,f4=47,f5=80,f6=90 at a different time instant t2 and at time instant t3, we assign the outcomes: f1=43,f2=32,f3=16,f4=89,f5=99,f6=56.

So, if I plot this variation on a Graph, I am getting something like this below:

Now does this set of sample functions (Ensemble) constitute a Random Process?


texts: What are some of the reasons why DC offset is considered undesirable for a timeseries signal other than wastage of power during transmission ? Specifically, How does it cause spectral leakage and out of band emissions ?


texts: Implementing a 1-D Kalman Filter Regression, Missing the smoothing action (getting the opposite)

texts:
I am implementing the 1D Kalman Filter in Python on a fundamentally noisy set of measurement data, and I should be observing a large amount of smoothing...but, instead, my Kalman Filter is doing the exact opposite of smoothing, so that my filtered results look like this:

And since my "model" is just the simple trend formula, $$x_{k} = 3*(x_{k-1} - x_{k-2}) + x_{k-3} $$
the model comes up with outrageous predictions when the expected smoothing doesn't happen:

I am working off an incompletely specified 1D Kalman Filter implementation, so I am having to fill in the gaps...but I should be seeing much softer filtered values, with the peaks and valleys being outliers, with my naive model.  Also, my predictions should actually be more modest than the actual data (which would follow from a smoothed filtered collection).

I am trying to figure out what I am missing.  Here is the pseudo code for what I have implemented so far:

 Step 1: Minimize T = -log(Q/R) in order to discover Q
         Here I basically run the kalman filter with given T values,
         calculating Q from the Residual std

 Step 2: Run Kalman Filter with T held constant

 xprior = [x1,x2,x3] (calculated from the first three values of data)

 #prediction
 x_{k|k-1} = model_estimate(xprior)

 #add measurement
 innovation = data[k] - x_{k|k-1}
 K = R / (R + Q)

 #filtered value
 x_{k|k} = x_{k|k-1} + K*(innovation)

 xprior <- update_xprior(x_{k|k})

 residuals <- update_residuals(measurement,x_{k|k})

 R = variance(residuals)
 Q = e^(-T)*R

 ...

I must be missing something...Also, I am not confident that what I am doing with T and the way I am calculating Q is correct...but I should be seeing smoothing regardless, since that is what they did in the document I am working off of.


texts: Relation of roots of a polynomial to amplitude distribution of its coefficients

texts:
Consider $g(x)=\sum_{k=0}^{p}a_kx^{k}$, and let $\{\tau_1,\tau_2,\cdots,\tau_p\}\in\mathbb{C}$ be the roots of $g(x)$. What is the relation between how $\{a_0,a_1,\cdots,a_p\}$ are distributed (in amplitude values) and the magnitude of roots of $g(x)$?
Is there a condition on amplitude distribution of $a_k$s and degree $p$ which can aid in making all the roots to lie inside the unit circle.

An example: Below plot shows roots of 1000 polynomials each for degree $p = 5,10,20,50,100$ in $g(x)$ above {depicted in 'r.','g.','b.','k.' colours, respectively}. In each trial p+1 random numbers are drawn from Normal distribution $\mathcal{N}(0,1)$. [Right sub-plot]: The random numbers are assigned as $a_k$s. [Left sub-plot]: The random numbers are sorted in descending order and assigned as $a_k$s such that $a_{k-1}<a_k$. In the Left sub-plot, roots tend to move towards unit circle specifically for higher degrees.


texts: Why discretize a continuous transition matrix in Kalman Filter?

texts:
In the Kalman filter toolbox at http://becs.aalto.fi/en/research/bayes/ekfukf/cwpa_demo.html the example code shows that a function lti_disc is called, which is essentially a matrix exponential function. What is the purpose of such matrix exponentiation? Given a really simple transition model such as, e.g.

F = [ 1 0 ;
      0 1 ] ;

and doing the matrix exponentiation with lti_disc will always give a different transition matrix A, depending on the values of the Q process noise covariance, which is also an input to the lti_disc function.


texts: should this be viewed as an impulse response or step response

texts:
I'm trying to teach myself the relation between simple discrete ODE's and the impulse response-step response concept.

Getting back to the question: I don't expect anyone to read the whole thing but I'm confused early on which is a good thing.

Right at the beginning of the link, the author provides an example of a leaky integrator model in which the water facet is turned on and the water goes into the bucket and there is no hole in the bucket.

Given the leaky integrator model, the author then shows how, if there is
no hole in the bucket, then by setting k = infinity in the model, one can model the water level for the bucket. He then provides

a plot of the input: the water going into the bucket and
a plot of the output: the water level of the bucket.

My confusion stems from the plot of the input. I think that the plot of the input, rather than being a spike as is shown at time t = 1, should instead be a step response where there are spikes all the way from zero to one ? so basically a step response ?

My thinking is that, if the water is turned out at time zero, then the spike for the input  is at t = 0 and remains a spike until time t = 1 ?

The output plot seems fine to me.

If I am correct about the input plot being incorrect as is, then it kind of tells me that I understand the basic impulse response-step response concept. If I'm not, then, even though I loved the example, I'm still not quite getting it.

I asked the creator of the link in an email and he-she hasn't gotten back to me yet and I would like to understand this and don't honestly know if I'll get a response ( no pun intended ).


texts: How to find spectral density of a signal whose correlation depends on time?

texts:
Suppose we have two uncorrelated, real white noise $\xi_1(t)$ and $\xi_2(t)$, with zero mean, and equal variance.

Now, I create the signal $\xi(t) = a\xi_1(t)\cos(\omega_0 t) + b\xi_2(t) \sin(\omega_0 t)$. ($a,b$ are real)

I want to find the spectral density of $\xi$.
The correlation $<\xi(t) \xi(t+\tau)>$ depends on both $t$ and $\tau$. Now, to calculate spectral density, should I calculate
$\int_{-\infty}^{\infty} e^{i\omega\tau} <\xi(0) \xi(\tau)> d\tau$?

Or should it be $\int_{-\infty}^{\infty} e^{i\omega\tau} (\text{average over all $t$ of} <\xi(t) \xi(t+\tau)>) d\tau$ ?

I have calculated the correlation $<\xi(t) \xi(t+\tau)>$

$=<[a\xi_1(t)\cos(\omega_0 t) + b\xi_2(t)\sin(\omega_0 t)]\times [a\xi_1(t+\tau)\cos(\omega_0 (t+\tau)) + b\xi_2(t+\tau))\sin(\omega_0 (t+\tau)]>$

$= a^2 \cos(\omega_0 t)\cos(\omega_0 (t+\tau))<\xi_1(t) \xi_1(t+\tau)>$

$+ b^2 \sin(\omega_0 t)\sin(\omega_0 (t+\tau))<\xi_2(t) \xi_2(t+\tau)>$

$+ ab \cos(\omega_0 t)\sin(\omega_0 (t+\tau))<\xi_1(t) \xi_2(t+\tau)>$

$+ ab \sin(\omega_0 t)\cos(\omega_0 (t+\tau))<\xi_2(t) \xi_1(t+\tau)>$

Since $\xi_1(t),\xi_2(t)$ are uncorrelated noise, this is equal to $= a^2 \cos(\omega_0 t)\cos(\omega_0 (t+\tau)) \cdot c\delta(\tau) + b^2 \sin(\omega_0 t)\sin(\omega_0 (t+\tau)) \cdot c\delta(\tau)$

(where $\delta(\tau)$ is the Dirac delta function, $c$ is a real number that depends on the amplitude of the noises)


texts:
I have to do two frequency shifts:

the first: to transform a lowpass filter into a highpass filter
I've implemented it in this way:

h_hpf=h_lpf.*exp(j*pi*(0:N))

the second: to transform a lowpass filter into a bandpass filter centered at $\pi/2$.
I've implemented it in this way:

h_bpf=h_lpf.*exp(j*pi/2*(0:N))
I used the fvtool() command to show frequency responses, and I've observed that the frequency response of h_hpf is also replicated in the negative frequency axis, while it doesn't happen for the frequency response of h_bpf.
Indeed, the frequency response of h_bpf appears only in the positive frequency axis, and the replica in the negative axis doesn't appear.
Moreover, its impulse response h_hbf has complex values, so I think that the correct shift is:
h_bpf=h_lpf.*(exp(j*pi/2*(0:N))+exp(-j*pi/2*(0:N)))
Which is the correct shift?


texts:
I try to understand audio signal processing from Deep Speech: Scaling up end-to-end speech recognition paper.
Quote:

Since the Switchboard and Fisher corpora are distributed at a sample rate of 8kHz, we compute spectrograms of 80 linearly spaced log filter banks and an energy term. The filter banks are computed over windows of 20ms strided by 10ms.

Am I calculating window_size correctly?
Is stride the same as hop_size between windows?
I'm not sure what does it mean

we compute spectrograms of 80 linearly spaced log filter banks

Example:
I took the audio from first sample data of The LJ Speech Dataset
:
import librosa

audiowaveform, sample_rate = librosa.load('Downloads/LJ037-0171.wav')
audiowaveform.shape, sample_rate
((167226,), 22050)

win_len = 20  # in milliseconds
n_fft = int(sample_rate * win_len / 1000) # 441
hop_size = 10  # in milliseconds
hop_length = int(sample_rate * hop_size / 1000)  # 220
librosa.stft(s, n_fft, hop_length).shape
(221, 761)
```


texts: cceps (complex cepstrum) function in Matlab

texts:
My doubt is in regards to the usage of the function cceps in Matlab.

The function is giving me the cepstrum for the entire speech signal of 3 or 4sec. This works well when I play back the signal after taking the inverse (using theicceps function). However, in order to analyze the given speech signal, we are supposed to only consider 20ms frame.

Is taking the complex cepstrum throughout the signal valid? Why is it working out here?


texts: Why is mean-square summation of the area under the curve?

texts:
Would someone please kindly help me to nail this "simple" issue? My input is a matrix of power PSD and frequency. I need to get RMS and I read this
http://vibrationresearch.com/wp-content/uploads/2018/04/RandomRMSCalculation.pdf
In page 2, under equation 5, it says "The total area under the curve will then be the sum of the individual area calculations between each pair of breakpoints, and this sum is the mean-square acceleration. We take the square-root of the result to get the RMS acceleration level".
Q1: Why the summation of the area under the curve is mean-square? I think it should be divided by the number of frequency samples. Am I missing something?
Q2: I tried with rms function in matlab, rms(matrix,2) and it returns the results much lower than following the instruction from the link. Is this way correct?
Thanks a lot.


texts: How to calculate bandwidth require for OFDM

texts:
I read LTE PHY profile and do not understand the relationship between number of FFT and bandwidth. For example LTE downlink channel bandwidth 1.25MHz and 5MHz have 128 FFTs and 512 FFTs accordingly. Another example is 802.11 a bandwidth 20MHz, FFT size 52.

I design an OFDM system system at carrier frequency 910MHz, with FFT size for example 64, how much coherent bandwidth that the system need? What is the formula? What other factors should be taken into account?

Add >> 18 Aug (I don't have enough reputaiton to add comment on Jason R reply)

When I test a simple software defined radio OFDM system with parameters as FFT = 256, data carrier per symbol = 125, number of pilot = 25, bandwidth = 5MHz, carrier frequency = 910MHz, the channel estimation wouks fine. If I change bandwidth to 1MHz, then the channel estimation goes wrong. In this case I reduce sampling frequency 5 times, so carrier spacing is reduced 5 times too.

I think the carrier spacing must be greater than a limit so that intersymbol interference (or inter-carrier interference?) is small enough to not cause channel estimation failture. How to find the carrier spacing limit?

Or delay spread $\tau \ll T$ where T is symbol interval may cause the problem. How to calculate the limit $\tau$?

Thanks


texts: From Delay Spread to Coherence Bandwidth?

texts:
The phase difference of two sinusoids is given as:
|  TX --(r)--->RX <------- (Reflection)|
d: whole path from sender to reflecting object
$ \Delta \theta = 2 \pi \cdot \frac{(2d - r)-r}{c} f + \pi $
The coherence bandwidth is given by: $ W_c := \frac{1}{(T_d)}$
$ W_c $ is $ f_1 - f_0 $, so basically I should somehow calculate that from the phase difference formula? But I'm missing the beginning.


texts: Do FFT-based filtering methods add intrinsic latency to a real time algorithm?

texts:
In a current audio processing application, I am working entirely in the time domain using IIR filters. This is so I can use an audio buffer of just 3 or 4 samples and can guarantee that there will be some form of audio output within that frame.

If I were to switch to FIR based filters (to reduce filter ringing time as outlined in this question), I would probably need to use an FFT-based overlap-add technique to maintain processing efficiency. As I see it, I would have to fill an FFT frame buffer with samples before I would be able to do the transforms. This would add an intrinsic latency to the algorithm proportional to the number of samples used to calculate the FFT. Am I missing a trick?


texts: Even more on: Kalman filter for position and velocity

texts:
As a follow-on to this great question...

I'm trying to implement the model in the referenced article (and referenced Wikipedia article), but I'm not getting the results I think I should be getting.  Ultimate I would like to take this example and marry it with intermittent GPS data so as to simulate accelerometer failover in situations where GPS is lost.

I've build the following model in Simulink:

UPDATED MODEL:

and power it using the following UPDATED Matlab code:

% Configure simulation
sim_length=10;
Ts=0.1;

% State space equations
A = [1 Ts;0 1];
C = [1 0];

% Noise characteristics
sigma_a = .1;
process_noise_power = 10;
measurement_noise_power = 10;

% Noise matrices
G = [(Ts^2)/2; Ts];
H = [1];

% Process noise
Q = sigma_a^2;
% Q = G*G'*sigma_a^2;

% Measurement noise
R = measurement_noise_power;

% Steps (not needed)
N = 0;

% Run simulation
sim('full_discrete_2inputs');

% Show results
close all;
figure(1);
hold on;
plot(tout,x_actual,'b-');
plot(tout(1:length(y_measure)),y_measure,'m-');
plot(tout(1:length(x_est)),x_est(:,1),'r-');
legend('True Position','Measured Position','Filtered Position','Location','Best');
title('Actual vs Kalman Filtered Position');
xlabel('Sample time [s]');
ylabel('Distance from origin [m]');
grid on;
hold off;
shg;

And here are my UPDATED results:

So with that, I have the following questions:

Did I inject the noise at the correct location? Is it appropriate to inject it prior to the double integration?
What's a good real-world noise power level to use in this application?  I picked a number at random, but I imagine whatever I pick has to be considered in my Q or R matrices.  I know if I bump the noise up to something like 50, the filter output and measurement value go wild, but still track each other.
To get the simulation to work I had to change my H matrix to 1, but [1 0] was specified in the example.  Does anyone know why?


texts: Why is the total signal response response of zero input + zero state, $y(t) = y_0(t) + h(t)\star x(t)$, not an LTI?

texts: For the zero input + zero state response in continuous time linearly time-invariant systems, why is the $y(t)$ equation not "technically" considered an LTI? I read this in a journal and there was no justification. How would you change the equation to make it a true LTI?


texts: Time delay estimation of two sensor signals

texts:
I have two sensors which are placed around 30  to 40 cm apart , both of them measure the same   signal say 'lambda signal'. so there is a delay in the measurement of the signal from the second sensor. During the travel there is some noise added to the signal. What are the procedure that I can use to calculate the time delay between them.
According to my knowledge if I cross correlate both signals, the peak gives me the delay value of the second, but this is not working. Is there any other methods to find the time delay between eachother. Please guide me as I am new to signal processing.![Part of the signals][1]![Another part signal][2]

link: https://docs.google.com/open?id=0BzT9Q3HfpmDKZERPSGJoTEFyTzA
link : https://docs.google.com/open?id=0BzT9Q3HfpmDKeVhMSnZLa3BfMjQ
Thanks in advance.
Naveen


texts:
I have a doubt regarding multi track/multi layer recording of audio signals. in recording studio each instrument or chorus,playback will use a independent layer/track..so after recording a all tracks ..how do they mix all of those into a single track ..any mixing scheme?..or it is the simple time division multiplexing??

I am asking this because while I read any mp3 or wav file in Matlab..I can see only a one channel or 2 channel signal...

any help would be appreciated...


texts:
I am trying to make my own demonstration in order to find what is equal to the FT of the Autocorrelation function.
The autocorrelation function in some book about  turbulence is defined as:

r($\tau$)=$\lim_{T \to \infty}\frac{1}{T}  \int_{t0}^{t0+T} u(t)u(t+\tau)dt$  (1)

where T is the time interval, t0 is the starting time of measurement of the parameter. The power spectral density in this case is defined:

S(f)=$\frac{1}{2\pi} \int_{-\infty}^{\infty} exp(-i\tau 2 \pi f)r(\tau) d \tau $  (2)

If we substitute (1) in (2) we obtain:

S(f)=$\frac{1}{2\pi} \int_{-\infty}^{\infty} exp(-i\tau 2 \pi f) [\lim_{T \to \infty}\frac{1}{T}  \int_{t0}^{t0+T} u(t)u(t+\tau)dt] d \tau $  (3)

Now I would like to do a similar demonstration like this the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal , but there are some differences and unconsistency, that I would like to help me to fix it (because I am getting mad trying to  demonstrate it)-
First of all, my autocorrelation definition is different and it hasn’t an integration from -$\infty$ to +-$\infty$ , and furthermore, in my definition of autocorrelation is divided by the total time interval.
So do you think that I should, redefine the integral domain in the power spectral density? If yes how?

Anyway I tried to continue this “demonstration” in order to explain all my doubts (maybe I made a mistake, but I hope that you can help me to fix my error).
I rearranged (3) in this way

S(f)= $lim_{T \to \infty} \frac{1}{T} \int_{-\infty}^{\infty}  u(t) [\frac{1}{2\pi} \int_{t0}^{t0+T}u(t+\tau)exp(-i\tau 2 \pi f) d \tau] dt $  (4)

I am not sure about that, but I thought that the part in the square brackets can be considered, as happend in the case in the link the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal, to be:

$\frac{1}{2\pi} \int_{\color{red}t\color{red}0}^{\color{red}t\color{red}0\color{red}+\color{red}T}u(t+\tau)exp(-i\tau 2 \pi f)d \tau =U (f) exp(i 2 \pi ft)$   (5)

where U(f) is the fourier transform of u(t). An other doubt about the equation (5) is the fact that the coefficient $\frac{1}{2\pi}$ does not appear in the demonstration of the book that I cited in the link.

If we substitute (5) in (4) we obtain:

S(f)= $lim_{T \to \infty}\frac{1}{T} \int_{-\infty}^{\infty}  u(t) [U (f) exp(i 2 \pi ft)] dt $  (6)

If you rearranged this:

S(f)= $lim_{T \to \infty}\frac{1}{T} U (f) [\int_{-\infty}^{\infty}  u(t) exp(i 2 \pi ft)] dt $  (7)

if we consider that the part in the square bracket is equal to U(f)* which is the conjugate of U (f):

S(f)= $lim_{T \to \infty}\frac{1}{T} U (f) U (f)*$ (8)

which is equal

S(f)= $lim_{T \to \infty}\frac{1}{T} |U (f)|^2$ (9)

Now in this demonstration there are a lot of uncertainities and "bug". I hope you can kindly explain in what I am wrong and which could be the best strategy to reach the right result.

If you find references of book where similar demonstration are done with this definition of autocorrelation, it will be amazing
 I need a formulation without citing the expected value and the $\delta (t-\tau -t')$


texts:
After studying z transform from different books and literature on internet I want to ask few which makes me confuse.

a) From the Discrete Time Fourier Transform we have drive equation for z transform. $$ X(z)= \sum _ {n=-\infty}^{+\infty} x[n]z^{-n}$$ where $z$ is represented in polar form $z=re^{j\omega}$
I want to know that why we represent $z$ in polar coordinates? as in some books it is written than $z$ is complex $z=\sigma + j \omega$

b) ROC for $z$ transform is same as Laplace Transform? In laplace transform we check that direction of $t$ (i.e. if we have $u(t)$ than the $Re[s]  > a$)?


texts: Is noise figure dependent on input noise power?

texts:
I was reading about Noise Figure on Wikipedia (https://en.wikipedia.org/wiki/Noise_figure). There, I saw the following definition:

  The noise factor is thus the ratio of actual output noise to that which would remain if the device itself did not introduce noise, or the ratio of input SNR to output SNR.

where noise figure is $10\log_{10} \text{(Noise Factor)}$.

That arises a question in me.

  If Noise Figure is given for an amplifier as a constant (for example Noise-Figure=$5.5$dB), since amplifier noise is additive (i.e. its power adds up to the input noise power), then the amplifier noise must be dependent on input noise power, which contradicts with my intuition.

As an example, is an amplifier has gain of $1$ (i.e. no amplification, no attenuation) with a NF of $5.5$dB, it preserves the input signal power, in a sense that $$P_{\text{out}}=G\cdot P_{\text{in}}=P_{\text{in}}$$and $$\sigma^2_{\text{noise,out}}{=G\cdot\sigma^2_{\text{noise,in}}+\sigma^2_{\text{amplifier}}\\=\sigma^2_{\text{noise,in}}+\sigma^2_{\text{amplifier}}}$$
hence $$\text{Noise Factor}={\sigma^2_{\text{noise,out}}\over \sigma^2_{\text{noise,in}}}=1+{\sigma^2_{\text{amplifier}}\over \sigma^2_{\text{noise,in}}}=10^{0.55}=3.55$$and $$\sigma^2_{\text{amplifier}}=2.55\sigma^2_{\text{noise,in}}\qquad\qquad(!!)$$

However, we know that Noise Figure is highly used in practice. Where is my mistake in understanding it?


texts: Wiener filter yielding negative MSE (and MMSE)

texts:
I implement a Wiener filter using the following code

function [h, yy, mse] = wiener_h(y,x,M)
yy=xcorr(y(1:M));
yy=toeplitz(yy((length(yy)+1)/2:((length(yy)+1)/2)+M-1));
debug = sprintf('Inverting Autocorrelation matrix..')
yy=inv(yy);
%%
debug = sprintf('Crosscorrelation vector')
yx = xcorr(y,x);
yx = yx((length(yx)+1)/2:((length(yx)+1)/2)+M-1); %consider only non-negative lags
debug = sprintf('Optimal Weiner impulse response..')
h =yy*yx;
mmse = var(x)-h'*yy*h;
mse = var(x)-2*yx'*h+h'*yy*h;

where y is the input signal, x the desired output, M number of coefficients, h the optimal (wiener  values) and mse the mean square error.
The input to the filter is are magnitude values of a pulse (NR, SV, SA)
The desired output is TX

The output of

[h, yy, mse] = wiener_h(y,x,15000)
plot(t,filter(h,1,y))

where the number of coefficients M is equal to the number of samples N, is mse = -1.14e+03
Whereas the output of

 [h, yy, mse] = wiener_h(y,x,1500)
 plot(t,filter(h,1,y))

Where the number of coefficients is 10% of N, is mse = -7.68e+04
The plotted results of the two above (left M=N, right M=0.1N)

Can someone help me understand what's going on please?


texts: Least Squares Fitting to Inverse Exponential Function

texts:
I have a time series of measurements that resembles the shape of an exponential function. The samples are a bit noisy and sometimes there is a weak sine like ripple signal ontop of it.

Simplified the function looks like this:

$$ y = a (1-e^{-bt + c}) + d $$

The data that I have lies between $t=$0 to 2, e.g. where the function is still curvy.

My goal is to find $a+d$, e.g. the value of the asymptote or the $y$ value of the graph above at infinite time.

My idea so far is to take the samples that I have, do least square fitting on the function above and then just solve it.

Q: Will this approach work or is there a easier way?

Q: I'm a bit lost how to do a least square fitting to a function like the one above. With a falling decaying exponential I would just take the log and least square fit a line, but how do I do this with my rising function?

Oh, one more thing: I don't need to do this as a one-off so I can't just plug the data into a general purpose solver. I have to implement this later for resource constrained microcontroller.


texts: Relating phase difference in a signal with its frequency

texts:
I have a 134 point complex signal. Using MATLAB, I have observed that the phase difference between any two adjacent samples of this signal is 2 radians i.e. the phase increases linearly from 0 to 268 radians over 134 samples.

How can I relate this phase difference between adjacent samples with the frequency of the same signal?

PS I am sampling at 422 Hz.


texts: How to construct an unknown signal from sparse periodic samples?

texts:
I have an unknown black-box that is generating a deterministic continuous periodic signal and the period is variable between 1 to 30 seconds.

When the black box is turned on, the period is fixed, but unknown.
Sampling rate and signal period are not synched, and it is highly unlikely to be integer coefficient of one another.

All that is available to me is sparse (contains missed) samples with period of 100 ms to 10 ms.  (The parameter of source is fixed every time, for example I should construct a 3s signal from periodic 10ms samples). There is no limit on the duration of sample gathering but samples may be missed. Signal is smooth and does not change fast and slowest sampling rate satisfies nyquist rate.

Is there a method to construct a single period of signal and enhance its resolution as time goes by?


texts: Whether $y(t) = x(t/3)$ is a memoryless system?

texts: I am a newbie in signal and system, and feeling confused about the memoryless property, it says at the definition, if a system gives output only depends on present input, then its memoryless. So if at $t=3$, $x(1)$ is needed, whether this count as a memoryless system? Or it only consider the time shift like $y(t) = x(t-1)$?


texts:
I wonder if anyone can point me to services or dev libraries that focus on or suit phoneme recognition? Target audience is children at learning to read stage.

We have trialed IBM Watson (not good) and Google (only a little better).

Hope this is a suitable question for this forum. Thanks in advance.