texts:
Let's say I have an audio clip of $N$ samples.  I want to describe every $10$ ms of audio with an FFT of size $1024$, without information bleeding over from other $10$ ms windows.  Given that the sample rate of this audio clip is a standard $44100$, I can't do that with a regular FFT.  A regular FFT would give me $441 / 2 = 220$ bins.

Is it possible to take many, many FFTs of differing sizes over different windows to infer a larger frequency resolution for a smaller timeframe?  Let's say I take a $1024$-FFT for each sample range increasing by $t = 1$ (samples $0$ through $2047$, then samples $1$ through $2048$ and so on), then a bunch of $512$-FFTs, then a bunch of $256$-FFTs.  Would this large dataset contain enough useful information to "subdivide" the bins of the smaller, $220$-FFT?


texts: Image processing research areas using GPU

texts:
I'm an undergraduate student and I am currently trying to find suitable topic for my bachelor thesis. The broader topic I've selected is "Image processing on GPU", so basically I'll be developing image processing algorithms specifically for GPU.

The problem is, I need to do a novel approach, not reinvent a wheel.  There's so much going on in image processing, lately on GPUs as well, that I'm having trouble finding a suitable area that hasn't been researched in this way and which isn't above my capabilities as an undergraduate.

The question is, can you point me in a direction what are current research areas in image processing that could benefit from GPU implementation? I'd love to do something for medical or astrophysical applications, but I'm open to other areas as well.

Thanks!


texts: Comparing with FIR, why does IIR always gives a bias error, when calculating signal DC?

texts:
I am trying the DC blocking function in matlab. The algorithm of DC blocking is to find DC of the signal first and then remove the DC value from the signals. Here are the codes I tried, which are sample codes from Matlab:

  t         = (0:0.001:10000)';
  x         = sin(30*pi*t) + 1;

  hDCBlock1 = dsp.DCBlocker('Order', 4);  % IIR is default
  hDCBlock2 = dsp.DCBlocker('Algorithm', 'FIR', 'Length', 100);
  hDCBlock3 = dsp.DCBlocker('Algorithm', 'Subtract mean');
  for idx = 1 : 10
    range = (1:1000) + 1000*(idx-1);
    y1 = step(hDCBlock1, x(range));
    y2 = step(hDCBlock2, x(range));
    y3 = step(hDCBlock3, x(range));
  end
  plot(t(1:1000),x(1:1000),...
    t(1:1000),y1, ...
    t(1:1000),y2, ...
    t(1:1000),y3);
  legend(sprintf('Input DC:%.3f',    mean(x)), ...
    sprintf('IIR DC:%.3f',           mean(y1)), ...
    sprintf('FIR DC:%.3f',           mean(y2)), ...
    sprintf('Subtract mean DC:%.3f', mean(y3)));

I tried to adjust the length of data and order of the IIR, the DC given by IIR always returns a bias error. However the FIR and mean does not, or gives very little errors.

Could anyone explain on this?


texts: I am doing my Bachelor's project on Road extraction. I have detected the line segments(edges) in the image using hough transform. I need to detect the curved lines in the road and link them(edge linking). help me with matlab coding.


texts: What would be the most significant difference when using either a sobel filter or a gaussian-derivative filter, apart from the obvious difference of the size of the 2 filters?


texts: Image processing: High pass filter to remove lowest 16 frequencies

texts: How can I remove the lowest n frequencies of an image? Are the lower frequencies around the centre of the FFT result?


texts: Implementing Optimal Crosstalk Cancellation with matlab

texts:
I am trying to implement a crosstalk cancellation in Matlab, more precisely part of : Optimal Crosstalk Cancellation for Binaural Audio with Two Loudspeakers (http://www.princeton.edu/3D3A/Publications/BACCHPaperV4d.pdf).

I am using the matrix H from equation (24) in the paper as my filter and use it as in equation (11). The output does not sound good and has lost all previous "binaural features". Here is a plot of the original binaural signal and the X talk cancelled one :
www.hostingpics.net/viewer.php?id=196157img.png

The binaural sound is a 1000Hz sound wav at the left ear (created with the CIPIC HRTF library). I would have expected the X talk signal to be kind of like the binaural one (where the right channel should be attenuated) but its not the case.

Here is the actual code :

%Create binaural sound
aSound = wavread( '1000hz.wav' );
binauralSound = CreateBinauralSound( aSound, pathToHRTF, 1, 9 ); % -90 degree, 0 elevation

%NFFT = 2^nextpow2( length( binauralSound ) );
NFFT = length( binauralSound );

binauralSoundF = fft( binauralSound, NFFT );

theta = 15; %Angle from nose to right speaker in degree
deltaR = 0.15; %Distance between the ears
deltaS = 2; %distance from  ear plane to speaker plane (going through nose);

%Compute ear and speaker position
[speaker1x, speaker1y, speaker1z] = sph2cart( deg2rad(90 + theta), 0, deltaS );
[speaker2x, speaker2y, speaker2z] = sph2cart( deg2rad(90 - theta), 0, deltaS );
lEarx = -1 * deltaR / 2;
rEarx = deltaR / 2;

%Compute distances
l1 =  sqrt( ( speaker1x - lEarx )^2 + ( speaker1y - 0 )^2 );
l2 = sqrt( ( speaker1x - rEarx )^2 + ( speaker1y - 0 )^2 );
deltaL = l2 - l1;
TC = deltaL / 340;
g = l1/l2;

%Filters and parameters
beta = 0.5;
H = zeros( 2, 2, NFFT );
alpha = zeros( NFFT, 1 );
final = zeros( NFFT, 2 );

%Compute value for each freq bins
for k = 1 : NFFT
    omega = exp( -2 * pi / k );
    %omega = k;
    H( 1, 1, k ) = ( g^2 * exp( i * 4 * omega * TC ) - ( beta + 1 ) * exp( i * 2 * omega * TC ) ) / ( g^2 * exp( i * 4 * omega * TC ) + g^2 - ( (g^2 + beta)^2 + 2 * beta + 1 ) );
    H( 2, 2, k ) = H( 1, 1, k );

    H( 1, 2, k ) = ( g * exp( i * omega * TC ) - g * ( g^2 + beta ) * exp( i * 3 * omega * TC ) ) / ( g^2 * exp( i * 4 * omega * TC ) + g^2 - ( (g^2 + beta)^2 + 2 * beta + 1 ) );
    H( 2, 1, k ) = H( 1, 2, k );

    alpha( k ) = exp( -1i * omega * l1 / 340 )/ l1;

    final( k, : ) = alpha( k ) * binauralSoundF( k, : ) * H( :, :, k );
end

resultSound = ifft( final, NFFT, 'symmetric');
resultSound( :, 1 ) = resultSound( :, 1 ) / max( abs( resultSound( :, 1 ) ) );
resultSound( :, 2 ) = resultSound( :, 2 ) / max( abs( resultSound( :, 2 ) ) );

I am really not sure I am computing H the right way nor I am sure I am applying it correctly. There is not a lot of implementation of crosstalk cancellation on the web so I find it quite hard to find my mistake. I am guessing its something really obvious for all the dsp gurus out there. Your help would be greatly appreciated !

Thanks in advance.


texts: How can I graph the amplitude for a given frequency window over time?

texts:
I have simple audio files of single plucked notes from a guitar:

I want to graph the amplitude of each partials over time, rather than just have this color-coded "dB" intensity graphing I have here.

I can delete the extraneous frequencies in Izotope RX easily enough which I think does the same thing as FFT windowing:

Is this the right approach? How can I next get the amplitude of this window over time?

I know basic C++/JUCE but mostly for synth design. I'm not a pro on DSP, FFT, etc.

I appreciate any guidance.


texts:
I am using emgucv with opencv.

I have an object I wish to detect based on colour.
I convert the image to Hsl colour space and perform a filter on it.

The object I am detecting is perceived as yellow so I cover that range and make sure I am using/filtering good saturation.

This has been working well.

But I have now upgrade my usb camera and I have changed the lighting in the room.

Now running the same code I can it is ‘missing’ a lot of the object.

So, I inspect the hue of the pixels that have been missed out.

I can now see that the hues of those pixels have ‘shifted’ to the red range.  This incidentally is more noticeable when the object is closer to the camera.

So my question is this;

Is yellow an unreliable colour to detect compared to say red, green, blue etc?

Is the object not truly yellow?  Though when the object is further away it does not appear to be a problem,

Is colour segmentation unreliable?

Should I consider using a different colour space?

I will post my existing code if people require but I do not think it is required?

Thanks


texts: Kalman filter equations/code to estimate velocity from known position and acceleration

texts:
I was reading this post and I have the same problem: I want to estimate the velocity using displacement and acceleration measures.

But the problem is that I don't know how I have to code/implement an algorithm using the matrices of Kalman filter presented in the post because there isn't equations with the matrices H, A, P, Q, R, ... And also the nomenclature from one post to other is different... It is the first time I need to use a Kalman filter, so I am not familiar with the use of this filter.

% data(:,1) Time
ts = data(2,1)-data(1,1);
ts2 = ts^2/2;
ts3 = ts^3/6;

nsamples = length(data);
x0 = mean(data(:,3));

ak = data(:,2)./(0.0991/9.81);       % Acceleration
dk = (data(:,3)-x0)./0.1266*2/1000;  % Displacement

% EXAMPLE MATRICES
HH = [1 0 0 0;
      0 0 1 0];

RR = [0.0100     0;
        0     0.0025];

AA = [1 ts ts2  ts3;
      0  1  ts  ts2;
      0  0   1   ts;
      0  0   0    1];

PP = [0 0 0   0;
      0 0 0   0;
      0 0 0   0;
      0 0 0 0.001];

QQ = [0 0 0   0;
      0 0 0   0;
      0 0 0   0;
      0 0 0 0.001];

% Current state estimate
xk = [];
xk_prev = [dk(1); 0; ak(1); 0];

%        Disp           Velocity     Acc                   Jerk
xtrue = [ dk;  zeros(nsamples-1,1);   ak;   zeros(nsamples-1,1)];

% Buffers for later display
xk_buffer = zeros(4,nsamples);
xk_buffer(:,1) = xk_prev;
Z_buffer = zeros(1,nsamples);

for k = 1:nsamples-1
    Z = xtrue(k+1);
    Z_buffer(k+1) = Z;

    P1 = AA*PP*AA' + QQ;
    S = HH*P1*HH' + RR;

    K = P1*HH'/S; % inv(S)
    P = P1 - K*HH*P1;

    xk = AA*xk_prev + K*(Z-HH*AA*xk_prev);
    xk_buffer(:,k+1) = xk;

    % For the next iteration
    xk_prev = xk;
end

(the language program used is MATLAB)


texts: Convolution : Anechoic Speech signal and Impulse response

texts:
I am trying to perform convolution in MATLAB of an anechoic speech signal (2.1 s) and a impulse response (20 s) using the following code:
[y_anech,fsSig]=audioread('MA03_03.wav'); %speech signal fsSig = 48000

[y_IR,fsIR]=audioread('IR.wav'); %reads the Impulse response, fsR = 44100

%%% Resampling the speech signal
[P,Q] = rat(fsIR/fsSig);
y_anech_resampled = resample(y_anech, P,Q);

sig_conv_meas = cconv(y_anech_resampled,y_IR)/fsIR;
sig_conv_meas = sig_conv_meas/max(abs(sig_conv_meas));

player1 = audioplayer(sig_conv_meas,fsIR);
play(player1);

Running the above code, it gives me the required convolved signal, however at 20s, I also hear the anechoic speech signal. I have the following questions:

What the reason behind that I hear the original signal after the convolved signal? How should I remove it?

See the attached plots. In the last plot, around 20 s, there appears the anechoic signal:

ETC plot of IR:

When performing convolution, why do I need to divide the convolved signal with sampling frequency as I do in my code (sig_conv_meas = cconv(y_anech_resampled,y_IR)/fsIR)?

Thank you!


texts: hough voting, turning transformation into a vote

texts:
I need help understanding how/what does it mean to "turn" a transformation parameters into a vote.
I know how to recover the transformation parameters (in my case: 6 DOF), but how this turns into a vote in 6D (??) space ?

Thank You !


texts: Time series forecasting on dsp embedded systems

texts:
I have a time series analysis algorithm that uses arima modeling and fits a model to the time series and gets the residual values. It is more of a statistical signal processing algorithm than a digital one. I have a matlab code for it which uses functions from econometric toolbox.

I need to program it on an embedded system. My first instinct was to do this on dsp embedded system, but I couldn't find anything about parametric methods in the dsp embedded books. Most of the books only cover the FIR, IIR and adaptive filtering. Is it possible to do something like this on embedded systems?

and if yes, assuming that I use one of ARM cortex M series microControllers, is there any specific header file that I can use?

I hope my question is clear enough but I can explain it more.

Thanks,

Giselle

clarification:

I was trying to say that what I have is more into parametric estimation part of digital signal processing but I guess I said it the wrong way.

I have a machine learning algorithm (more specifically an outlier detection algorithm). my input is a time series data, lets say from an ADC and it is continuous but I am doing segmentation on the data. I think the segmentation will also help with the computational complexity and the limited RAM resources (If that was what you meant by estimating required RAM resources).

The limited resources of microcontrollers is my exact problem. I do have libraries in both C++ and python that will provide me with functions that I need for my code, but adding them all in the code will make it so huge. and that is why I thought maybe dsp systems are my best choice.


texts:
Suppose we use a uniform linear array for DOA estimation. Assuming neither of the signals and noise are correlated, and the noise power on each sensor is equal, it's well known that the covariance of the received signal can be decomposed as:
$$\boldsymbol{R}=\boldsymbol{A}\boldsymbol{P}\boldsymbol{A}^*+\sigma^2 \boldsymbol{I}=\boldsymbol{S}\boldsymbol{\varLambda_{\mathrm{s}}}\boldsymbol{S}^*+\sigma^2\boldsymbol{G}\boldsymbol{G}^*,$$
where the eigenvectors in $\boldsymbol{S}$ and $\boldsymbol{G}$ are selected to be orthogonal and normalized.
My question is how to prove the equality $$\boldsymbol{G}\boldsymbol{G}^*=\boldsymbol{I}-\boldsymbol{A}(\boldsymbol{A}^*\boldsymbol{A})^{-1}\boldsymbol{A}^*$$ in the paper Maximum Likelihood Methods for Direction-of-Arrival Estimation. It seems very basic but I can't come up with the proof.


texts: I want to do frequency analysis of audio data, basically trying to figure out what the notes are in a song algorithmically.  The standard approach is to decode the MP3 into PCM data and run it through an FFT.  However, notes below around middle C require too much precision for an FFT to work well.  Since the MP3 already has frequency information in small time increments, is it a) precise enough to identify frequencies within about 2 - 5 Hz, and b) is there any code out that there does this conveniently?  Most code (unsurprisingly) seems to convert it to the time domain, but I'd like to have just the frequency data.


texts: Interpretation of Clarke's Doppler power spectral density

texts:
What I understand of Doppler spread is that the relative motion between Transmitter (TX) and Receiver (RX) change the exposing time of signal. In rapport to a constant-distance TX-RX, a moving toward each other TX-RX "compresses" signal in time (signal takes less time to propagate), then signal is "expanded" in frequency domain. Similarly, a moving away RX-TX "expands" signal in time and "compresses" its spectrum. In short, that is scaling Fourier Transform. These two extreme cases set the left and right bounds of spreading an original frequency between $-f_d$  and $+f_d$ where $f_d$ is max Doppler spread.

In looking at the Clarke model, it is just multiple propagation model with rich scattering environment and equal angle of arrival. (link for more details
Clarke model)

If I understand well, there are two assumptions which are reasonale in urban environment:

Rayleigh fading
equal angle of arrival, or equal receiver sensitivity

I have followed the math from the original article, it seems ok. The final Doppler power spectrum is then
$\displaystyle S(f) = \frac{1}{\pi f_d \sqrt{1 - \left(\frac{f}{f_d}\right)^2}}$

What I don't understand is that why energy is concentrated to the two extreme spread frequency $-f_d$ and $f_d$ while angles of arrival are uniform. Is there any physical interpretation ? What am I missing from the famous Clarke model ? Personally, this model seems well-model the typical urban environment.

R. H. Clarke, A Statistical Theory of Mobile-Radio Reception, The Bell System Technical Journal, July/Aug 1968, p. 957ff

  Answers Although the answer of Carlos captures the most fundamental mathematical part, the real answer is in his comment about "mapping
  between angle and frequency". Moreover, the answer of Maximilian is
  interesting too.


texts: How accurate is Wavelet packet decomposition for emg signals with frequencies from 10 to 40 Hz?

texts: I have an electromyographic signal sampled at a frequency between 10 and 40 Hz in a time window of 300ms. Is wavelet packet decomposition adequate in order to perform features extraction for subsequent pattern recognition?


texts: Estimate mean amplitude in frequency range using fft

texts:
I would like to roughly estimate the mean amplitude over a frequency range (2-3hz). My impression was that I could transform the data using FFT, and just take the mean of the sub part of the output (absolute values).

However, the input I provide to numpy.fft.fft is in the range -1.5 to 2.5, while the output (absolute values) are in the range of about 0-6500. Do I need to scale the y-values somehow?

Any other suggestions on how to roughly estimate the mean amplitude of a frequency range?

Update:

The problem is that I am doing peak detection and want to discriminate out (as noise) all peaks that are smaller than 10% of the rest of the peaks. The majority of the peaks are expected to have rather equal amplitudes, and it is 10% of this amplitude that all peaks must at least be.


texts: Biquad Treble Shelf filter gain being weird

texts:
Trying to compute Biquad Coefficients for use in a codec chip TLV320AIC3100.

Bass Shelf(boost and attenuation) work as expected. Treble Shelf attenuation works as expected but Treble Shelf boost is being weird in the sense that the amplitude difference between non treble frequencies and treble frequencies is as expected (Say 10dB) but instead of the overall gain being +10dB at treble frequencies, it ends up being +3dB(ish) and non treble frequencies end up being -7dB (instead of 0dB).

Used the TI BiQuad Calculator tool and also the tool here

Don't know where things are going wrong. Any guidance is much appreciated. Thank you.

Used the formulas from p3 of this PDF

// Coefficients from the Calculator
N0  2.336268076
N1  -2.773236549
N2  1.01094234
D1  -0.680740702
D2  0.254714569

N0_Real 2.336268076
N1_Real -1.386618274
N2_Real 1.01094234

Scaling Factor  2.336268076 (Max of the above 3 values)
Range   32767 (Fixed for the Chip - 16 bits)

Digital Coefficients
N0  32767
N1  -19448
N2  14179
D1  11153
D2  -8346

Some information:

Sampling Frequency: 44100
Treble Shelf Center Frequency: 7kHz (to get required boost/attenuation at 10kHz)


texts: What is the phase response of SRRC/NRZ/RZ pulse shapping filters and does it matter?

texts: So I would like to know if filters usually used for pulse shaping (like SRRC) have a linear phase response, and if not, what this response looks like and does it matters in any circumstances?


texts: Can a digital image be considered as a tensor ? or have a tensor like structure?

texts:
I am wondering if one can assume that a digital image has a tensor structure. I mean for example when we talk about isotropy in image processing, we are talking about the intensity arranged in a tensor like structure; having coordinates and values as opposed  to a matrix which carries only information about values.

Why is this important ? because I would like to explain some processes in comparison with tensor functions. For example be able to explain isotropic operations

Thank you,


texts: Most efficient way to extract frequency components of audio signal with near-human accuracy in real time

texts:
I'm trying to work out how (if possible) to extract the frequency components of an arbitrary audio sample (typically music) in an FFT-like manner, but in my research on the FFT algorithm, I'm learning that it suffers some severe restrictions for this purpose.

There are 3 problems which the FFT is presenting:

Because the FFT bin-resolution is equivalent to your window size, to achieve a pretty reasonable accuracy (say 1 Hz), you need an unreasonably long window (say 1 second). This means you can't detect transients or newly introduced frequencies quickly. It also means that the problem can't be solved with a faster CPU and a higher sample-rate -- the restriction is intrinsically tied to time.
Humans perceive frequency logarithmically, but FFT bins are spaced linearly. Eg a difference of 20hz at the low end of our hearing is huge, whereas a difference of 20hz at the high end is imperceptible. So to get the accuracy we require at low frequencies, we have to calculate far more than we require at high frequencies.
Some of these problems can be resolved by interpolating between FFT bins. This may work for much musical audio, because frequencies will often be spaced quite far apart and so no more than 1 frequency will leak into a pair of bins. But this will not always be the case, particularly for inharmonic sounds like percussive instruments. So interpolation is really just guesswork.

From what I understand of the DFT/FFT algorithm, the outputs (bin amplitudes) are effectively the correlation of the sine/cosine at each bin's frequency. It strikes me that if the algorithm could be redesigned so that the bin frequencies are spaced non-linearly (i.e. we correlate a different set of sines/cosines), then we could achieve pyschoacoustically-equal resolution at all frequencies. Is this possible, or is it a pipe-dream based on my incomplete understanding of the maths involved?

I guess I could also solve the problem with brute-force, by correlating sines/cosines at every single frequency I'm interested in. I'm not too clued up on the maths here. Is this possible? What sort of efficiency? Would it solve my problem?

Is there a different way to achieve more accurate, real-time, frequency decomposition of a signal? CPU efficiency is a concern, but not the major concern -- I'm partly interested in whether it can theoretically be done at all. However, something that's feasible in realtime on a modern desktop machine would be ideal.


texts:
If we plot for example the relation between the probability of False Alarm (pfa) and the threshold. The threshold is in x-axix, while the $log _{10}(pfa)$ in y-axis.

If I get for example the following plot:

Which is better? the green one or the yellow one?
Thank you if you can give me some explanations.


texts: I'm writing a simple application in Python that sends and receives data over binary PSK using audio signal. I used cos() for 0 and -cos for 1. On the receiver side I figured I can demodulate by multiplying the signal by 2*cos so this will essentially shift the signal below and above x axis for corresponding bits. It works great. Now, I'm thinking about using low pass filter or maybe averaging the data in the time frames to get the output "square wave"? What would be the best/most efficient approach here? I'm might be a pro software engineer but tinkering with this as a hobbyist, thus probably very basic question :)


texts:
1)Suppose that we have low-pass filter with cf equal to 0.1. How frequency response would look like  if I multiply filter kernel by sine with frequency say 0.3?  How this is related to converting
low-pass filter algorithm into  bandpass filter?

2)What is the difference between multiplying filter kernel by sine and cosine?

Thanks in advance.

EDIT

I'm still don't quite understand the difference between multiplying  filter kernel by sine and cosine. It is said that with help of cosine it is possible to transform lpf to hpf. But I don't understand why?

Frequency spectrum for cosine is similar to sine (is it?) So what's the difference?

Thanks in advance again.


texts:
I'm trying to use the native functions of OpenCV HoughLines and HoughLinesP to detect the lines from an image. The Canny edge output of that image is the folllowing:

The problem is that the OpenCV functions I mentioned above can't detect any of those lines. I have fiddled with the parameter values, trying to get a better result, but nothing. I think the issue may be that there are too many contiguous lines in the image, but is just a guess.

Can anybody help me with this?
And if the algorithm is not suited, can you point out an alternative?

EDIT:

I mainly lowered the threshold value, but, unless I lowered it to much (case in which the returned image is fully covered with miss detected lines, i.e. the output image is a red rectangle), the output image is exactly the same as the input image, without any line detected. I also tried increase the angular resolution and use subpixel distance resolution, so instead of

HoughLines( dst, lines, 1, CV_PI/180, 100 );

I tried something like

HoughLines( dst, lines, 0.5, CV_PI/360, 20 );

Without any success.

EDIT II:

Here is the code:

void showLines(Mat src){
    Mat dst, cdst, color_dst;
    // Image only have 9 pixels vertically, so it needs to be resized
    resize(src, src, Size(), 2, 2, INTER_CUBIC);
    //Canny -> generates the image above
    Canny(src, dst, 50, 200, 5, false);

    vector<Vec2f> lines;
    // detect lines
    HoughLines( dst, lines, 0.5, CV_PI/360, 20 );

    // draw lines
    for( size_t i = 0; i < lines.size(); i++ ) {
        float rho = lines[i][0], theta = lines[i][1];
        Point pt1, pt2;
        double a = cos(theta), b = sin(theta);
        double x0 = a*rho, y0 = b*rho;
        pt1.x = cvRound(x0 + 1000*(-b));
        pt1.y = cvRound(y0 + 1000*(a));
        pt2.x = cvRound(x0 - 1000*(-b));
        pt2.y = cvRound(y0 - 1000*(a));
        line( cdst, pt1, pt2, Scalar(0,0,255), 3, CV_AA);
    }

    imshow("source", src);
    imshow("detected lines", cdst);

    waitKey();
}

The source (src) image is the following:


texts: Can I compute Fourier series without looping through all frequencies?

texts:
I need to compute Fourier series of an audio stream.
But DFT/FFT is slow.

Are there any ways to compute Fourier series of a signal without using the Fourier transform to check if whether a frequency is present or not?


texts: When is the sum of two (parallel) minimum-phase filters also minimum-phase?

texts:
Say I have two minimum-phase filters:
$$\frac{A(z)}{B(z)} \: \text{ and } \: \frac{C(z)}{D(z)}$$
That is, the roots of $A(z)$, $B(z)$, $C(z)$ and $D(z)$ are all in the stable region.
If add them together:
$$\frac{A(z)}{B(z)} + \frac{C(z)}{D(z)} = \frac{A(z)D(z) + C(z)B(z)}{B(z) D(z)}$$
The result is causal, and stable (because the roots of $B(z) D(z)$ are also in the stable region).  But it might not be minimum-phase, because the numerator $A(z)D(z) + C(z)B(z)$ can (in general) have zeroes outside the stable region.
What conditions are required for this summed filter to also be minimum-phase?
EDIT: I found a sufficient but not necessary condition here, but that's not quite what I was hoping for


texts: Get histogram of edge lengths/orientations in an image using OpenCV

texts:
In the process of trying to estimate whether a ROI contains [Japanese] text, I need to get a rough count of edges, hopefully binned by length and/or orientation... I think this would give me strong-enough features to judge if the ROI is likely to contain text or not.

However, I can't seem to figure out an efficient (this eventually has to work on a video stream) way to get such a histogram (or even just a count of edges for certain length/orientation threshold): Hough transforms will tend to focus on certain shapes (lines etc.) while Canny and Sobel will detect any type, but their output is not particularly useful for counting...

I would appreciate any suggestion or hint at some algorithm I can use, to quickly extract a list of edges from an image...


texts: I'm looking for a DSP, around 500MHz, with an in-chip ADC (preferably 2 but 1 will suffice) that samples around 20MSPS. Does any such DSP exist or do I have to use external ADCs?


texts: What is the expression for intensity of Laser beam focussed at a point at an angle theta?

texts:
If the laser beam is focussed directly at a point then the simplified expression is:
$$I_1(i,j) = 8\frac{P}{\pi d^2} e^{-\frac{8}{d^2}\left(x^2(i)+y^2(j)\right)}$$

If the laser is focused at same point but directed at an angle theta then what will be the intensity expression?

I am asked to work on a project related to quadrant detector. In this I have to plot the current variations of the four photodiodes. For finding the current I need the intensity of laser beam reaching the four photodiodes. So in case I focus the laser beam directly at a point then intensity is same in all four quadrants and that is given by expression above. If I focus the laser beam at an angle theta then what will be the expression of intensity. I want to know the intensity expression, so that I can plot it in Matlab.

I tried to replace d in above expression with d/cos(theta) but that is not correct. Next, I replaced x(i)^2+y(j)^2 with (x(i)^2+y(j)^2)*cos(theta)^2 but that is also not correct.

Please help me out with the correct result.


texts:
I have short-wave digital receiver that send RP packets with IQ components over network depend of the RF frequency and passband. I need to work with the telephone passband (0.3 - 3.4 kHz). The task is:

Using ordinary notebook receive RTP packets in the real-time, extract IQ-data, build FFT, visualize it and then I need to playback this data on the ordinary soundcard.

Please, help me out. How can I do this quickly (in the short time)? Maybe prepared solution exists.


texts: Machine learning with only positive labels

texts:
Suppose I have a binary classification problem with 10 features and about 1000 samples. In the training set, most of my data is unlabeled (75%). The rest of the data is labeled but contains only positive labels.

In the test set, I have both negative and positive labels. How should I approach this classification problem?


texts:
I've been working with an FMCW Fast Chirp radar transceiver.  Mainly, I've been looking at the Doppler range map but really it is not suitable as there is only one point per range/ speed and I want to build a better picture of the surroundings.

I came across the following article.  The below picture struck me as excellent:

Does anyone know how I could generate a picture as shown on the right using radar?


texts:
I have done the computation of the autocorrelation of a sequence (vector), here called "frame", and I have plotted it. Here the code:

The resulting graph is:

I do not understand the position of the peak in this graph: I know that the peak of the autocorrelation should be in the origin, why is it there?


texts:
Sorry that I created a new question about OFDM, but I donät have enough reputation to comment on other questions.

I have read all the topics here and I have some doubts.

1) I found in the book that If a system has a low symbol rate then --> negligible ISI. Why?
I have thought about it. Maybe because of that the system, in this case, has a longer symbol.

2) A cyclic prefix redices data rate. How can it be improved?


texts:
I have a question concerning BPSK and the bandpass filter that seems to be sometimes placed before the receiver circuit (the receiver is a mixer and lowpass filter).

Basically, from what I've noticed, there is sometimes a bandpass filter before the receiver of BPSK (and other types of receivers too I think).  Is the only purpose of this bandpass filter to band limit the signal and basically just improve the carrier to noise ratio?  It doesn't actually need to be there does it?  It's optional right?

Thanks.


texts:
I have a 3D-printer z-axis noisy signal. What I am trying to achieve is to detect each printed layer according to the z-value (e.g. z = -53.5 -> 1st layer, z = -52.5 -> 2nd layer ...).
However, due to the small measurement inaccuracies (marked in red), I'm having a hard time differentiating layers from each other.
Ideal would be replacing all errors by its previous value (e.g. all errors above or below 1st layer = -53.5 -> set them back to the 1st layer value)
How can I approach this issue using DSP techniques?


texts: Are the FFT coefficients symmetric in image processing?

texts: On page 11 of Fundamentals of Image Processing by Ian T. Young, Jan J. Gerbrands, Lucas J. van Vliet (pdf) the results of the Fourier transform are shown (Figs 4a and 4b) and it appears to me (please correct me if I am wrong), that the transform may exhibit redundant data in four quadrants. If this is so, would it not be possible (especially for highly symmetrical images) to take just the top left quadrant, flip it left and then take the top half of this result and flip it downwards. Would this not preserve much of the information when a reverse FFT is performed? If so, it would allow as 4x reduction of transmitted information before a compression algorithm such as JPEG algorithm was applied.


texts: Quantifying the frequency shift of a delay with LP filtered feedback

texts:
I am working at a Karplus-Strong sound synthesis technique where an excitation impulse is subjected to a delay with feedback and a one pole LP filter in the feedback line. The process is simply
y[0]=(x+ky)[-n]
Where x is the input signal,  y is the output signal,  n is the delay length in samples in the circular delay buffer and k is the feedback factor from 0 to 1. If no filtering is applied at the feedback line,  the resulting frequency will be

sample rate / n.

If I apply a simple one pole LP filter in the feedback line,  which is very effective for simulation of plucked strings,  of form

y[0]=x+f*(y[-1]-x)

where f is the filtering factor from 0 to 1,  I get a frequency shift which is in some manner proportional to f (the higher f,  the higher the shift downwards in frequency).

My problem is,  I have no means to quantify this shift,  the underlying math is likely too complex for me. But I would like to discover the exact relationship linking f and n so I can compensate n in function of f in order to avoid detuning (n is not integer but a real number becsuse I am using interpolation to read my samples from the circular buffer).
The only thing I discovered is that such function D
m=D(n,f)
always passes by 0 for f=0 and by 1 for f=0.5 and likely goes to infinity for f=1 but it is neither a suitably scaled hyperbole nor exponential but something else..
Thanks in advance


texts: How do I remove "jitter" when detecting frequency of Sine Wave

texts:
Input: Sine wave of any frequency centred at 1.66V
Output: Square wave with a frequency equal to input
I have been working on implementing a common DSP problem -- Detect the frequency of a Sine Wave. I have used a handful of posts from stack exchange to get me started, and found This answer best suited for my problem. However, implementing the solution and measuring both the input signal (sine wave) and the output signal (square wave) with an oscilloscope, it appears I am getting some "Jitter" on my output signal.
Does anyone know what could be causing this? Or is this just... as accurate a detection I can get?

Here is the code I am using to detect the frequency of the sine wave:
volatile int period;             // equal to the number of samples taken between zero crossings
volatile float frequency;        // the calculated frequency of the sine wave
int numSamplesTaken = 0;         // used in frequency calculation formula

int currValue = 0;               // the current sampled value of sinewave input
int prevValue = 0;               // the previous sampled value of sinewave input

int zero_crossing = 32767;       // ADC range is 0v - 3.3v, so the midpoint of the sine wave should be 1.65v (ie. 65535 / 2 = 32767)
int threshold = 500;             // for handling hysteresis
int isPositive = false;          // whether the sine wave is rising or falling

// interupt occuring at 8000hz
void sampleSignal() {

  currValue = input.read_u16();   // convert analog voltage input (sine wave) to a 16 bit number

  if (currValue >= (zero_crossing + threshold) && prevValue < (zero_crossing + threshold) && isPositive) {

    output.write(0);               // write digital output pin LOW
    isPositive = false;

  } else if (currValue <= (zero_crossing - threshold) && prevValue > (zero_crossing - threshold) && !isPositive) {

    output.write(1);               // write digital output pin HIGH
    period = numSamplesTaken;      // how many samples have occurred between positive zero crossings
    frequency = 8000 / period;     // sample rate divided by period of input signal
    numSamplesTaken = 0;           // reset sample count
    isPositive = true;
  }

  prevValue = currValue;
  numSamplesTaken++;
}


texts:
Stability of the below function was asked.

$$
T(x[n]) = \sum_{k=n_0}^n x[n]
$$

I understand that $x[n]$ shouldn't exceed $M$ value. $x[n]\le M$. I mean, we are looking for if the input is bounded, what will the output be.

However, i didn't understand from where $|n-n_0|M$ came from in the below function.

$$
\left |T(x[n]) \right| \le \sum_{k=n_0}^n \left| x[n] \right| \le |n-n_0|M
$$


texts: Removing DC Offset (via removing mean) of input and output signals affects their coherence

texts:
I have input $x$ and output $y$ signals (fs=10 Hz, length N=10800). I have removed DC offset via the following equation:
$x_{normalized} = 100 \frac{(x-\bar x)} {\bar x }$

where $\bar x$ is the mean of x, so as to get a normalised percentage of change

The code for this:

function x_norm = mean_norm(x)
  m = nanmean(x);
  x_norm = 100*(x-m)/m;
end

I have also attempted other methods such as high-pass filtering and using the detrend(,'constant') function.

I'm taking 1800 sample long epochs (so 3 minutes) at a time and performing system identification to get an estimate of the IR via the inverse FFT of the transfer function $H$, which I can then use to convolve with $x$ to find predictions of $y$.

So $X_i$ and $Y_i$ can be said to be $x$ and $y$ data for the epoch $i$.

The coherence (calculated using PSDs obtained using Welch's method) between $X_i$ and $Y_i$ before and after removing DC offset is significantly different to the coherence between them before removing it, and I cannot work out why. There is also a change in the spectra at non zero frequencies, which I did not expect.

What could be causing this change in coherence and spectra? I was under the belief that removing DC offset should not affect it, but perhaps I am wrong.

Thanks!

As an update, I was originally using a Hanning window, and leakage from that was having a big effect on the data. I've now changed to a Rectangular window, which has reduced the differences but there is still some, so it may be that any differences are just due to leakage.


texts:
I have a very simple code for ploting frequency spectrum of a signal in MATLAB. My center frequency is at 5e6 but I see a much higher peak at zero frequency despite the fact that I remove DC component from my signal. Why this peak happens?

Fs = 40e6 ;
[freq,z_freq] = freq_spectrum(RF_ch,Fs) ;
figure;
plot(freq,abs(z_freq));xlabel('Frequency');title('Frequency Spectrum')

function [f,z] = freq_spectrum(x,Fs)
N = size(x,1) ; %size(x)
f = -Fs/2:Fs/(N-1):Fs/2;
x = detrend(x) ; % or %x = x - mean(x);
z = fftshift(fft(x))./numel(x) ;
end


texts: Should the correlation check should happen at the start of the receive chain?

texts:
Should the correlation check should happen at the start of the receive chain.
does it act like a start of burst data? But does phase and freq offset has any effect on this?
or should it happen at bit level after the data has been recovered as the freq and phase will be already recovered?


texts: Create distortion from basic linear (and non-linear if neccessary) DSP elements

texts:
I'm studying mechatronics and I'm intrested in DSP basics. My lecturer said that there are four basic linear DSP elements:

Adder (and other mathematical operations)
Amplification (shown on diagrams as triangle)
Delay ($z^{-1}$)
Junction

I know from electronics that I can create distortion using for example op-amps and clipping circuits - and I can do it practically. But is there a way to create distortion using four DSP elements theoretically?
I know distortion is non-linear - so if not this four elements, what kind of elements existing in theoretically DSP should i use?


texts:
This paper refers to the concept of "spectrum perumutation" on page 4. I'm having a difficult time understanding what it means to permuate a signal. This is mainly due to the fact I don't understand the "invertible mod n" means.

What exactly happens to a signal when we permutate it? What does it look like in the time domain and what does it look like in the frequency domain?


texts: signal uncorrelated with delta in the origin

texts:
I have a stochastic process completely uncorrelated. Why the autocorrelation function has a delta of dirac in the origin?
Which is the reason of that?

$R_{XX}({\tau})=A{\delta(\tau)}$

where $R_{XX}$ is the autocorrelation function of $X$ process and $A$ is a constant dependent by process spectrum


texts: How to plot magnitude and phase response by hand if I have the Transfer Function?

texts:
I have the transfer function of the system, which is:

$$H(z) = \frac{1-z^{-1}}{5(1+2z^{-1})}$$

How do I sketch the magnitude and phase response?

I'm sorry for the bad formatting, it's my first time post a question.

Thank you very much in advance for the help!


texts:
Back around 5 BC, where C stands for the programming language, engineers would photograph their oscilloscope's display to record their test results.

This appears to be such a photograph, taken from this Wikipedia article on Elliptic filters.

The photo is described as follows:

Log of the absolute value of the gain of an 8th order elliptic filter in complex frequency space. The white spots are poles and the black spots are zeroes. There are a total of 16 poles and 8 double zeroes. What appears to be a single pole and zero near the transition region is actually four poles and two double zeroes as shown in the expanded view below. In this image, black corresponds to a gain of 0.0001 or less and white corresponds to a gain of 10 or more.

I have seen plots like this many times and always wondered how they were created. In other words, what methods and equipment were used to do this? They even show poles in the RHP, which doesn't make much sense.


texts:
I have heard that it is possible to detect a low power (ie not visible on a spectrum analyzer) DSSS-BPSK (direct sequence spread spectrum) signal even when you don't know the spreading code by computing the power spectral density of the incoming I/Q sampled datastream squared (ie PSD(x^2)). The dominant peak is the center frequency of the signal.

Why does this algorithm work and what is it called?


texts: False negative probability of cross-correlation

texts:
Suppose $D$ is a string of length $L$ where its values are uniformly distributed in $[-a,a]$. Also, all values of $D$ are independent (i.i.d).

$X$ is a noisy version of $D$ in this way: $X=D+N$ where $N$ is a white Gaussian noise with mean and variance of $0$ and $\sigma_N^2$.

$P$ is a sub-string of $D$ where its length is $l$  $(l<L)$.

To detect $P$ in $X$, I use the cross-correlation of $X$ and $P$:
$$ C_n = \sum_{i=1}^l  X_{n+i-1}P_i \space\space , \space\space 1<n<L-l+1 $$
and a threshold of $T=K\sigma_{C_n}$ where $K$ is a constant value and $\sigma_{C_n}$ is the standard deviation of the correlation signal $C_n$.

A false negative occurs If no points of the correlation signal reaches the threshold.

Now the question is: how can I theoretically calculate the false negative probability? It tried to use Q-function, but I couldn't achieve the answer.


texts:
What will be the new Fourier series coefficients when we shift and scale a periodic signal?
Scaling alone will only affect fundamental frequency. But how to calculate new coefficients of shifted and scaled version.
I tried searching, but couldn't find an answer where both properties are used.
Please help.


texts: What is the "Energy" equivalent of the Magnitude Spectrum

texts:
I have a semantic issue regarding the discrete spectrum $X[k]$ of a time domain signal with limited length $x[n]$, so we're talking about the FFT. Now, when I what to compare the bandwidth and amplitude of sinusoidal (a peak) in $|X_1[k]|$ vs. another $|X_2[k]|$, am I allowed to talk about the sinusoidal having more or less energy or having the energy spread over a wider frequency range?

As far as I understand, the energy spectrum is defined as $|X[k]|^2$ and when I would directly use values, I would need to take these in order to talk about energy. Though In practice I'm working with the decibel magnitude spectrum $20 \log_{10}|X[k]|$ anyway. My argument would be that even though the graph I'm describing doesn't directly display the energy, it displays a value directly proportional to energy, especially in the  decibel case were as far as I know the difference is simply a factor of $2$, as in

$$
20 \log_{10}|X[k]| = 10 \log_{10}|X[k]|^2
$$

If this wouldn't be technically correct, is there an equivalent to "energy" in the magnitude spectrum? I've never heard someone explain that the magnitude of a sinusoidal is spread over a larger frequency band.


texts: Is there any new scientific capability to build non-causal filters in real world?

texts:
According to this post:

  In discrete-time systems, causality is a requirement only when
  processing (filtering) signals in real time; i.e. when the
  index nn corresponds to a physical time n×Tsn×Ts. In this case, a
  non-causal system is impossible to implement since calculating the
  current output would require inputs from a future time.

  In practice, however, many systems are not real time.

  Image processing, for example, has indexes for horizontal and vertical
  coordinates. So it's OK to apply a non causal filter to an image
  (usual examples being 0-centered FIR filters for blurring or
  sharpening, or IIR filters applied in both directions to obtain 0
  phase distorsion).

  In time sequences, many times the signal is stored and then processed
  in non real time, so you may apply non-causal filters there also.

I like to know is there any new scientific theory which could yield to making none causal FIR Filters in the real world. Like quantum theory in physics or other scientific fields?

Thanks for your attention.


texts: Problem Involving Aliasing and Critical Fusion Frequency Concepts

texts:
I have been trying to study a course from MIT on Biomedical Signal Processing and I thought that the best way to learn would be to solve some problems. Although they might be out of my depth, I am not looking for so much as quick answers as ways to understand the problem and how to remedy my knowledge gaps in this area. Specifically, I would like to solve this problem:

If you have ever seen a Western, you must have noticed that, as the stagecoach changes
speed, the wheel appears to reverse direction even though the coach moves forward. This
is an example of aliasing. The sampling mechanism is the movie camera capturing images
at a frequency of 24 frames/second. The eye acts as a lowpass filter, in the sense that
successive images reaching the eye at a rate exceeding 10 images/second appear as if they
are superimposed. This is called the critical fusion frequency.
Answer the following questions assuming that the wagon wheel is 1.6 m in diameter and has
12 spokes.

(a) Express the true angular velocity of the wheel (in radians/second) as a function of the
coach velocity.

(b) At what set of coach velocities does the wheel appear to reverse direction?

(c) Determine an expression for the apparent angular velocity of the wheel as a function of
the coach velocity.

(d) How many spokes does the wheel appear to have the first time it appears to reverse
direction? The second time? Explain the difference between these two cases.

I think I understand aliasing, whereby because of the low frequency of the sampling-rate it causes the output signal to be of lower frequency than it actually is, and in order to solve this problem a low-pass filter is required. However, I am struggling to begin to tackle this problem as I am stuck at the first question of expressing the true angular velocity of the wheel as a function of coach velocity.

Therefore, I would be much obliged if someone could provide me with some hints or some sort of a guide as to how I should tackle these problems.

Edit: I think I overthought (a) and it looked to be pretty simple. So that would mean angular velocity of the wheel is simply the velocity of the coach divided by the radius of the wheel which in this case is 0.8 m.


texts: Sampling interval $T$ as a multiplier in digital processing of a continuous-time signals

texts:
[from: Discrete-Time Signal Processing, Oppenheim and Schafer, p.224]
Q:  Why do we have $T$ as multiplier in $TY_a(j\Omega)$ in Eq.155?


texts: Symbol Synchronization for GFSK/FSK signals in GNU Radio

texts:
I am trying to demodulate GFSK signal, and extract the packets using GNURadio.

TX information:

It is basically BLE transmission but customized. BW: 2MHz, data rate: 1MSymbols/sec, 1 bit per symbol. Each packet is 160 bits (which equals to 181us duration). And the packets are repeated every 100ms, hence approx. 10 packets per sec

RX information:

I am using a BladeRF with 4Msps and storing it to a .complex and read back in GNU Radio.
The samples are given through an LPF performing decimation by 2 which reduces the sample rate to 2Msps. Then I use a quadrature demodulation block following which I use a Symbol Sync block. After this I use a custom OOT module, to count the packets.

So, everything is perfect until the demodulation, after demodulation I can clearly see the data stream. Only thing is that we can't make a hard decision here  (we have 2 samples per symbols). In order to fix that we pass it through the Symbol Sync block. Then we do packet detection

Problem:

My main issue is not every packed is synced/recovered by this block. As most of the symbol timing recovery is a kind of estimation/control loop algorithm. Setting the correct parameters is very important. I played a lot with the parameters. And whatever I have set now works best for me.
Even in this best-case scenario, I am losing some packet. The capture file I have is for 32.03 seconds, which approx. contains 319 packets, but I can only recovery 311 packets.

I have been experimenting with this for nearly a month, with different packet rates, like 10, 25, 50, 75, 100 packets per second. All the test capture is done in a very low noise environment, and RX close to TX. So, any channel effects can be ignored. I even tried inserting a matched filter before demodulation, use a polyphase clock recovery, etc. but still, packets are lost.

And every time packet loss rate is always 5% to 7%

Questions:

1) why exactly this happens? is it because that the limit of symbol recovery can do or Am I doing something wrong?

2) is there any other alternative for this problem, to recovery 100% or at least with less than 4% loss? I have seen a technique which reverse engineer like M.ossmann's clock recovery use (https://github.com/mossmann/clock-recovery). But, I want to do it in real-time.

Any other suggestions to modify the algorithm or parameter or even approach are highly welcomed and appreciated....!!! I stuck with this problem for soo long and it's getting on my nerves.

Please, I am calling out experts in the field for help and looking for discussions.

Thanks in advance

P.S: I am attaching my screenshot of the GNU Radio flow graph. And also a plot showing how some packets are distorted after passing through the symbol sync block. The red is the bits of a packet after the symbol sync block and blue is before that. You can see that it's been lost many bits. Hence, the packet detection block is unable to detect it.

UPDATE:

Based on the suggestion from the comments, apparently the loop bandwidth and damping factor value are not very optimal. As my preamble is 8 symbols the loop should settle to linear operation within this time. This corresponds to the loop bandwidth of ~40KHz (normalized loop BW 20m). After these changes also I was still getting a 5% packet loss. In order to analyze a bit deeper, I created a bit more detailed graph with the debug signals from the Symbol Sync block. I have superimposed the recovered packet with error, instantaneous period, average period debug signals. I am attaching the graph for both cases, where the packet is recovered/detected good and the failure case. The first image corresponds to a good packet, the second image is for a bad packet.


texts: Are there libraries for extraction of sound wave features?

texts: I am looking help for my project for which I need C++ (or any other language) libraries useful for extraction of sound wave features like frequency, loudness, pitch and orientation. Using these features I want to make attention in a humanoid robot simulator .


texts: Energy definition for Autocorrelation lag 0 and lag 1 for complex signals

texts:
I am studying the role of an auto-correlation matrix for random signals and the difference of energy between a lag 0 and lag 1 matrix.

Consider a complex input signal $x(k)=[x1,x2]^T$ and $x(k-1)=[x0,x1]^T$, as column vectors with auto-correlation matrix of lag 0, $R(0)$, with diagonal $r(0)$ and auto-correlation matrix of lag 1 $R(1)$, with diagonal $r(1)$.

where $R(0)=E[x(k)x^H(k)]$ and $R(1)=E[x(k)x^H(k-1)]$

what is the relationship between the energy resulting from $R(0)$ and $R(1)$? can we thus say that eigenvalue max of $R(0)$ is greater then that of $R(1)$
and the difference in diagonal entries $r(0)$ and $r(1)$?

given that $r(0)$ is always positive real and $r(1)$ can be complex.

What can we comment of the real and imaginary part of the $r(1)$ entry?
and where does the energy resides i.e. real portion or imaginary?
Additionally to a wide sense stationary signal $x(k)$ can we assume that $r(1)$ is approximately equal to $r(0)$?

Thus if we have $B = R(0) + R(1)$, we know that $R(0)$ is positive definite, symmetric and hermitian however $R(1)$ is not. Can we say that eigenvalue max of $B$ is eigenvalue max of $R(0)$ + the norm of that of $R(1)$ assuming its complex?

My apologies for the successive questions however I am trying to understand the differences between the auto-correlation matrix at lag 0 and lag 1 in terms of eigenvalues and energy because I have to work with a matrix $B = R(0) + R(1)$.


texts: Identifying an Event based on a Jump in Rolling Variance in Python

texts:
I am trying to figure out a way, in Python, how to identify when an event occurs, based on a jump in the rolling standard deviation.

As shown in the plot below, around the 12000th sample, an event occurs.  In my Python script, I am currently using a threshold of 0.00051 to signify when that event occurs.  However, sometimes the event occurs at 0.0005 and other times the event occurs at 0.000495.

My question is - how can I algorithmically, in Python, detect this jump in the rolling standard deviation to create the event altert?  Because if I have a threshold set too low, I don't want it to be triggered too early.  And if I run another test, and the threshold is too high, then I don't want the event to not be triggered at all.

Any advice is greatly appreciated!

def animate(i):
    data = pd.read_csv("C:\\Users\\Desktop\\data.txt", sep="\[|\]\[|\]",engine = 'python', header = None)
    data = data.iloc[0, ::4]
    data = data.astype(str).apply(lambda x: x.split(',')[-1]).astype(float)
    data.pop(0)
    xar = range(len(data))
    yar = pd.DataFrame(data)
    # Starting from sample 1050 to get rid of any initial noise
    yar = yar[1050:12500]
    xar = xar[1050:12500]
    std = yar.rolling(window=2500).std()

    if (np.any(std>.00051)):
        choices = ["Confirm Event"]
        reply = easygui.buttonbox("Event Alert!, image, choices)
        if reply == "Confirm Event":
            sys.exit(0)

    ax1.clear()
    ax1.plot(xar,std)
    ax1.set_title('Rolling Standard Deviation')

fig, (ax1) = plt.subplots(1, sharex = True)
ani = animation.FuncAnimation(fig, animate, interval=.01)
plt.show()


texts: What is the difference between Spectrogram and Scalogram?

texts:
Could someone please explain the difference between the two?

Spectrogram:

  A visual representation of the spectrum of a sound changing through time.

Scalogram:

  (signal processing) A visual representation of a wavelet transform, having axes for time, scale, and coefficient


texts: Why do we use conjugate transpose IFFT in OFDM

texts:
I have a question regarding the OFDM and using the conjugate transpose of vector (matrix).

As you know, if we have a block of data,let's call it $X_n$, consisting of $N$ data symbols, we should apply for IDFT for that block, and let's say that we get $Y_n$ which is $Y_n = \text{IDFT}(X_n)$.

After that step, the conjugate transpose for the IDFT must be taken.  So what's the benefit from taking that? If you have documents which explains that in simple way, it will be appreciated.


texts: Derivation of fixed-point $\tt atan2$ with self-normalization

texts:
I'm trying to understand the maths behind this Fast fixed-point $\tt atan2$ calculation with self-normalization. In particular, equation $(2)$ for theta1 appears to provide a first-degree expansion of some series. Similarly, equation $(2a)$ under "QUICK NOTE" seems to expand this series.

What is the series that's being expanded?
How exactly is equation $(2)$ derived?


texts: Finding the Sections of Signals that are the Most Similar

texts: I have two discrete signals, $A$ and $B$ (both with $L$ values). I want to look for the section (of a certain length $R$) on each signal ($[A_1:A_1+R]$ and $[B_1:B_1+R]$) where they are the most similar (using some similarity metric, say Manhattan distance) (technically, I just need to know how similar the most similar sections of the two signals). Right now, I have a double loop that compares each section of signal $A$ to each section of signal $B$, and choose the $A_1$ and $B_1$ that lead to the most similar sections. However, that gets really slow as I scale up. Is there some other, more efficient way of doing this?


texts:
So I have been learning about how to test systems for causality and stability but I am confused about the implications on their unit circle representation.

Would it be safe to say that a causal and stable allpass filter does not have a zero inside the unit circle?


texts: Superdirectivity Beamformer for Circular microphone array

texts:
The matrix below represents the inverse of correlation matrix assuming diffuse noise field. I am unable to figure out the values of the elements of this matrix.

how to calculate $c_{M}(\omega)$ ?


texts: Highest frequency of collected force signal

texts:
I have collected force at the footrest during paddling at two occasions. During one occasion I accidently collected the force signals at 150 Hz instead of 1500 Hz. During the other time the data was collected at 1500 Hz. I now want to calculate the highest frequency of the signal collected at 1500 Hz to be able to know if the data collected at 150 Hz is within the Nyquist criterion.
How do I calculate the highest frequency in my signal?


texts: Low Level Audio Engine Audio Data Reformatting - Sample Rate and Bit Depth Conversion

texts:
First off, there seems to be a sparse number of low level audio programming tutorials out there so if anyone knows of some good resources I would be very appreciative. Now on to the topic at hand!
I am using WASAPI for the audio engine that I am trying to develop and I am having difficulties reformatting read audio file data into a format that the audio device is prepared to receive. Specifically I am wondering how to on the fly convert the sampling rate and bit depth of a read audio file into the requested format of the device. I saw an example where linear interpolation was used but read elsewhere that such a filter introduces too many artifacts such as aliasing, so again I'm wondering what practical algorithms would be for reformatting audio data in a real time application.


texts: What is effect of fractional sample delay effect on FIR Filters?

texts:
In the DSP lecture slides we have this slide:

And

So why this part of FIR Type 2 and Type 4 filter impulse response equation bolded in this slide?

Thanks for your attention.


texts: How to make sure that a sinusoid wave with different subsequent frequencies matches up?

texts:
A very basic formula to transform a sequence of notes into a signal would be:

$$
x(t) = A\cos( 2  \pi  (2^{current\_note}/12) t)
$$
where $current\_note$ is $notes[\lfloor t / period\_per\_note\rfloor]$ and $a[k]$ is the indexing operation into a string; in this case into a sequence of numerical values representing notes on a chromatic scale.

Listen to an example

However, when making a signal like this, whenever the frequency changes (when moving from one note to the next), the function is discontinuous, often starting at a completely different amplitude because in general, the periods of the two sinusoids with different frequencies end up being very different (not being a multiple of the same 'base frequency')-

This can also clearly be heard in the example, as 'clicks'.

Now the question: How can we make sure that the sinusoid with the 'new' frequency starts at the same value of x as the one with the old frequency?

I am guessing that this somehow involves keeping track of the current 'phase', so that we can time-shift the new frequency to make sure their amplitudes match up. But I have no idea how to approach this, neither from a mathematical nor from a programmatic perspective.


texts: Help with Time shifting a discrete signal

texts:
I have the following problem statement which I have no idea how to approach,

I have a IR depth sensor, that produces a distribution for it's depth value. So it kind of looks something like this:

As you can see, it shows that there are two possibilities, one at 10 and one at 20. Taking the expectation of this will give me the real depth.

What i want to do is apply some signal processing to model the motion of this camera in the Z axis. Basically i want it to shift say 5m and spread out to represent the increasing uncertainity. Hopefully to look something like this

I guess for the spreading out part, applying a 1D gaussian convolution would work, but i am not sure about the timeshifting part. Can i convolve the signal below with another signal to actual timeshift it somehow?


texts:
Say we have a sampled signal for which the inputs can be out of range. For instance, say a sampled sensor has a range of 0 to 1, and is presented with a time-varying stimulus where a peak occurs with max amplitude 1.2. In the data we would see a saturated or clipped peak.

What are some of the approaches to making estimates about the peak value?

How good of an estimate can be made by interpolating the value of the peak from the slope on either side of the clipping? This becomes difficult because if analysis relies on the dataset being at a specific sample rate, interpolation will result in upsampling, unless interpolation is done and the new signal is downsampled with an offset so that the peak becomes a data point. I'm sure this isn't an ideal approach.

What about calculations in the frequency domain?

Has anyone ever been presented with this problem?


texts: References and resources on hyperspectral (multispectral) image processing

texts:
Due to my work, I just started studying hyperspectral or multispectral image processing.

I think I have enough programming background and studied mathematics for my degree - it was long time ago though. Since early this year, I have been exposed to image processing in my development work for manufacturing industry. And I got involved in a research project dealing with the subject above.

I wonder if anyone can recommend good references or resources for the subject. Especially, I would like to be familiar with that kind of image data. Any advice would be appreciated.


texts: What is the most reliable shape representation for the purpose of detecting humans shapes?

texts:
The requirement of the algorithm is to detect human shapes. While the standard HoG or Harr cascade can be used, I am exploring whether the human shape - which primarily includes face+neck+upper-torso shape as some kind of template and identify human shape.

It will facilitate human detection in extremely low light condition where practically all i one can see is a "shadow like" shape of the human.

I am looking at images that looks .

Any suggestions?

PS: I do not have computational power (GPU) to run deep learning. So looking for traditional vision based algorithm.


texts: Blended use of DoH and LoG for blob detection, is there any paper I can cite for this?

texts:
tl;dr: For blob detection, I determine locally maximal DoH responses in image regions which correspond to negative LoG responses. Is there publication I can cite for such blended use?

Two of the most prominent approaches to blob detection in image processing are Laplacian of Gaussians (LoG) and Determinant of Hessians (DoH). Both are based on analysis of local second-order image statistics. One could say, they consider a second-order Taylor expansion of the image intensities $g\left(x\right)$. Lets say, the Matrix $H_x = \nabla^2 g\left(x\right)$ is the corresponding Hessian matrix for an image point $x$. Then, LoG corresponds to
$$
\Delta g\left(x\right) = \frac{\partial^2 g}{\partial x_1^2} + \frac{\partial^2 g}{\partial x_2^2} = \operatorname{tr} H_x
$$
where $\operatorname{tr} H_x$ is the trace of the Hessian (computed using Gaussian derivative operators). So obviously, LoG corresponds to the sum of the eigenvalues of the Hessian.
For DoH, lets also use Gaussian derivative operators for the computation of the Hessian. The determinant of a matrix corresponds to the product of its eigenvalues. Thus, the only difference of LoG and DoH is that LoG considers the sum of the eigenvalues, whereas DoH considers the product.
Lets say the two eigenvalues are $\lambda_1, \lambda_2$, ordered such that $\left|\lambda_1\right| > \left|\lambda_2\right|$. Using LoG to detect bright blobs, we seek local maxima of $-\lambda_1 -\lambda_2$. This is prone to saddle points, which occur when $\left|\lambda_1\right| \gg \left|\lambda_2\right|$ and $\lambda_1 \cdot \lambda_2 < 0$. DoH on the other hand aims to determine local maxima of $\lambda_1 \cdot \lambda_2$ and is thus not prone to saddle points. The disadvantage of DoH is that, it is impossible to tell whether a high response of $\lambda_1 \cdot \lambda_2$ corresponds to a pair of positive (dark blob) or a pair of negative eigenvalues (bright blob).
So it is quite obvious to combine the best of the two: First, use LoG to determine the set of image points $M = \left\{x : \operatorname{tr} H_x < 0 \right\}$ which cannot correspond to dark blobs. Second, use DoH to determine the blobs from the remaining image pixels, i.e. determine the local maxima of $\operatorname{det} H_x$ subject to the constraint that $x \in M$.
I expect this is more robust than discarding false-positive DoH detections in post-processing, since the false-positive DoH detections might prevent true-positive detections due to non-maximum suppression.
Question: Does this blended use of LoG and DoH have a name? Is there any publication I could cite when using it?


texts: Could multiple echo reflections be combined to improve signal to noise ratio?

texts:
I am working on a system processing ultrasonic echo reflection. The transmitter creates a short ultrasonic pulse and the receiver analyses the reflection.

If I send two pulses (back to back or with some small defined delay) I receive two superimposed reflections.
I can process these superimposed reflection with the IIR filter described in this post:

Trivial echo separation

to restore the single reflection.

My question is, could I get a signal with better signal to noise ratio by using two or more transmit pulses compared to the signal from a single pulse?


texts: AWGN interpretation of bit error rate in time domain

texts:
I'll consider a 20 MHz bandwidth between 1.99-2.01 GHz. The carrier has been modulated by an OFDM signal (I/Q modulation) and is passed through a narrowband filter at the receiver.

The noise that is actually received after the input filter is narrowband. So, it is not white anymore. The noise looks like a sinewave of 2 GHz with varying amplitude and phase, which has a Rayleigh distribution. In time domain, something like this:

This is something like the aliasing measured when undersampling a signal, but it has a Rayleigh distribution.

This noise is overimposed on the OFDM signal, upconverted on the carrier at 2 GHz, which in the time domain looks like this:

BER is expressed in the simplest case for an AWGN channel. I've seen that the procedure is to add AWGN of a certain SNR to the signal. The representation of this in time domain is not how the signal measured after the filter would appear, in my understanding.

An example of AWGN addition is this (https://www.mathworks.com/help/comm/ug/analog-passband-modulation.html)

Looking at the picture above, considering that the signal is sampled at the Nyquist rate, we can observe that there is an uncertainty in the values that are sampled. Could the same uncertainty appear with a bandlimited gaussian noise of the same power? Is that the baseband received signal that is represented above?
Are there cases when the white gaussian noise is not equivalent to the real bandlimited gaussian noise?

BER in AWGN appears to be the most common metric of noise performance, so there must be an explanation for this method. What is it?


texts: extracting audio and video from video file

texts:
I tried the following code in MATLAB R2013b

clc;
clear all;
VideoFreader=vision.VideoFileReader('video7.mp4','AudioOutputPort',true,'AudioOutputDataType','int16');
i=1;
while ~isDone(VideoFreader)
    [I,Audio]=step(VideoFreader);
    a(i)=Audio;
    i=i+1;
end

There are no errors, but warning as follows:

Warning: The AudioOutputPort property is not relevant in this
configuration of the System object.
> In C:\Program Files\MATLAB\R2013b\toolbox\matlab\system\+matlab\+system\SystemProp.p>SystemProp.parseInputs at 802
  In C:\Program Files\MATLAB\R2013b\toolbox\matlab\system\+matlab\+system\SystemProp.p>SystemProp.setProperties at 294
  In C:\Program Files\MATLAB\R2013b\toolbox\vision\vision\+vision\VideoFileReader.p>VideoFileReader.VideoFileReader at 133
Warning: The AudioOutputDataType property is not relevant in this
configuration of the System object.
> In C:\Program Files\MATLAB\R2013b\toolbox\matlab\system\+matlab\+system\SystemProp.p>SystemProp.parseInputs at 802
  In C:\Program Files\MATLAB\R2013b\toolbox\matlab\system\+matlab\+system\SystemProp.p>SystemProp.setProperties at 294
  In C:\Program Files\MATLAB\R2013b\toolbox\vision\vision\+vision\VideoFileReader.p>VideoFileReader.VideoFileReader at 133

But the Audio samples are of type logical and values are all 0

What is wrong in this? How do extract audio samples in double or int16 or uint8 type?

Please help me


texts: I have searched this question but couldn't find the answer in this network. I know this is very confusing question for DSP beginners. Both DFT and Z-transform work for Discrete signal. I have read that "Z-transform is the general case of DFT, when we consider unit circle then, Z-transform becomes Discrete Fourier Transform (DFT)". What does this mean? Ok, I can understand the mathematical verification but what is the physical meaning of this and how this affect the analysis in DSP?


texts: How Can I Add noise to a Frequency Domain Input Signal?

texts:
I am new with Matlab and doing some work on Amplifier and Bipolar Shappers Amplifiers. Adding Noise to a Time Domain Signal is quite easy and I have done it multiple time as:

T=10/wu;
N=100;
time=1:1:N;
time=T*time/N;
noisevec =0.55e-3*sqrt(0.707)*randn(size(time));
y22t=y2t+noisevec;

But Now I want to add a particular bandwidth Limited Noise to INPUT SIGNAL in S-DOMAIN which is :

x1s=(exp(-s*to1))*1/s-(exp(-s*to2))*(1/s);

Can somebody help me with it? It may a simple or dumb question, sorry for that, BUT i just dont know how to add noise to a signal in S domain
I tried to take Laplace of "noisevec =0.55e-3*sqrt(0.707)*randn(size(time));" but that does not help.


texts:
I'm doing research on how we can do matching between 2 audio files using Android to see if they have the same content or not. After doing lots of research, I have concluded that I must do 2 main things:

FFT
Correlation

Am I on the right track?


texts: Calculating an impulse response in a stable environment

texts:
I would like to estimate the impulse response of a simple environment (room, two speakers, one microphone). I am aware of the issues/limitations of working with "near" LTI systems. That said, I have a rather basic question regarding the mechanics of performing cross-correlation using audio which has been recorded in my environment.

A: my original audio (swept sine) - 1 second in length, 16 bit, 44.1K

B: the result of recording A in my environment

Again, my goal is to calculate the environment's impulse response, h[n], by calculating the cross-correlation of A and B.

Here is my question: I don't know if I should be doing anything to B (trimming on either end, shifting towards the left, etc.) before attempting to estimate the impulse response. The goal will be to use the IR to model what arbitrary audio will "look" like when played in my environment.

Additionally, I'm not sure how long (in # of samples) my IR should be. I've read that using a shorter IR can speed up processing without too much loss in estimating y[n] (the output of the system).

In the final version of this experiment the recording of B will be automated thus I'll need to know how long the recording needs to last to fully capture A's signal.

Thank you in advance.


texts: What's a good physical description of the following matlab code?

texts:
I've inherited some matlab code and I can't really understand what the resultant signal really means. The input is audio and the code first performs an STFT/spectrogram (sliding window FFT of size 1024 with 50% overlap, hanning window applied). The output is then filtered to the range (using Matlab notation) to the range (2:N/2). In other words, the first bin is dropped and the negative half frequencies too. Just to be clear, for each 1024 window of audio data (audio_chunk below) the code does:

y = abs(fft(hann_win.*audio_chunk));
filt_y = y(2:length(y)/2);

The value filt_y is then summed and that becomes the output signal:

output = sum(filt_y);

So, for every window's FFT we get a single output value. The output sampling frequency is audio_sampling_frequency/512 because we slide the 1024 sample window along by 512 samples at each step.

My question is: what would you call that output signal? What sort of physical significance does it have, or what use is it for analysis?


texts:
I record a sample of 2000 data points, with a sampling rate of 40khz, so freq. resolution is 20hz. I multiplied the chunk of data with a Hanning window prior to doing the FFT.

I recorded a 1khz test tone, and the result of FFT showed a peak at 1khz, as it should. But the magnitudes at 980hz and 1020hz is also significantly higher than the rest, although not as high as the magnitude at 1khz (about 50%).

Is this really bad? Did I apply the Hanning window wrong or is this unavoidable? Will using finer frequency resolution help?


texts: Computing Energy of sparsly sampled spherical function

texts:
I'm trying to understand the following paper: http://www.cse.cuhk.edu.hk/~ttwong/papers/shfit/shfit.pdf
On page 8, the energy of the function Ec is mentioned. How do I estimate the energy of a spherical function based on a handful of samples?
Wikipedia says on the energy of a discrete function:
$E_{s} \ \ = \ \ \langle x(n), x(n)\rangle \ \  =  \sum_{n=-\infty}^{\infty}{|x(n)|^2}$
but isn't it necessary to average over the area of the domain somewhere?
if I have sparse samples on the spheres, do I approximate the energy simply by summing the squares of my samples, or do I sum the squares of the samples, divide by the number of samples and multiply by the area of the sphere?


texts:
What happens if the controller of a closed-loop system is placed in the feedback path before the summation point, instead of between the summation point and the plant to control? Are we allowed to do so? If so, what are the main differences?
Thanks!


texts: Importance of Discrete Cosine Transform over High Pass filter?

texts:
I have recently studied few work where the low and high frequency components of an image is accessed using the discrete cosine transform (one such example is presented here).

Low and High frequency components of an image can be easily accessed by applying a high pass filter in transformed domain.

My question is that what extra importance does DCT have over the simple high pass filter in accessing high and low frequency components of an image. In other words instead of applying DCT can't we directly access these components using High Pass filter?


texts: Initial value Problem with difference equation

texts:
How to solve the initial value problem with

x(n + 1) + x(n) = n^2;
x(0) = 1

I do know how to solve for initial value theorem in differential equation, but not in difference equation.

Please help.Thanks.


texts:
Is there a canonical or analytic expression for the probability distribution for the circularly-symmetric complex random variable $Z$:
$$
Z = e^{j\theta},
$$
where $\theta \sim \mathcal U(0, 2\pi)$?

Side notes:

It is known that the real and imaginary parts, i.e.:
$$
\Re(Z) = \cos \theta \\
\Im(Z) = \sin \theta
$$
have marginal densities given by:
$$
f_{\Re(Z)}(z) = f_{\Im(Z)}(z) = \frac{1}{\pi\sqrt{1-z^2}}, \quad -1 < z < 1,
$$
but because they are not independent, computing their joint PDF is nontrivial.

EDIT: $Z$ is different from a complex normal in that here, the amplitude $|Z|$ is deterministic and identically 1, whereas if $Z$ were complex normal, $|Z|$ would be Rayleigh distributed.


texts:
My text book (Communication Systems analysis and design P.334, 335) have a FM signal is

  $s_{FM}(t) = A cos[2\pi[f_c+ks(t)]t+\phi]$ where k is a constant and $\phi$ is the arbitrary phase angles

then the textbook suggest

  while in FM, instantaneous frequency minus carrier frequency is a linear function of s(t)

So, I try to prove the $s_{FM}(t)$ is a FM signal

The instantaneous frequency is $f(t)= \frac {d\psi}{dt}=2\pi f_c+2 \pi s'(t) t +2 \pi s(t)$

and the instantaneous carrier ferquency is $f_{carrier} (t)=\frac{d\psi_{carrier}}{dt}=\frac {d}{dt} [2 \pi f_c t + \theta]=2 \pi f_c$. Assume the carrier signal is $Asin(2 \pi f_c t + \theta)$

So, $f(t)-f_{carrier}(t)=2 \pi s'(t) t +2 \pi s(t)$.

And my problem is why the text book suggest that signal is a FM signal? Clearly, I cannot prove the signal is linear.


texts:
In Kaldi, "povey" is a window  made to be similar to Hamming but to go to zero at the
edges, it's:

$$g(n) = \left(\frac{1-\cos\left(\frac{n}{N2\pi}\right)}{2}\right)^{0.85}$$

Please, can anyone mathematically explain why we are considering raising to the power of 0.85? Why not 0.8 or 0.9?


texts: Is there a method to compute an average amplitude for 2D sinusoids using the F.C.'s?

texts:
I am working on FFT analysis AFM like the one shown below. I am currently at a point where I understand fairly well what an FFT does, and how to think about the information it gives. I currently have an $M\times N$ matrix of Fourier coefficients (FC) for a 2D-image like the one below. Choosing to look at the real and imaginary part of the FC's gives the amplitudes for sine and cosine terms in the Fourier expansion. I use the notation $X[k,l]$ to index the FC's. For some of my images the real and imaginary parts of the FFT look something like what is shown below.

Sample Image:

Real part of FC's:

Imaginary part of FC's:

Note that the above FC plots are not for the image shown above. I am looking now to pull out the dominant functional forms. In general with these experiments it is expected that there will be 2 two-directions along which there are dominant sinusoids.

Are there standard ways of computing an average amplitude for 2D-FFT's? I was thinking that there is probably a way to perform a weighted average using the magnitudes of the Fourier coefficients. References and suggestions appreciated.

Additional:

I am thinking something like this in order to compute the average amplitude:

$\langle X \rangle=\sum_{k=0}^{M-1}\sum_{l=0}^{M-1}\frac{X_R^2}{X_R^2+X_I^2}X_R(k,l)-\frac{X_I^2}{X_R^2+X_I^2}X_I(k,l)$

The terms in front of $X_R$ and $X_I$ are my crude way to weight the importance of each sinusoid to the total signal. I have no justification for my choice, which is something I would like clarified. Perhaps there is an accepted way to weight the amplitudes. Clearly, my average is bad because of when $X_R=X_I=0$ I get a 0/0.

Cheers.


texts:
Can anyone explain what are Autoregressive Coefficients? What is their meaning that is.

Consider a method:

public double[] calculateARCoefficients(double[] inputseries, int order)

When this method is called on 256 values, lets say these values represent some kind of signal, I will get an array of n numbers, where n=order. What do these numbers represent?

I mean for example variance tells me how spread out the values of the signal are. The higher variance value is, the more spread out values are. The smaller variance is, the more unified signal is.

Are you able to say, in a similar fashion to above paragraph about variance, what value of AR coefficients says about the signal?


texts: Determining height of object from images

texts: Given two images of an object taken from the same distance but different vertical angles, is there a way to determine the physical height of the object?


texts:
I am trying to find the PSD of a 2D set of data and representing it as a function of frequency in both these directions. I start with a 2D FFT and generate the 2D hanning window for a $M\times N$ sized dataset $X$ by

w_1 = hanning(M);
w_2 = hanning(N);
w = w_1*w_2';
y_m = fft2(X.*w);

Now, I am having difficulty in applying the scaling to this 2-D PSD. For a 1-D dataset the scaling is given by:

\begin{align}
S_1 &= \sum_{i=0}^{M-1}w_i\\
S_2 &= \sum_{j=0}^{N-1}w_j\\
\textrm{PS} &= \frac{2\lvert y_m\rvert^2}{S_{1}^{2}}\quad\text{for a 1D PSD}
\end{align}

I am not sure how the scaling should be modified for a 2D PSD. Should it be $PS = \frac{2|y_m|^2}{S_{1}^{2}S_{2}^{2}}$?


texts: Computational Complexity of Polyphase Resampling

texts:
does someone know the approximate computational complexity of a polyphase resampling algorithm?

I am trying to find a good argument to prefer a polyphase resampling of Scipy over the FFT resampling method. In its documentation for resample_poly() it says:

  This polyphase method will likely be faster than the Fourier method in scipy.signal.resample when the number of samples is large and prime, or when the number of samples is large and up and down share a large greatest common denominator.

My data meets those criteria. I resample 15 channels * 30k-50k samples from 10kHz to 300Hz and practice shows it is way faster than the FFT method.

Yet I would like to have a referable source on those criteria.

I found a few papers on the computational complexity of FFT algorithms to be around O(n*log(n)) with n=n_channel*n_samples.


texts: Inverse Fourier Transform of the real part of fourier transform, and inverse transform of the imaginary part of fourier transform

texts: How can I calculate the inverse fourier transform of the real part of the fourier transformation?


texts: What is the right way to extract particular frequency phase information from FFT?

texts:
Signals acquisition:

Two proximity sensors are positioned Orthogonal to each other. These would be measuring the vibrations of a rotating shaft.

Problem:

Each of these 2 signals has two dominant frequency components (2441 Hz & 762.9 Hz ) in it. I don't know the phase and amplitude information. I'm interested in finding the phase of the signal at 2441 Hz and 762.9 Hz. I did the following process to achieve amplitudes and phases in MATLAB

Applying FFT to the raw signal and calculating Amplitude and Phase.

%fft for X signal
fftx=fft(X,NFFT);
% absolute value of the signal
Xn = abs(fftx);
A = mean(Xn);  % I took the mean to get single value.
phase = unwrap(angle(fftx));
phase_mean = mean(alpha);

Filtering raw signal to specific frequency prior to applying FFT using Butterworth filter. ( In this example I filtered to 2441 Hz)

[b,a] = butter(2, [ 2440/(fs/2),2441/(fs/2)],'bandpass');
X_filtered = filter(b,a,X);
% Applying FFT to filtered signal
fft_filtered = fft( X_filtered,NFFT);
Xn_filtered = abs(fft_filtered);
A_filtered = mean(Xn_filtered);
phase_filtered = unwrap(angle(fft_filtered));
phase_filtered_mean = mean(phase_filtered);

I did the same process as point 2 to filter to another frequency say 762.9 Hz.
My sampling frequency and other things are defined as follows.

L = length(X); % 50050 is length of signal
NFFT = 2^nextpow2(length(X)); % Zero padding to nearest N power 2
%Define frequency axis
fs = 1e7; % Sampling frequency
df = fs/NFFT; % frequency resolution
dt = 1/df; % time resolution
X_detrend = detrend(X,0); % Removing DC Offset

I've been reading a lot of other posts on the FFT phase and Amplitude.But I'm confused with most of them. In my case, I need to extract frequency particular information (amplitude and phase) to do further  calculation. I'm still unsure if using the filters is right way? Will FFT and filters alter the signal information in any way? Could someone enlighten in understanding the efficient way to extract amplitude and phase information from a signal of particular frequency component.


texts: Different representations of frequency space of 2D image FFT

texts:
I'm learning images processing using FFT.  In my test example provided below the input pixel values are clamped 0-1 (0-255), but I do eventually want to process floating point heightfield pixel values.

The software I'm using (Houdini v18) provides forward and inverse FFT functions and as a base line I can successfully convert an image to frequency space and back.

However when I look at the 2D representation of the frequency space it looks different to any examples I've found online.

Test image:

This is the result of Houdini's FFT with zero frequencies at center (height offset is pixel intensity):

The function returns 2x "images" representing real an imaginary values.

This however is the frequency space representation I find everywhere online:

From what I understand the radial type FFT images represent frequency between 0-2PI in u and v and the pixel value is the magnitude.  And that most images are offset so 0 is centered.

Do I need to convert the real and imaginary components to frequency and magnitude and plot those?  If so how?

EDIT:  My end goal is to apply radial pass filters to the FFT so I just need to apply any spacial transforms to get the FFT to that state.


texts: How do I implement an adaptive thresholding algorithm for underwater sonar

texts:
I want to implement an adaptive thresholding algorithm in MATLAB for filtering data received by an underwater sonar receiver. The data received has an interactive noise component resulting from underwater noise and specular reflection. CFARD method is close, but it does not serve my purpose. I have to image the data so that i am able to see the object, on a screen, that is placed underwater within the scandepth of the sonar. Any help will be greatly appreciated.

EDIT:

It is an underwater environment. I am trying to threshold a signal that has been received from a sonar transducer after it has been reflected by a solid target, that is situated in the same environment as the transducer. The problems belongs to Underwater Acoustic Imaging sonar domain. Problem is, I haven't been able to model the underwater environmental noise. From what I have read till now about this topic, the noise model follows a $K$-distribution. Also the environmental noise is not additive in nature, rather it's interactive. Hence the threshold has to be adaptive. I also mentioned CFARD method in my question. That is useful for signal processing in radar applications as we're just interested in finding a single point in a large area which has high energy. Same cannot be said about underwater acoustic imaging sonar, where we try to display the target on the screen as a video. I hope I have made it more clear now.


texts: How to perform a gaussian blur using FFT?

texts:
I am trying to implement a gaussian filter with matlab.Here is my implementation:
So far All I get as output is  a black image.
Any hints?

clear all;
close all;
I = double(imread('Put your path here'))/255;
I=rgb2gray(I);
sigma=1;
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[fx,fy]=meshgrid(f1,f2);
x=exp(-2*pi*pi*sigma*sigma*(fx/M.^2+fy/N.^2));
If=fft2(I);
If=fftshift(If);
If=If(:,:,1)*x;
If=ifftshift(If)
I=ifft2(If);
imshow(I)


texts: Calculation of an autocorrelation function

texts:
A sample of a random process is given as:

$$
x(t) = A\cos(2\pi f_0t) + Bw(t)
$$

where $w(t)$ is a white noise process with $0$ mean and a power spectral density of $\frac{N_0}{2}$, and $f_0$, $A$ and $B$ are constants.
Find the auto-correlation function.

Here's my attempt at a solution:

Let $a = 2\pi f_0t$, and $b = 2\pi f_0(t+\tau)$

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{x(t)x(t + \tau)\right\}\\
& = E\left\{\left(A\cos(a) + Bw(t)\right)\left(A\cos(b) + Bw(t+\tau)\right)\right\}\\
& = E\{A^2\cos(a)\cos(b) + AB\cos(a)w(t+\tau) + AB\cos(b)(wt)\\&\quad + B^2w(t)w(t+\tau)\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{AB\cos(a)w(t+\tau)\right\} + E\left\{AB\cos(b)(wt)\right\}\\&\quad + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(R_w(\tau)\right)\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

The expectation terms with the noise in them all equal $0$ (the last is just the auto correlation of white noise ... hence the simplification above. Using trigonometric identities:
$$
\cos(a)\cos(b) = \frac 12\left[\cos(a + b) + \cos(a - b)\right]
$$

we have:

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = E\left\{\left(A^2\right)\frac 12\left[\cos(a+b)+\cos(a-b)\right]\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = \left(\frac{A^2}{2}\right)\left[E\{\cos(a+b)\} + E\{\cos(a-b)\}\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

We're dealing with constant terms, so expectation term goes away and subbing in our initial conditions we get:
$$
\frac {A^2}2 \left[\cos(2\pi f_o(2t + \tau) + \cos(2\pi f_o\tau)\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))
$$

For some reason I can't help but feel I did something incorrectly calculating that autocorrelation ... it's supposed to be a function of $\tau$, but has a $t$ is in there ... I would very much appreciate it if someone could point me in the right direction, or explain what I messed up. I don't know whether it matters, but in this class we're dealing with only wide sense stationary processes.


texts: Approximate a System Frequency Response with a Filter in MATLAB

texts:
suppose I know the frequency response of a (linear) model approximating a real physical system but only at a specific frequency $f_0$ (so basically I have a complex number whose module is the magnitude response and the phase is the phase response).

In Matlab I want to filter a sinusoidal signal whose frequency is $f_0$ by a filter that approximates the response of my system at $f_0$ without altering the other frequencies. The filtering process must be done simulating the on-line behaviour, assuming that the next samples of the signal are unknown.

In other words my question is this: how can I build a filter whose behaviour at $f_0$ is the same of my system and does not modify the other frequency components?

The motivation I am asking this is that I need to perform some transformations on the frequency response at $f_0$ of the real system exploiting the results of the Singular Value Decomposition performed on the frequency response at $f_0$ of the identified (linear MIMO) model that approximates the real system.

Thanks in advance


texts:
I'm trying to compare the estimation of the power spectral density (in this case of a fake cosine signal) by two different methods: 1st, MATLAB's fft function, and 2nd, a direct calculation from the definition, $\int f(t)\cdot e^{-i\omega t}\mathrm{d}t$. The problem I find is that, for whatever reason, the FFT result is slightly shifted towards the right in frequency, and I'm not able to see why.

I'm adding the full code so that if anyone wants to try it, it should directly work and plot the results in MATLAB, something like

>> compare_FT(10);

As is clear from the code, the peak should be exactly at $\omega = 6$: that is indeed the case for the calculation from definition, but not exactly so for the FFT.

function[] = compare_FT(fs)

%% signal parameters (now fake, then from lab)
f = 6/(2*pi);
w_d = 2*pi*f;
T = 100;
num_samples = T*fs;

%% vector definitions
t = linspace(0, T, num_samples);
xx = cos(w_d*t); % the "fake" signal
dt = t(2) - t(1); % or, equivalently, dt = 1/fs
% omega vector for DTFT in a specific range:
ww_size = num_samples; % the bigger this size, the more points in the estimated FT
w_start = 5; % start calculating the FT at this omega
w_end = 7;  % stop at this omega
ww = linspace(w_start, w_end, ww_size);
Sx = zeros(1, ww_size); % prealocate space for integration

% DTFT method
for i = 1:ww_size
    int_value = trapz(xx.*exp(-1i*ww(i)*t))*dt;
    Sx(i) = int_value.*conj(int_value)/T; % estimate PSD
end

%% FFT method
n_fft = 2^14; % number of points in FFT (more points -> better FT discretization)
signal_fft = fft(xx, n_fft); % the fft itself
X2_fft = abs(signal_fft).^2; %signal_fft.*conj(signal_fft);
L = length(X2_fft);
Sx_fft = X2_fft(1:L/2+1)/(fs*num_samples); % proper scaling FFT to estimate PSD
w_fft = 2*pi*fs*(0:L/2)/L;

%% Plot results to compare both methods
figure(3);
clf;
box;
hold on;
% Plot theoretical value at top
plot(ww, ones(1,length(ww))*T/4, '--', 'Color', 'g', 'LineWidth', 3.5);
% Plot "direct" calculation
plot(ww, Sx, '.', 'Color', 'b', 'LineWidth', 3.5);
% Plot "FFT" calculation
plot(w_fft, Sx_fft, '-', 'Color', 'r', 'LineWidth', 3.5);
set(gca,'FontSize',18,'FontName', 'CMU Sans Serif');
xlabel('\omega', 'fontsize', 24, 'FontName', 'CMU Sans Serif');
ylabel('$S_x$', 'fontsize', 24, 'Interpreter', 'latex', 'FontName', 'CMU Sans Serif');
title('DTFT');
xlim([w_start, w_end]);
grid on;
end


texts: Correcting outliers using the standard deviation of a two-point different function?

texts:
I recently read a research paper that discussed how to correct lightcurve data (flux vs. time) by using the method as explained below:

  We have considered as outliers in the datasets the individ- ual
  measurements showing a point-to-point deviation in the two-point
  dierence function of the light curve greater than 3 for the SC and
  5 for the LC, where  is dened as the standard deviation of the
  two-point dierence function.

SC and LC are just different versions of the data.

Could someone explain to me what they mean by the two-difference function and how they are using it to correct the time-series data?

Edit:

Here's the link to the paper. You can find this quote on page 3, section 3.1.
http://arxiv.org/abs/1103.0382


texts: How to compute a sweeping cutoff filter in python?

texts:
I would like to compute a sweeping cutoff filter.
I wish to use SciPy.
I sliced my wave into x samples slices, in order to apply a different cutoff to each such created "frame".
The problem is that the filter, when applied to a slice filters the "clic" at the start of the slice, which creates an artefact for regathered result.
I searched stackexchange all over for sweeping filters questions and found nothing.
Any Idea on how to get a smooth sweeping without clic?
    #for filter
    from scipy.signal import butter, buttord, lfilter, freqz
    from scipy.signal import freqs

    def truncate_if_outside(data, bitwidth):
        if bitwidth==2:
            max_sample_value=32767
            min_sample_value=-32767
        else:
            print "Only 16 bit accepted by truncate_if_outside function"

        for i in range(len(data)):
            if data[i]>max_sample_value:
                data[i]=max_sample_value;
            if data[i]<min_sample_value:
                data[i]=min_sample_value;
            if math.isnan(data[i]):
                data[i]=0
            data[i]=int(data[i])
            if math.isnan(data[i]):
                data[i]=0
        return data

    def lowpass_filter(integer_data, sample_framerate, sample_bitwidth,threshold):#works with cutoff from 400 and higher

        def butter_lowpass(cutOff, fs, order):

            nyq = 0.5 * fs
            normalCutoff = cutOff / nyq
            b, a = butter(order, normalCutoff, btype='low', analog = False)
            return b, a

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

        result=butter_lowpass_filter(integer_data, threshold, sample_framerate, order=4)#what is order?
        result=truncate_if_outside(result, sample_bitwidth)
        return result

    def sweep_filter(integer_data, start_freq, end_freq, sample_framerate, sample_bitwidth, slice_width):
        num_slice= int(len(integer_data)/ float(slice_width))
        print "nombre de tranche dans cette durée :",num_slice
        new_data=[]
        diff_cutoff=end_freq-start_freq
        for i in range(num_slice):
            current_cutoff=start_freq+((i/float(num_slice))*diff_cutoff)
            print "cutoff",current_cutoff
            print "start:",i*slice_width
            offset_first_sample=(i*slice_width)-1
            offset_last_sample=((i*slice_width)+slice_width)-1

            slice_data=integer_data[offset_first_sample:offset_last_sample]
            slice_data=lowpass_filter(slice_data, sample_framerate, sample_bitwidth, current_cutoff)
            for j in range(len(slice_data)):
                new_data.append(slice_data[j])
        return new_data

The here filtered signal is a "smoothed noise" which means a random numbers smoothed with exponential curves. There should be no gap in it.