texts:
I know how to calculate the mean of a data set. In DSP, I think zero-mean should means that the error fluctuates around zero, and the mean of errors at any interval is zero. So what about the error in the picture that I show below? The noise in the picture doesn't look like a zero-mean error at all. But the overall mean of all the error values is zero (near zero). How do I mathematically say that the error in the picture is not zero-mean? Usually, which statistical tool is used for testing if a error has zero-mean?

(I think my problem related to the interval at which the mean is calculated. With the error I show in the figure, the overall mean is zero, but the mean at smaller interval is not zero. But I don't know if the official definition about zero-mean care about such thing.)

Edit: to be more correct, I'm asking about zero-mean error, not noise. The figure I show below is the error of my measurement (only error, not include the true value). The error is not stationary, not ergodic. The mean and variance of error change over time but following a repetitive patterns. However, the mean of error over a long period is near zero.


texts: Matlab: xcorr: Having trouble when I try to calculate delay between two sinusoids

texts:
I have a problem when I try to use Matlab built-in function for cross-correlation 'xcorr'.
Background: I created two sinusoids, one of them delayed in time. I choose center frequencies of these signals slightly different than another. In this way, I am trying to analyze the effects of frequency mismatch on cross-correlation results. Problem: The problem is when I try to calculate delay using xcorr function in MatLab, I could not find the exact lag. The result is quite far away from my offset-time. I try to increase my signal length or sampling frequency. But nothing changes.
Do you have any idea? Here is my simple code, below.
clc
clear all
close all

Fs=48e3;
dt=1/Fs;
t=0:dt:(4*Fs-1)*dt;
CarrierFreq=5e3;
Dev=50;
CarrierFreqDev=CarrierFreq +Dev*(2*randi(2,1)-3);

Offset=10000;
OffsetTime=Offset*dt;
NoisePower=-95
DesiredSNR=20;

% Signal Power
ReceivedSignalPower=NoisePower+DesiredSNR;
ReceivedSignalPowerWatt=10^(ReceivedSignalPower/10);
NoisePowerWatt=10^(NoisePower/10);

% Generate Noise @ specified power
Noise1 = wgn(1,length(t),NoisePowerWatt,'linear');
Noise2 = wgn(1, length(t),NoisePowerWatt,'linear');

% Generate Signal1 @ specified power
Signal1=cos(2*pi*CarrierFreq*t);
pSignal1 = sum(Signal1.^2) / length(Signal1);
Signal1 =  Signal1/sqrt(pSignal1);
Signal1_Norm=sqrt(ReceivedSignalPowerWatt)*Signal1;
Signal1=Signal1_Norm+Noise1;

% Generate Delayed Signal @ specified power
Signal2=cos(2*pi*CarrierFreqDev*t);
pSignal2=sum(Signal2.^2) / length(Signal2);
Signal2_Norm=sqrt(ReceivedSignalPowerWatt)*Signal2;
Signal_Delay=circshift(Signal2_Norm, -Offset);
Signal2=Signal_Delay+Noise2;

[Coeff, Delay_Time] = xcorr(Signal1,Signal2); %% with dev.
[ ~,Indx]  = max(abs(Coeff));
Delay_Time = Delay_Time(Indx)*dt;


texts: What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF?

texts:
I`m reading about Bilateral filter with images, I knew that it uses Gaussian filter (which is a LPF) as a domain filter plus a range filter and it is used to preserve edges(high freq component).

Can we consider Bilateral filter as a LPF ? .

I searched some articles:
1- Which approach is better for decomposing a image into high frequency and low frequency component?

The answer was that bilateral filter attenuates ""medium frequency "

2- How to Extract High Frequency and Low Frequency Component Using Bilateral Filter?

The answer was to consider bilateral filter as a LPF.


texts: Finding optimal detectors by numerical optimization (Canny)

texts:
In "A Computational Approach to Edge Detection", John Canny explains that the optimal "detector" for finding step edges profiles in an image via convolution is a non-closed-form function closely approximated by the first derivative of a Gaussian. He derives this by defining "optimal" as the function which maximizes the product of signal-to-noise ratio and localization (ability to detect the edges at the right location), as well as another criteria concerning multiple responses.

On p.5 of the paper, he also shows in a picture the optimal detectors, under that those criteria, for ridge and roof profiles, which also look Gaussian-ish in shape.

I'd like to compute the optimal detector for other, arbitrary profiles. Unfortunately, I'm not quite well-versed in numerical optimization enough to figure out how to do that based on the formulas in the paper, and I'm not able to find any resources that explains how to do so explicitly. Does anyone know of a guide or an implementation that will do that - finding the optimal detector for an arbitrary feature in an noisy image?


texts:
I'm trying to extract the code of my garage remote, but i have some difficulty to find witch modulation has been used.
At first I thought it was a simple digital modulation, but there is some strange situation that it seems that the signal changes phase.
Here are some pieces of signal extracted with my SDR dongle:
f.Carrier = 868.25Mhz
T.rectangle = 400us

Here is the spectrum of the signal:


texts: How to determine the order from Auto-Regressive model in Matlab?

texts:
I have time series vibration signal recorded for 10 seconds. I segmented the vibration signal into a number of frames. I need to use an Auto-regressive model for this signal. AR model is used to extract the data from the vibration signal. The goal of an AR model is to estimate the AR coefficient  that can fit the original data as much as possible through an optimization process. I am going to classify the vibration signal using the AR coefficient as the features. In order to find the effective classification accuracy, frame analysis has been done. AR model coefficients are extracted from each frame.

My problem is, how can I develop the AR model? How do I determine the order of polynomial from the AR model? And also how do I find the coefficient?

Perhaps, can anyone provide me the Matlab code for this?


texts: How to cut off the initial disturbance automatically?

texts:
I would like to find a way to cut off the initial disturbance (e.g. the part tf_n<20). Though it is easy to do manually for single figure, I need a way to automatically identify the steady point, as I need to do it 100 times otherwise.

I would prefer to do it by Matlab.


texts: Designing a 2d low pas filter in fourier domain for an image

texts:
Firstly this is not a homework question. I am working on a project and my professor gave me this task so I can also learn about filters and fourier domain.

I have a picture that is 2800 x 4500 pixels. I need to make a filter that will only give an arbitrary number of "frequencies".

My thought process is this:
1. Transfer the picture in fourier domain using

F = fftshift(fft2(img));

2: Multiply it by a filter that has ones around the centre and zeroes elsewhere.

3: Transfer the image back into x-y domain

The problem is that I don't know how to make the filter since the picture is non-square. I tried making a circle but that didn't work. How do I make an ellipse shape?

In 1-d it is easy because the filter is a vector with ones around the centre.

I tried looking for answer to this but the pictures in question are always square pictures so I cannot apply the same thing. Also nowhere does it say how it works but rather it is just done.


texts: Coloring a waveform with spectral centroid or by other means

texts:
Update 8

Unhappy of having to upload tracks to the service and looking at the new release of RekordBox 3 I've decided to take another look for an offline approach and a finer resolution :D

Sounds promising though still in a very alpha state:

Johnick - Good Time

Note that there's no logarithmic scale nor palette tuning, just a raw mapping from frequency to HSL.

The idea : now a waveform renderer has a color provider which is queried for a color for a specific position. The one you are seeing above gets the zero crossing rate for the 1024 samples next to that position.

Obviously there's still a lot to do before getting something robust but it seems like a good path...

From RekordBox 3:

Update 7

The final form that I'll adopt, much like in Update 3

(it's been a little Photoshopped to get smooth transitions between colors)

Conclusion is I was close months ago but did not consider that result thinking it was bad X)

Update 6

I've unearthed the project recently, so I thought about updating here :D

Song : Chic - Good Times 2001 (Stonebridge Club mix)

It's way better IMO, beats have a constant color etc ... it's not optimized though.

How ?

Still with http://developer.echonest.com/docs/v4/_static/AnalyzeDocumentation.pdf (page 6)

For each segment :

public static int GetSegmentColorFromTimbre(Segment[] segments, Segment segment)
{
    var timbres = segment.Timbre;

    var avgLoudness = timbres[0];
    var avgLoudnesses = segments.Select(s => s.Timbre[0]).ToArray();
    double avgLoudnessNormalized = Normalize(avgLoudness, avgLoudnesses);

    var brightness = timbres[1];
    var brightnesses = segments.Select(s => s.Timbre[1]).ToArray();
    double brightnessNormalized = Normalize(brightness, brightnesses);

    ColorHSL hsl = new ColorHSL(brightnessNormalized, 1.0d, avgLoudnessNormalized);
    var i = hsl.ToInt32();
    return i;
}

public static double Normalize(double value, double[] values)
{
    var min = values.Min();
    var max = values.Max();
    return (value - min) / (max - min);
}

Obviously there's a lot more code needed before you get here (upload to service, parse JSON etc) but that's not the point of this site so I'm just posting the relevant stuff for getting the above result.

So I'm using the first 2 functions of the analysis result, there is certainly more to do with but I still need to test. If I ever find something cooler than above I'll come back and update here.

As always, any hints on the topic are more than welcome !

Update 5

Some gradient using harmonic series

The color smoothing is ratio-sensitive otherwise it looks bad, will need some adjustment.

Update 4

Rewrote the coloring to happen at source and smoothed colors using an Alpha beta filter with 0.08 and 0.02 values.

A little better when zoomed out

Next step is to get a great color palette !

Update 3

Yellows represents mediums

Not so great when unzoomed, yet.

(the palette needs some serious work)

Update 2

Preliminary test using second 'timbre' coefficient hint from pichenettes

Update 1

A preliminary test using an analysis result from EchoNest service, note that it's not aligned very well (my fault) but it's much more coherent than the above approach.

For people interested in using this great API, start here : http://developer.echonest.com/docs/v4/track.html#profile

Also, don't get confused by these waveforms as they do represent 3 different songs.

Initial question

So far this is the result I get using a 256-samples FFT and calculating each chunk's spectral centroid.

Raw result of the calculations

Some smoothing applied (form looks much better with it)

Waveform produced

Ideally, this is how it should look like (taken from Serato DJ software)

Do you know what technique/algorithm I could use to be able to split the audio when average frequency changes over time ? (like the above image)


texts: Creating uneven sidebands with AM + PM modulation?

texts:
I am trying to modulate an optical signal such that I create a pair of sidebands with arbitrary phase and amplitude with respect to each other. For example I would be like to be able to modulate my carreir signal, centered at 100, to have one of the two resultant spectra after my modulation.

(Blue is the real component of the spectrum and yellow is imaginary).

So that the 110Hz component has a different phase and amplitude form the 90Hz component.

If I take the inverse fourier transform of these two spectra I get real signals which look like this.

Both of these two signals aren't unphysical, however I don't know how to analyse them in what way I'd need to modulate my carrier to create them. Is it possible to recreate these time signals with Amplitude and phase modulation alone?

Thanks! :)

Edit: In response to  Dan Boschen. This is what my spectrum looked like when I did phase modulation. I am using a large modulation "depth" of 0.5. When I reduce it the harmonic components do become negligible, so I guess in that case I do get just two sidebands.

Mathematically what I'm doing. Where fc is the carrier frequency and fm is the modulating frequency.

My python code to generate that signal:

timeStep =0.0001
t = misc.stepRange(0,1,timeStep)

#carrier freq
carrierFreq = 100

#modulating freq
modulatingFreq=10

#carrier= 10*np.cos(t*2*np.pi*carrierFreq)

modulatingSig= np.cos(t*2*np.pi*modulatingFreq)

modIndex = 0.5
resultantSig1 = np.cos(t*2*np.pi*carrierFreq+modIndex*modulatingSig)


texts: How Do I Measure the Time Duration of a Finite Length Discrete Sequence?

texts: Assume I have a five-sample time-domain sequence (none of the five samples are zero valued) and the time period between each pair of samples is one second. Measured in seconds, what is the time duration of my five-sample sequence?


texts:
Can you simulate the "response time" of a lower window size FFT on a higher window size?

That is, is it possible to calculate the equivalent lower window sized FFTs from a higher one in order to iterate over the equivalent e.g. 128 window sized FFTs, but plot the higher one e.g. 8192 samples window size?

This is a problem in e.g. when one is sampling at a higher window size for the FFT, but some signal parts (e.g. strong amplitude peaks in audio) appear at a faster pace and are thus "missed" by the FFT. Or not necessarily the FFT, but any processing that uses that higher window size FFT.

What I therefore need is the visible resolution of a higher window size and the faster response of a lower window size.

My application (in audio) requires me to be able to track frequencies that follow closely to the beats of the audio. If e.g. a drum loop plays at a high bpm, then a 8192 samples FFT might be too slow to be done reading form since it changes slower than the drum loop's frequencies.


texts: determining type of filter given its pole zero plot

texts: How can I classify a filter given its pole-zero map. For example I've got my zero's located at $\pm j$ and my poles located at $\pm\frac{1}{2}j$.


texts: Is there a relation between the dead zone of sigma delta modulator and filter order?

texts:
This question refers to the link https://pure.tue.nl/ws/portalfiles/portal/1546321/521243.pdf. In chapter 4 it discusses the topic of dead zone, which appears in the delta sigma modulator. Dead zone means that when input signal to the modulator is very small than the output of the modulator will not follow the input, for small input signal the output of the modulator depends on the state of the modulator. The input must be higher than the dead zone to break the repetative pattern of the modulator.

The question: How the the dead zone is related to the filter order of the modulator? In the link it shows for first order filter but not much is written about higher order modulator and dead zone.


texts:
I have an observed signal which is the norm of the sum of two 3-vectors:
$$
S(t) = || \mathbf{A}(t) + \mathbf{B}(t) ||
$$
Now, I have full knowledge of the vector $\mathbf{A}(t)$,
$$
\mathbf{A}(t) = \pmatrix{A_1(t) \\ A_2(t) \\ A_3(t)}
$$
i.e. I know each of the components $A_i(t)$. The signal $\mathbf{B}(t)$ is unknown, and $S(t)$ is what I observe. Now the signal $\mathbf{B}(t)$ has peaks, particularly in the $B_1(t)$ component, and my goal is to determine the position of these peaks. I don't care about the amplitude of the peaks, I only want the position (time).

So what I need to do is somehow filter out the signal $\mathbf{A}(t)$ from the observed signal $S(t)$ and then apply a peak detection algorithm on the residual signal. My question is whether it is possible to do this?

As a simple example, I made a matlab script to construct synthetic data as follows:
$$
\mathbf{A}(t) = \pmatrix{ X(t) \\ Y(t) \\ 0}
$$
$$
\mathbf{B}(t) = \pmatrix{ G(t) \\ 0 \\ 0}
$$
where $X(t) = \sin{(2 \pi (f/4) t)}$, $Y(t) = \sin{(2 \pi f t)}$, and $G(t)$ is a single Gaussian peak. The signals and matlab script are given below. As you can see, it is not easy to tell that there is a peak in my observed signal at $t = 500$ as there is in the unknown signal $G(t)$.

So is it possible to filter $S(t)$ to remove $\mathbf{A}(t)$ and clearly find the peaks in $\mathbf{B}(t)$?

N = 1000;
t = [1:N];

f = 20;
X = sin(2*pi*(f/4)*t/N);
Y = sin(2*pi*f*t/N);

G = gausswin(N, 100.0);
Z = zeros(1,N);

A = [ X ; Y ; Z];
B = [ G' ; Z ; Z ];
C = A + B;

signal = vecnorm(C);

figure

subplot(4,1,1);
plot(t,X);
title('Known X');

subplot(4,1,2);
plot(t,Y);
title('Known Y');

subplot(4,1,3);
plot(t,G);
title('Unknown signal G');

subplot(4,1,4);
plot(t,signal)
title('Observed signal');


texts:
Does extreme clipping white noise colour the audio spectrum?

I thought (and hoped!) not on the basis of this from research.net

  Clipping is mathematically equivalent to the addition of impulse
  noise, (at the point of clipping you are effectively adding a negative
  impulse to the original waveform); consequently the effect on the
  spectrum is similar to additive impulse noise, with a fairly uniformly
  flat spectrum (depending on the degree of clipping). The spectrum is
  obviously going to be coloured if the clipping occurs with some
  periodicity

Emphasis mine.


texts:
I'd like to know which conditions must be satisfied in order to get a signal with a finite bandwidth.

Thank you for your time.


texts: Calculating maximum signal frequency knowing its sampling rate?

texts: I was curious as to how I could calculate the maximum frequency of a signal knowing only its sampling rate (for example 125 nanoseconds). Does the sampling theorem have anything to do with it (sampling rate is at least two times the maximum frequency)? So is the max frequency half of the sampling rate?


texts: inquiries for writing continuous wavelet transform codes manually

texts:
I want to write continuous wavelet transform  codes manually by matlab. And I want to use complex morlet function. Here are some background:
Continuous wavelet transform definition: $C(S,T;f(t),\psi (t))=\frac{1}{\sqrt{S}}\int_{S}^{b}f(t)\psi^{\ast }(\frac{t-T}{S})$
S is scale vector.for example 1:60;$T$ is time slid.
And $\psi(t)= \frac{1}{\sqrt{\pi f_{b}}}e^{2\pi f_{c} t}e^{\frac{-t^{2}}{f_{b}}}$

psi = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*t).*exp(-t.^2/fb); % for example fb=15;fc=1;

my discrete signal has $N$ points. discrete version of that integral for the first element in the scale vector  is:  for s(i)
$C_{T}(s)=\sum_{n=0}^{N-1}f(n)\psi^{*}(\frac{T-n}{S})$ this must calculate for all scales,1:60;
at the end it must return me a complex matrix with $N\times S$ size.
I am confused about wrting this codes. Anyone can help?
P.S I don't want to uses matlab function conv2 to calculate that convolution.
Thanks in advance.
 here is my first attempt, but it's not working at all.

%% user CWT
clear all
N=300;                %sample point numbers
t=linspace(0,30,N);
%% signal
x=5*sin(2*pi*0.5*t);  % signal with freq of 0.5 HZ
%% cwt
    fc=1;fb=15;
% psi=((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
%     t).*exp(-t.^2/fb);
%% convolution Psi([N-n]/S)*x(n) so we calculate convolution(psi(n/s),x(n))
    for s=1:60   %scale vector s=[1:1:60]
for i = 1:N      % number of discrete times
       for k = 1:i
               if ((i-k+1)<N+1) && (k <N+1)
                PSI = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
       (t/s)).*exp(-(t/s).^2/fb);
                        c(i,s) = c(i)+ x(k)*PSI(i-k+1);
                end
       end
end
    end


texts: How to compare two 802.15.4 O-QPSK modulated digitals signals?

texts:
I am working on wireless sensors project and I am trying to simulate a transmission between two nodes.

So let's say I have two nodes which transmit the message to each other, i.e. node 1 transmits to node 2 and then node 2 transmits to node 1.
Let's say the message being transmitted is exactly the same and the difference between the packets is just Source/Destination Address as well as Frame Check Sequence Field(FCS) of MPDU. The message is being generated and encoded into signal wave.

Now, without demodulation process and  is there any way to distinguish two signals and see if it comes from node 1 or node 2 by looking only at the waveforms of transmitted signals? Observer sees signals separately, i.e. they are not mixed and come to receiver one after another with short interval in time between them. So it is basically to classify them as "signal 1" and "signal different from signal 1". Is there also a way to extract bits from the signal ?

For example possibly  by converting it to frequency domain by looking at the amplitude or phase spectrum of the signal?
Or may be it is achievable after certain number of iterations so enough data can be sampled and thus pattern can be drawn?

Please share your ideas or suggestions.


texts: How to make the $\ell_2$ norm of all columns and rows of an $n \times n$ matrix equal to $\sqrt{n}$?

texts: I have an $n \times n$ matrix and I would like its columns and rows to have $\ell_2$ norm equal to $\sqrt{n}$. Is this possible?


texts: how to find frequency response of microphone

texts:
hello I want to find the frequency response of a microphone.  I give the input signal to my speaker and it produces a specific SPL with specific frequency. on the other side, I read the microphone output voltage.
but I wonder which value of SPL is required to find the frequency response of microphone?
should I fix the SPL and just change the frequency and read the microphone output? if so, which SPL value I must choose? should I fix the SPL to 94 dB SPL?
you can assume that I have a really precise reference microphone and I can find out the SPL of speaker by means of this reference microphone. I also don't really care about direction and distance because I use the reference mic so the only thing that matters to me is SPL. also assume that I afford very precise speaker too and don't care about it's response. I appreciate your help
thanks


texts: Area under waveform in t-domain including 2 delta functions

texts:
question:

my attempt:

$(a)$ $\displaystyle\int_{-\infty}^{\infty} f(t) dt=2 u (t)+r(t)-r(t-T)+Ku(t-T)$

then we draw it's diagram and found it was unbounded so,we doubt how area can be zero for any value of $K$ so, we need help here finding value of K (which we think  will come negative )

$(b)$ $\displaystyle\int_{0_{+}}^{\infty} f(t) dt=r(t)-r(t-T)+Ku(t-T)$

same diagram as part $(a)$ , and same story as above ,except  we pull above curve down by two units .

so , we're unable to find value of $K$. any help will be appreciated

thanking you !

edit : $r(t)$ is ramp function $= t u(t)$


texts:
I am reading a book on filter design and the author mentions to evaluate a transfer function, let 2 * PI * f (w) vary from 0 to PI to get the first half of the response.

How exactly is this done in code and what is meant by vary in terms of actual values?

Should f not change as well to get the response of different frequencies?

How do the calculations look if there is a varying f and a varying 0 to PI?


texts: What does a "vector" in a hidden Markov model mean?

texts:
I know that a Hidden Markov Model (HMM) is used in speech recognition and understand it to some degree. However, what I don't know is how input (speech) is "transformed" to a vector which in later used in HMM.

How do you get a vector from a sound input? Is this vector readable by a human?


texts: Digital filter coefficients from low-pass to high-pass

texts:
Given I have coefficients a0, a1, a2, b1, and b2, defining the difference equation for a digital filter as:
y[n] = a0 * x[n] + a1 * x[n - 1] + a2 * x[n - 2] - b1 * y[n - 1] - b2 * y[n - 2]
Which defines a low-pass filter with particular cutoff frequency, how can I obtain the coefficients A0, A1, A2, B1, B2, which similarly define a high-pass filter with the same cutoff frequency? I'm aware there are so-called "bandform transformations" for converting a prototype low-pass into a high-pass, but to my knowledge, these are not directly applicable to discrete-time/digital filters, so I am unaware of any way to apply them to this problem.
If these coefficients are derived from a complex-conjugate pair of zeros and/or of poles given the discrete transfer function for the low-pass filter would be:
H(z) = (z - Zero[0]) * (Z - Zero[1])/[(Z - Pole[0]) * (Z - Pole[1])]
Is there then a way to transform this function to the corresponding high-pass filter I'm looking for to get the poles and zeros from the new transfer function?


texts: Whats the technical term for a filter that allows a non unity stride when sliding it's window?

texts:
So a 1D mean filter of size two with a sliding window defined as W below could be considered as having a stride of one.

Wn=(Xn + Xn+1) / 2
Vn=W2n

where as if I were to define a 1D 'stepped mean filter' with a stride of two and labeled U below.

Un= (X2n + X2n+1) / 2

where the number of samples in U (and also V) would be half that of W.
My questions are the following

Is U still a filter (I'm guessing it is even though not fitting with the traditional kernel convolution implementation I associate with filters)
What is the technical language one could use to communicate and research questions regarding filter/non-filter such as just given?
All though not directly related but since I have the floor is, if  V and U would produce equivalent results (my intuition says no)?


texts: OFDM Modulation: what's the output of the FFT?

texts:
I'm trying to understand how OFDM is used in technologies like Wi-Fi.
On the Wikipedia page, an idealized transmitter is presented with this schema:

Data is converted from serial to parallel, mapped to some complex symbols using constellation mapping, then, using the inverse fast fourier transform, is converted to what seems one complex number with a real and imaginary part, which is then converted to some continuous signal using DAC.
On some other diagrams, the output of the IFFT is the same size as the input, and all of the outputs are multiplied by some (orthogonal) carriers before being transmitted, e.g. on this diagram (in French, but it should be clear what is what):

I'm having trouble understanding how these two diagrams relate with each other: are they different systems ?
As the FFT takes as an input a sampled continuous signal (so, let's say $N$ complex numbers), the IFFT should output $N$ complex numbers, but this does not appear on the Wikipedia diagram. Am I missing something?


texts:
Are there any "good" low-computational-cost approximations for parameterized Kaiser window generation suitable for small systems/languages that do not include Bessel functions in their standard library?

(maybe only to 8 to 12-bit accuracy, and something that can be run in Basic on an Apple II, et.al.)


texts:
I have a 3D binary image that I must fill all of the holes in the image.

"Fill holes 3D" in this case, is basically a 3D version of

BW2 = imfill(BW,'holes') from MatLab,

but I am not focusing on Matlab here. I am focusing on the algorithm in general.

Currently, I simply perform a "flood fill 3D" algorithm on all voxels at the border area of the 3D matrix. Anything not filled are holes so take this as a mask and remove all voxels on the real image with this mask.

Sure, the algorithm works, but somehow, I feel this is inefficient.

What is the fastest way to "fill holes" in 3D in this case?


texts: Group delay of $H(\omega)= 1- re^{j \theta}e^{ - j \omega} $

texts:
I'm studying chapter 5 of Discrete-Time Signal Processing 3rd edition by Alan Oppenheim and I'm having serious difficulties understanding how he obtained equation 5.57. For those who don't have this book I tell you that in this part it is analyzing the frequency, phase and group delay of

$$ 1- re^{j  \theta} e^{-j \omega} $$

which could either be a pole or a zero depending on whether this factor is in the denominator or numerator of the frequency response. Here r is a random magnitude variable and theta is a random phase variable

So far I have been able to understand how is the phase expression obtained as it is defined as $$\arctan  \frac {\Im}{\Re}$$

and

$$ 1- re^{j  \theta} e^{-j \omega} = 1- r[ \cos (\theta - \omega) + j \sin(\theta - \omega) ]
$$

and as cosine is even and sine is odd

$$ = 1- r[ \cos (\omega- \theta) - j \sin(\omega- \theta) ]
= 1- r\cos (\omega- \theta) + j r\sin(\omega- \theta)  $$

the resulting phase expression is:

$$ \arctan  \frac {r \sin(\omega- \theta)}{1- r \cos (\omega- \theta)}$$

which coincides with equation 5.56 in the book

But when it comes to finding the group delay (which is the negative derivative of this expression) I'm not obtaining what it says in the book. Moreover I introduced the expression in Matlab and I'm obtaining the following result:

According to the book the group delay is;

$$ \frac{r^2 - r\cos (\omega- \theta)}{1 + r^2 - 2r\cos (\omega- \theta)}$$

How did they get there? Can you help me?


texts: Positive or negative sign on Fourier transform formula

texts: I have seen both the formula of Fourier transform with positive and negative sign on exponential as  $$ X(\omega)=\int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt$$ and $$ X(\omega)=\int_{-\infty}^{\infty} x(t)e^{j\omega t}dt$$ I am confused which one is the correct formula. I also solved for Fourier transform by taking the following example $$x(t)=\begin{cases} 1, \hspace{5mm} \text{for} \hspace{2mm} |t|<1 \\0, \hspace{5mm} \text{for} \hspace{2mm} |t|>1  \end{cases}$$ and got the same result as $$ X(\omega)=\begin{cases} 2\frac{\text{sin}\omega}{\omega}, \hspace{5mm} \text{when} \hspace{2mm} \omega \neq 0 \\2, \hspace{13mm} \text{when} \hspace{2mm} \omega = 0\end{cases}$$ Can anyone explain whether both the formula for Fourier transform are correct or not?


texts: Modeling channel response as complex vs splitting I and Q

texts:
If I have a communications channel where a digital QAM signal is transmitted, I can either model the channel's impulse response as a complex $h(t)$, or as 2 separate impulse responses $p(t)$ and $q(t)$ for the I and Q channels, respectively. Let $s(t)$ be the transmitted signal, with $s_I(t)$ and $s_Q(t)$ denoting its I and Q components.

For the complex case, the filtered signal would be $$[s_I*\operatorname{Re}(h)-s_Q*\operatorname{Im}(h)]+j[s_I*\operatorname{Im}(h)+s_Q*\operatorname{Re}(h)]$$

For the 2 separate case, it would be $$[s_I*p]+j[s_Q*q]$$

Clearly, they look rather different. Is either one a "better" model than the other? Better in the sense of, for example, using an adaptive filter to equalize the channel. We can train for a complex response with the cross coupling of the former, or run the channels separately. If not, can we talk about some sort of mathematical equivalence between them? Is there a mathematical mapping?


texts:
I am working on an accelerometer on a mobile device. When an event occurs, I record around 30 acceleration values. The accelerometer, being extremely noisy is very unreliable. What makes it worse is that my motion during the acquisition time is also of the same order as the noise. Is there a good method to estimate the curve the mobile goes through with a decent accuracy?

Edit: I am trying to estimate the position of the mobile. Further, the amplitudes of acceleration are of the same order as the noise(I logged the accelerometer values when the mobile was stationary and made this estimate). Further, The noise is approximately gaussian.


texts: How to interpret the signal in frequency domain

texts:
I'm having difficulties understanding some basics regarding FFT. I was hoping you could answer a few questions I have. These questions will probably seem silly to you, sorry for that. I'm a complete beginner in DSP.

Let's consider the following script:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, np.pi * 2, 256)
y = np.sin(x) + np.sin(2*x + 1)
Y = np.fft.fft(y)

freq = np.fft.fftfreq(x.shape[-1])
plt.plot(freq, Y.real)
plt.savefig('Y_real.png')
plt.close()

plt.plot(freq, Y.imag)
plt.savefig('Y_imag.png')
plt.close()

y2 = np.fft.ifft(Y)
plt.plot(x, y2.real)
plt.savefig('y2_real.png')
plt.close()

The signal in time domain (y2_real.png):

The signal in frequency domain, real part (Y_real.png):

Imaginary part (Y_imag.png):

I was expecting that the real part will represent amplitudes of two sinusoidal signals and the imaginary part will represent a phase shift of them. However the amplitude 109 doesn't make much sense to me, nor does the phase shift of +/-128.

Where did these numbers came from? Also I'm not sure whether everything is right with OX axises in the frequency domain. How exactly can I figure out which frequency represents one peak or another?


texts: Phase difference between two sinuisoidal signal

texts:
Currently, I have two signals, the main components of both signals are 60Hz, but both also have weaker response at 180Hz + small amount of noise. As shown in the photo below, I want to find the phase difference between the 60Hz component of first signal, and the 180Hz component of the second signal.

As you can see the second waveform of the figure, the phase lagg is of the 180Hz signal is 1.39ms, which coresponds to a 90 degree phase shift.

I am wondering, if I want to implment this phase measurement, could I simply take the fft of two signals, and then find the phase at 60Hz of signal 1, and subtract the phase at 180Hz of signal 2?

Or if there is better way to do so?


texts:
I'm trying to make a program that catch and audio signal of a bass and show the notes that have been played. My program:

Gets the sound and put it in a double data type array. (48 000Hz, 16bits) (1 frame in array per sample)
Splits the array into many arrays, with 4~5 arrays in a second approx.
Applies FFT algorithm to every array and left the same number of array's that have been created in the last step with frequency domain instead of time domain.
Obtains the highest number of each array. That's the note that I have played with the bass.

This works and get the frequency when I don't split the array, for most of the times, but when I split it doesn't work even for the half of the results, maybe I have to use another way for get several notes instead of splitting. But the results are good too, so I think the problem is that I have to apply some filters.

The first one has to be a filter for the range of the bass (41 - 784), assuming that I'm not playing anything that time if it's out of that range.

Many times, if I play 41Hz (E1), the program interpret that I'm playing 82Hz (E2), or even a higher harmonic, and I don't know what to do with that. I could turn into 0Hz every note that has been played after the same note or one of its harmonics, but then how could I play the same note several times.

Some advice for deal with this problems?
Should I split the original array into more/less array's for get more accuracy?
Which filters should I apply on my frequency or time domain array for get all the notes that have been played, with the right frequency and avoiding their harmonics?


texts: Correct Radar Range Equation for SNR at Detection

texts:
Assuming you'd like to estimate the probability of detection of a target for a monostatic radar.

For estimating detection SNR of a monostatic radar for given target rcs signature, do you use the energy or power form of the SNR equation?

Budge's Radar Equation Notes

Energy form of the SNR equation:

Power form of the SNR equation:


texts: How to find a 'square root filter' such that $d*\overline{d}=f$ given $f$?

texts:
I have a desired frequency response $F(\omega)$. I know how to create a FIR filter (say, with linear phase) $f(n)$ with this frequency response.

But how can I get a filter $d(n)$ such that $D(\omega)D^*(\omega)=|D(\omega)|^2=F(\omega)$, or equivalently in the time domain, $(d*\overline d)(n)=f(n)$, i.e. the 'square root' of $f$? (here $\overline d$ is the filter $d$ with reversed order of coefficients)

I would like to get a real filter with compact support, preferably using MATLAB.

Thanks.


texts:
In most of the textbooks, the motivation for the usage of DFT in OFDM goes as follows. Let $X[0],\cdots,X[N-1]$ be complex symbols, and
$$s(t)=\sum_{k=0}^{N-1}X[k]e^{j2\pi\frac{k}Tt}\qquad (0\leq t\leq T)$$
be the complex baseband signal with duration $T$. Then, we sample $s(t)$ at points $t=mT/N$ to obtain
$$s[m]=s(mT/N)=\sum_{k=0}^{N-1}X[k]e^{j2\pi\frac{km}{N}}\qquad m=0,1,\cdots,N-1$$
Then, most of the textbooks usually say that since $\{s[m]\}$ is, up to some constant factor, the inverse DFT of $\{X[k]\}$, we may use IDFT to modulate the symbols $\{X[k]\}$ and then pass the result to D/A converter. At the receiver side, if there is no channel distortion, DFT can be used to demodulate the signal and get back the transmitted symbols $\{X[k]\}$.
The following are my questions:

Although we can sample any signal we want, usually one of many goals is to reconstruct the continuous signal back later from the samples. However, $s(t)$ is not even bandlimited, since it is timelimited. Therefore, the sampling theorem does not apply. Then, what is the point of sampling here? Is it just for motivation for the usage of DFT?

After the IDFT and the D/A converter, what does the resulting signal look like? Ideally and mathematically, I think we should obtain

$$s'(t)=\sum_{m=0}^{N-1}s[m]\mathrm{sinc}(t/T-m)=\sum_{k=0}^{N-1}X[k]\Bigl(\sum_{m=0}^{N-1}e^{j2\pi\frac{km}{N}}\mathrm{sinc}(t/T-m)\Bigr)$$
At least, mathematically, this is the most powerful interpolation method and should give me back the continuous signal under some conditions. However, I don't see why it gives me back $s(t)$

Let us say I just pretend IDFT to be some coding algorithm, and just transmitted the (complex baseband) signal $s'(t)$. I think it is okay. Then, what is the spectrum of $s'(t)$? Whatever it is, it is bandlimited. However, usually, in most of the textbooks and pictures, spectrum of OFDM signals are shifted sinc functions, which are not bandlimited. What is going on here?

I am sorry if I have asked three questions. I do think they are intimately related. Thank you very much.


texts: Amplitude of the Signal in Frequency Domain different from Time Domain

texts:
How does the Amplitude of the Signal is changed when taking the FFT of a Signal,

Have a look,

The amplitude had changed from 2 to 30,

and Here is my code for generating the above output,

f=1e3;   %Frequency of Wave
A=2;    %Amplitude
Fs = 1e6;  %Sampling Frequency

Ts = 1/Fs;   %Sampling Rate

t = 10/f;       %Time period of 1 Oscillation = 1/f

n = 0:Ts:t;  %Generating Samples

x=A*sin(2*pi*f*n);
subplot(2,1 ,1)
plot(x

% 100,000 = Fs
% 10,000 = Length of the Signal
% 100,000/10,000 = 100Hz <- First point in FFT = 100Hz
% 2nd Point = 200Hz
% 3rd Point = 300Hz
% 4th Point = 400Hz
% .
% .
% 10th Point = 1KHz <- Original Signal Frequency

subplot(2,1 ,2)
F=fft(x);
plot(real((F))),grid on
xlim ([0 20])


texts: Bandpass butterworth filter in python is not working

texts:
I've implemented this bandpass filter in matlab for the audio signals and it's working as expected. Now, I tried to implement the same thing in Python but it doesn't produced expected results, it just sort of discard the whole signal.

Here is the Matlab code:

[f,fs] = audioread('abc.wav');
n=6; %order
beginFreq = (300/(fs/2));
endFreq = (2300/(fs/2));
[b,a] = butter(n,[beginFreq,endFreq],'bandpass');
%filter the signal
filtered = filter(b,a,signal);

Here is the Python code (using scipy.signal)

[originalSignal, sampleRate] = sf.read('abc.wav')
lower = (300/(sampleRate/2))
higher = (2300/(sampleRate/2))
n=6 #order
[b, a] = signal.butter(n, [lower, higher], 'bandpass')
filtered = signal.lfilter(b, a, originalSignal)
plt.plot(originalSignal, 'b', filtered,'r')

There is one more weird thing with the output! I tried to change the values of the lower and higher frequencies. As the frequency range increases; only the amplitude of output signal increases from 'almost zero' to higher value. For example, when I give the lower and higher frequency values as '100.0' and '13000' respectively, the output signal is "same" as the input signal but with almost half the amplitude. So, apparently the filter is not affecting the frequencies at all but rather affecting the amplitude.....
I've no idea what's going on... :(

Filter results:


texts: Determining the fundamental frequency/pitch of a note

texts:
What I am trying to achieve is getting the fundamental frequency of a note played by an instrument. What I have already done is performing an FFT on a samples of audio file, and here's what I get:

After that I just tried to find the maximum value in that dataset and I get "4138". Unfortunately, the note being played in the audio file is an E and from internet I know that the E note should have 5274.04 Hz or 2637.02 Hz.
I continued reading some additional info on the internet and found interesting info that in some cases fundamental frequencies could be missing from the sound, so I have tried doing some autocorrelation on the samples I got from the audio file to find the fundamental frequency and I got something like this:

It looked like a frequency domain for me so I performed FFT on that dataset but got similar results.
Also I have read that some people suggest using HPS (Harmonic Product Spectrum) for finding the fundamental frequency, I have tried that but obviously the result was the same, as the FFT itself didn't change.
P.S. Of course I use Hamming windowing function for samples. The library I use for FFT is JTransforms.


texts: Decompose a Morlet wavelet in a sine wave to its gaussian and sine components

texts:
I have a sine wave where there is a morlet wavelet inside:

Now, first I'd like to decompose my signal (func) to 2 components which are the sine wave and the wavelet.
After this, I'd like to decompose the wavelet to its gaussian and sine components.
I tried to do this with a continuous wavelet transform manually but I'm having troubles the recover the correct amplitudes of the wavelet and also with the logarithmic scale of the wavelet transform.
I also don't want to do this manually and would like to find an automatic way to extract those properties.
How is it possible to do this in another way?
I didn't try any sophisticated mathematics (which I don't know right now...)
I'm working with MATLAB.
So the code that creates the wavelet and the wavelet transform is:
clearvars;
close all;

j=sqrt(-1);
f_sys=50;
SR=1024;
cycles = 10;
Fs=SR*f_sys;

dt=1/(SR*f_sys);
eps=0.01; %half power bandwidth

miu2=0.15; % center move over t
f0=f_sys;

dF = 50;
fn = 0:dF:10e3;
wn = fn*2*pi;
t = 0:1/(SR*f_sys):(cycles/f_sys)-(1/(SR*f_sys));

for i=1:length(wn)
    phin = sqrt(1/0.0125)*(exp(-((t-miu2)/(eps)).^2));
    psin(i,:) = phin.*exp(j*wn(i).*(t-miu2));
end

sin2= 50.*exp(j*2*pi*f0*t);
sig=sin2;
wvlt = psin(8,:);
func = sig+wvlt;

subplot(3,1,1)
plot(real(sig))
title("sine wave");
subplot(3,1,2)
plot(real(wvlt))
title("wavelet");
subplot(3,1,3)
plot(real(func))
title("func = wavelet+sine");

plot_min_freq = 0;
plot_max_freq = 1000;

fb = cwtfilterbank('SignalLength',length(t),'SamplingFrequency',(SR*f_sys),'FrequencyLimits',[plot_min_freq plot_max_freq],'Wavelet','amor');
figure;
cwt(func,'FilterBank',fb);
[wt1,f]=cwt(func,'FilterBank',fb);

This is an example of a signal which I create but I'd like to do the same on a real signal.
The signal is composed of a sine wave at 50Hz and a morlet wavelet at 350Hz.
I thought about some first steps:

finding the 2D peaks in order to find the wavelet - I didn't succeed yet with the noise that I get from 50Hz and it is also possible that I'll get more than one wavelet in a specific time.
take a vector of the row where the peak of the wavelet is located - then I can see a kind of gaussian.
fit #2 to a gaussian.

Now I'm stuck and I don't believe that this is the correct way to solve it.
Thanks.


texts: Relationship between Fs (the Nyquist frequency), and the frequency used with a sine wave

texts:
Say you are given a vector t =[...] for the time to be used to generate a sine wave of a certain frequency say FHz. I have all the data for t so I can conclude what dt is (my dt = 1d-5). Now I am to do an FFT analysis and plot magnitude vs frequency so I expect to get my original amplitude and also get back the original phase
This what I do:

frequency = ...; %This in Hz

dt = t(2) - t(1); %which is equal to 1d-5

fs = 2*frequency; % should this be 1/dt? or something else?

Nfft = fs;

intensity = sin(2*pi*frequency*t);

FFT_intensity =fft(intensity,Nfft)/Nfft;

magIntensity = abs(FFT_intensity);

f=linspace(-fs/2, fs/2, Nfft);

plot(f,fftshift(magIntensity))

The results I get are sometimes correct and sometimes incorrect, it seems like when I chop the length of the t vector or change fs (say fs = 1/dt instead) I get the wrong amplitude ( and wrong phase, not shown here)... In general how to relate FHz, Fs and t?
I can actually include the numbers I am using if that makes any difference.
Thank you for any input.


texts:
Property:

The frequency of oscillation of discrete time sinusoids sequence increases as $\omega$ increases from $0$ to $\pi$. If $\omega$ is increased from  $\pi$ to $2\pi$ then frequency of oscillation decreases.

My question:

What is meant by "frequency of oscillation" ?

How can a frequency of a signal vary (increase) from $0$ to $\pi$ and (decrease) from $\pi$ to $2\pi$. Isn't it should be same from $0$ to $2\pi$ ?


texts:
How do you estimate the detection SNR for a radar transmitting an LFM (chirp) waveform? I'm looking for an SNR equation that is a function of chirp bandwidth and compressed (or uncompressed) pulsewidth.

Do you use the uncompressed pulsewidth or the compressed pulsed width? Do you add a term for pulse compression gain (proportional to the time bandwidth product = bandwidth*pulsewidth?) to the radar range equation?

According to Mahafza, "For a given set of radar parameters, and as long as  the  transmitted  pulse  remains  unchanged,  then  the  SNR  is  also  unchanged regardless of the signal bandwidth. More precisely, when pulse compression is used, the detection range is maintained while the range resolution is drastically improved by keeping the pulse width unchanged and by increasing the band-width."


texts: Signal Processing and Emotion Recognition

texts:
I have to write a paper on emotion recognition.
More specifically on how different parameters affect the algorithms effectiveness. I have already read a lot about this and done research using OPENSmlie, (extracting mostly MFCCs) and WEKAs clasificators to get the best results. I use EMO DB.

However what i need is more theoretical knowledge, and information about already published papers on that subject. Especially newer ones, say from the last couple of years. I want to know, how peaople approach this problem, and how good their latest results are.

I was wondering if maybe there were some experts on that here.
Could you recommend any good recources / links to such papers. Or an effective way of finding the latest ones?


texts: Is multirate signal processing involved in ADSL?

texts:
Telecom companies offer different data rates for home/office purposes depending on the money they pay.

How are the telecom people able to control the flow of data of different cables?
Say supposing a cable can be used till $20\textrm{ MBps}$. But he paid only for $5\textrm{ MBps}$ scheme, why is the user not able to steal the remaining bandwidth?
Is there a multi rate theory behind this?
Or it is a simple filtering?


texts: How to derive the results that averaging $N$ signals yields a $\sqrt{N}$-fold increase in signal-to-noise ratio?

texts:
If I have understood correctly, averaging $N$ noisy independent segments or signals increases the signal-to-noise $\sqrt{N}$-fold.

How does one derive this result?


texts: WiFi and the sampling rate of spectrum analyzer

texts:
I'm trying to understand the correlation between RF analyzesrs and WiFi.

Looking at the device's datasheet I see that it has a sample rate of (max.) 100msec.
[to my understanding, the device can measure the signal strength/power 10 times  (or less) per/each 1 second].

I understand that WiFi signals are being transmitted in bursts.
Hunderds/thousands transmissions per 1 second.
Each transmission last for a very short time, 3-5micro-second.
That is to say, that each 1 second there is a signal transmission-pause-signal transmission-pause and so on.

Since the analyzer measure only 5-10 times each second, is it a theoretical possibilty that the device can miss [shall not measure] most of the signals?


texts:
I'm working on adding MFCCs to my kNN classifier.
Originally my kNN classifier was using the zero-crossing-rate, centroid, bandwidth, skew, kurtosis.  This was about 80% accurate for drum tracks using identical instrumental to the data it was trained with (there were at least 10+ training tracks for each part of my musical piece, ie. bass kick, snare, hihat, one of each were used in a 10 second clip repeating each 3 times).

Now I added 12 MFCCs (dropping the 0th) to this list of 5 features, and my kNN classifier has broke.

What am I doing wrong?  Below is the part of my code, I have edited out the repetitive parts for extracting features, they are identical to the first one.

snareDirectory = ['/Snares/'];
snareFileList = getFileNames(snareDirectory ,'wav');

kickDirectory = ['/Kicks/'];
kickFileList = getFileNames(kickDirectory ,'wav');

HiHatDirectory = ['/Hihats/'];
HiHatFileList = getFileNames(HiHatDirectory ,'wav');

currentFrameIndex=1;
for i=1:size(snareFileList,2)
    [x,fs]=wavread([snareDirectory snareFileList{i}]);
    frameSize=.1*fs;
    currentFrame=x(1:frameSize);
    featuresSnare(i,1)=zcr(currentFrame);
    [centroid,bandwidth,skew,kurtosis]=spectralMoments(currentFrame,fs,8192);
    featuresSnare(i,2:5)=[centroid,bandwidth,skew,kurtosis];

    [mfceps] = mfcc(currentFrame ,fs)' ; %note the transpose operator!
    delta_mfceps = mfceps - [zeros(1,size(mfceps,2)); mfceps(1:end-1,:)]; %first delta
    % Calculate the mean and std of the MFCCs, MFCC-deltas.
    MFCC_mean(currentFrameIndex,:) = mean(mfceps);
    MFCC_std(currentFrameIndex,:) = std(mfceps);
    MFCC_delta_mean (currentFrameIndex,:)= mean(delta_mfceps);
    featuresSnare(i,6:17)=MFCC_mean(currentFrameIndex,2:13);
end

for i=1:size(kickFileList,2)
    {...}
end

for i=1:size(HiHatFileList,2)
   {...}
end

[trainingFeatures(:,1:5),mf,sf]=scale([featuresSnare(:,1:5); featuresKick(:,1:5); featuresHiHat(:,1:5)]);
[trainingFeatures(:,6:17),mf,sf]=scale([featuresSnare(:,6:17); featuresKick(:,6:17); featuresHiHat(:,6:17)]);
labels=[[ones(size(featuresSnare,1),1) zeros(size(featuresSnare,1),1), zeros(size(featuresSnare,1),1)]; [zeros(size(featuresKick,1),1) ones(size(featuresKick,1),1) zeros(size(featuresKick,1),1)]; [zeros(size(featuresHiHat,1),1) zeros(size(featuresHiHat,1),1) ones(size(featuresHiHat,1),1)]];
model_snare = knn(17,3,1,trainingFeatures,labels);

[y,fs]=wavread('together.wav');
framesize=.1*round(fs); %.1 second
features=zeros(floor(size(y,1)/framesize),5);

for i=1:floor(size(y,1)/framesize);
    currentFrame=y(i*framesize:i*framesize+frameSize-1);
    features(i,1)=zcr(currentFrame);
    [centroid,bandwidth,skew,kurtosis]=spectralMoments(currentFrame,fs,8192);
    features(i,2:5)=[centroid,bandwidth,skew,kurtosis];

    currentFrameIndex=i;
    [mfceps] = mfcc(currentFrame ,fs)' ; %note the transpose operator!
    delta_mfceps = mfceps - [zeros(1,size(mfceps,2)); mfceps(1:end-1,:)]; %first delta
    % Calculate the mean and std of the MFCCs, MFCC-deltas.
    MFCC_mean(currentFrameIndex,:) = mean(mfceps);
    MFCC_std(currentFrameIndex,:) = std(mfceps);
    MFCC_delta_mean (currentFrameIndex,:)= mean(delta_mfceps);
    MFCC_delta_std(currentFrameIndex,:)= std(delta_mfceps);
    features(i,6:17)=MFCC_mean(currentFrameIndex,2:13);
end

soundsc(y,fs)
featuresScaled(1:5) = rescale(features(1:5),mf,sf);
featuresScaled(6:17)=rescale(features(6:17),mf,sf);

[voting,model_output]=knnfwd(model_snare , featuresScaled );
output = zeros(size(model_output,1),3);
output(find(model_output==1),1)=1;
output(find(model_output==2),2)=1;
output(find(model_output==3),3)=1;

I have read this post
https://stackoverflow.com/questions/16546524/how-to-use-mfcc-vectors-for-classifying-a-single-audio-file
which is close to what I am trying to figure out but in my eyes I am already doing what the answer is saying.

Any help is greatly appreciated, I have spent several days trying to figure this out and have only resulted in creating a model that is 50% less accurate than without the MFCCs.

Thanks in advance!

-scarlson


texts: What is QEEG technically, and how is it computed?

texts:
I've been working with EEG signals for a while. I wonder how to compute QEEG (or Brain Map) representation. I think it is extracted from raw EEG and I want to know what is technical method for QEEg extraction? Is it simply a Fourier transform or something complicated is utilized?

In a nut shell:
 - What is qEEG exactly and how to compute it technically to come up with something like in figure above?

I feel it should be related to short time-Fourier transform.


texts:
I have a floating-point digital signal processing system that operates at a fixed sample rate of $f_s = 32768$ samples per second implemented using an x86-64 processor.  Assuming that the DSP system is synchronously locked to whatever matters, what is the best way to implement a digital oscillator at some frequency $f$?

Specifically, I want to generate the signal:
$$y(t) = \sin(2\pi f t)$$
where  $t=n/f_s$ for sample number $n$.

One idea is to keep track of a vector $(x,y)$ which we rotate by an angle $\Delta\phi = 2\pi f/f_s$ on each clock cycle.

As a Matlab pseudocode implementation (the real implementation is in C):

%% Initialization code

f_s = 32768;             % sample rate [Hz]
f = 19.875;              % some constant frequency [Hz]

v = [1 0];               % initial condition
d_phi = 2*pi * f / f_s;  % change in angle per clock cycle

% initialize the rotation matrix (only once):
R = [cos(d_phi), -sin(d_phi) ; ...
     sin(d_phi),  cos(d_phi)]

Then, on each clock cycle, we rotate the vector around a little bit:

%% in-loop code

while (forever),
  v = R*v;        % rotate the vector by d_phi
  y = v(1);       % this is the sine wave we're generating
  output(y);
end

This allows the oscillator to be computed with only 4 multiplications per cycle.  However, I'd worry about phase error and amplitude stability.  (In simple tests I was surprised that the amplitude did not die or explode immediately--maybe the sincos instruction guarantees $\sin^2+\cos^2=1$?).

What is the right way to do this?


texts: We are given a set of grayscale image patches obtained from images after edge detection. Each patch is 10x10 pixels with intensity varying between 0 and 255 for each pixel. This set may contain a very few (maybe one or two) or a large number of patches. We want to represent the set by a single 10x10 pixel patch or “feature” such that the feature captures the invariance or common properties of the set. Intuitively speaking, the feature should be the intersection of the edges in the image patches in the set. For example, if all the patches contain a vertical edge at a fixed position and also contain other edges in different orientations and positions, then our feature should contain only the vertical edge at the fixed position and no other edge. Unfortunately, the operation of intersection is not well-defined for images. How can we approach to compute such a feature for a set of image patches?


texts: Implementing a bandpass filter in R - this code's logic confuses me

texts:
I am reading this "how to" in R: http://www.inside-r.org/howto/time-series-analysis-and-order-prediction-r

I am confused about the tutorial's implemented band pass filter. In particular, here's the code I'm wondering about. The goal is to extrapolate a time series based only on certain ranges of frequencies determined by peakind.

subsignals <- lapply(c(peakind$freqindex, midindex+1), function(x){
    upperind <- x
    fsub <- f
    notnullind <- ((fsub$freqindex >= lowerind
                    & fsub$freqindex < upperind)
                   |
                    (fsub$freqindex >  (lindex - upperind + 2)
                     & fsub$freqindex <= (lindex - lowerind + 2)))
    fsub[!notnullind,"coef"] <- 0
    lowerind <<- upperind
    Re(fft(fsub$coef, inverse=TRUE)/length(fsub$coef))
})

My question is, why this condition:

 (fsub$freqindex >  (lindex - upperind + 2)
                         & fsub$freqindex <= (lindex - lowerind + 2))

This is basically saying to include (for computation of inverse FFT) those frequencies that are within the narrow range of the top frequency minus upper band frequency +2 and top frequency minus lower band frequency + 2.
So basically just some tale at the top of the frequency. This doesn't make any sense to me. Can someone tell me where I should read up on this or what an intuitive hand-wavey sort of justification is?


texts: Checking Parseval's Theorem for Gaussian Signal by Using Scipy

texts:
I'm trying to check Parseval's theorm for Gaussian signal. It's well known that fourier transform of $\exp(-t^2)$ is $\sqrt{\pi}\exp(-\pi^2 k^2)$. So I implement it by using quad and simps. I think that the one is continuous integral, the other is integration which uses sample.

Code

from scipy import integrate
import numpy as np

def f(x):
    return (np.exp(-x**2))**2
def F(k):
    return ((np.pi**0.5)*np.exp((-np.pi**2)*(k**2)))**2

a=integrate.quad(f,-np.inf,np.inf)
print(a)
#>>>(1.2533141373155017, 4.4674503165883495e-09)
b=integrate.quad(F,-np.inf,np.inf)
print(b)
#>>>(1.2533141373155066, 1.592743224014847e-08)

a=b so above code follows Parseval's Theorm. The problem is on following code.

import scipy.fftpack as fft

N=50000
t  = np.linspace(-1000,1000, N)
h=np.exp(-t**2)
H=2*np.abs(fft.fftshift(fft.fft(h)/N))
freq=fft.fftshift(fft.fftfreq(H.shape[0],t[1]-t[0]))

S_h=integrate.simps(h**2,t)
print(S_h)
>>>1.25331413732

S_H=integrate.simps(H**2,freq)
print(S_H)
>>>1.25326400525e-06

S_h is not equals to S_H. What is the problem? How can I fix it??


texts:
How to find poles of transfer function by looking at the step response?

Given a step response graph like such:

How would I find the sketch for its poles on the complex plane? The only thing I can gather are the rise time, settling time, overshoot, etc. But how do I derive the poles from that?


texts: ssbmod single sideband modulation artifacts

texts:
clc; clear all; close all;
[c,fs]=audioread('60hz.wav')

r=ssbmod(c,160, fs)
r=r/max(max(abs(r)))

audiowrite('artifacts.wav',r,fs,'BitsPerSample',24)

I am trying to understand what causes the artifacts such as the one example in the picture below...This happens when the audio sample I am processing is 142157 samples in size, but if I use multiples of 44100 sample chunks of 60 hz sinewav I won't get that upper side band that is present in the picture below, the obvious Hilbert Transform is in play but what aspect of it is causing the artifact is what I am trying to find out...
this is the code for ssbmod:

y = x .* cos(2 * pi * Fc * t + ini_phase) + ...
        imag(hilbert(x)) .* sin(2 * pi * Fc * t + ini_phase);


texts:
I hope this question is not too simple, I just started learning digital image processing. The 1D binomial filter of size 2 is defined by $B_2 = \frac{1}{4}\begin{bmatrix}1 & 2 & 1\end{bmatrix}$. The DFT should be $\frac{1}{2} + \frac{1}{2}\cos(\omega)$, but I have some trouble deriving this solution. I found it here.

My calculations:

$$X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N}nk} = \frac{1}{4} + \frac{2}{4}e^{-\frac{2\pi i}{3}1k} + \frac{1}{4}e^{-\frac{2\pi i}{3}2k}$$ Then I took the real part $Re(X_k)$, and found after simplifying $\cos^2(\frac{\pi k}{3})\cos(\frac{2\pi k}{3})$, which is incorrect. Only if I use two times the second element in $B_2$, I get the same result.

I'm also trying to understand the reason for applying the DFT in this case. My line of thought is: By taking the Laplace/z-transform of the impulse response, one obtains the system's transfer function. Because we are only interested in the steady-state response, we can ignore the imaginary part and calculate instead the fourier transform. Since the input $B_2$ is real, the real part of the DFT is enough. Is this correct?


texts: I have a pretty good understanding of what chroma is when represented against time, but I don't understand what it means when one represents chroma against beats


texts:
I read a lot about Wiener Filters (focusing on discrete time case). I understand the math, but I am quite disconnected from the real life assumptions behind using such a filter. Unfortunately textbooks are just laying out the formulas without putting forward much intuition in it.

The task is to extract x(t) given y(t), where y(t) = x(t) + n(t). Great but first we need to build an appropriate filter. At this point:

Are we assuming that we have access to some training samples of x(t) and/or n(t) in addition to y(t)?
Or are we just assuming that we have access to statistical properties(auto correlation, cross correlation) of x(t) and/or n(t) in addition to y(t)? If so since x(t),n(t) are not explicitly available, are these statistics computed from other signals available to us with similar statistical properties?
Or is it a combination of both, such as having access to some training samples of x(t), but having only access to second order statistics of n(t) (or vice versa)?
Is the noise n(t) always assumed to be uncorrelated to x(t) in Wiener filtering, or are the textbooks going for the simplified version? What is the default scenario?
Why are all the examples coming up with all-zero filters but not using pole-zero filters? Is it just a matter of simplified textbook treatment or is the problem definition not amenable to using transfer function modelling approach? Yes the first approach is computationally easier, but these days we have significant computing power and I don't think it is a good justification anymore, if that's the only reason.

Looking forward to hear from practitioners.


texts:
We have several audio tracks, all of which are different versions of the same track under one type of distortion: the speed of the track is increased by a constant factor K (hence frequences are multiplied by K, timesteps are divided by K), K is relatively small (0.8< K < 1.2).

I want to know the best way to get a good spectrogram out of these tracks.
Before defining what I mean by 'good', let me state my final goal: I want to retrieve points, à la Shazam, based on the maximum amplitudes (ie amplitudes > a certain threshold). A good threshold is to be defined too, but I was thinking of taking the average local RMS of the signal and its standard deviation and make the treshold = mean(localRMS) + (constant factor)*std(localRMS). Then I'd like to apply computer vision algorithms to the resulting image to measure similarity.

A good spectrogram is then defined as such:

the frequencies are logarithmically spaced the points of interest
for one track are a translation along the frequency axis of the points of interest of
another track.

Do you know what type of spectrogram would be the best? (I'm thinking of using either a logarithmically spaced fourier transform or a constant Q Fourier transform) Would my threshold work or are there better ones? Any remark on the whole design would be truly appreciated.


texts: Application of notch filter in radio receivers

texts:
So according what I have read so far notch filter (or stop-band filter) can be used to get rid of the narrowband noise when it is present in the frequency spectrum of the received signal.

For example stray sine waves may impinge on the receiver as we listen to radio.

My question is, how would the notch filter detect the particular noise frequency in first place so that it could remove it?


texts: Can every type of linear filter be modelled by a convolution?

texts:
I have an input time series going through a filter that creates another time series as output.

If I assume in first approximation that my filter is linear, does it necessarily mean that I can model this filter by a convolution?

In other words is it only certain types of linear filters that can be modelled by a convolution? (LT1?...).

I have never seen any clear statement on this in any publication. It is as obscure and vague as definitions of time invariance, ergodicity, causality etc in filtering textbooks.

My question comes from wondering if there is any point trying to attempt an already risky deconvolution if I am not even sure that my linear filter is necessarily akin to a convolution.

Why classes of linear filters are not clearly explained in publications (that I know of) is a mystery to me, as it should perhaps be the first question that arises when trying to model real world phenomenons.


texts: I have a system with limited memory. So basically the question is this: If I have (say) 24KB of data space and 4096 data samples of single precision floats (16KB) can I do an FFT on the data? More generally, what is the minimum data space required to do an FFT?


texts: Which filters are the most commonly used and easier to implement for ECG noise removal?

texts:
ECG signal is affected by interferences such as the Powerline Interference, Baseline Wandering and Muscle Noise.

There are many filters available which can reasonably cancel these types of noise.

I've found that Butterworth, Chebyshev, Wiener and Kaiser filters, along with wavelets could be used for this purpose.

Now, my main questions are:

Which filters are the most commonly used and easier to implement
(both from the above list and from other filters you have in mind)?
Can you provide me some pseudo-code samples on how to implement some simple filters, maybe Butterworth or Wiener?
In which references should I look at concerning ECG digital
filtering, I mean which provide an overview of the various filters instead of foccusing on a specific one?

Thanks in advance


texts: What is the difference between t=n*Ts sampling vs impulse train sampling?

texts:
I know that if I sampling with impulse train so I get in the frequency plane
X(f)*h(f)
(when x(f) is my signal, * means convolution and h(f) is fourier transform of impulse train).

what the difference between this way to x(t=n*Ts) sampling way?

If I can get explenation in the "frequency plane" it will be great.


texts: Light Detection and Ranging (LIDAR) Range Bias Compensation with Light Dependent Resistor (LDR)

texts:
I have a distance sensor (VL53L0X), which measures the distance to an object. It is placed outside. I currently do 100 measurements in the best available range mode, then I take the median and insert it to a database.

The current distance is about 1m.

This is how it looks like during a period af about 12 hours. As can be seen there is a light dependent error in the data.

I also have a LDR (light dependent resistor) output (log scale on y-axis), which shows this more. At "peak" time the error is about 15cm.

Green: distance in meters, yellow, Voltage from LDR (light dependent resistor). The darker it gets, the higher the voltage of the LDR. At night it's about 2.5V.

Now, the distance hasn't changed in those 12 hours. My question is: how can I compensate this error? I am looking for a term on how this procedure is called.

My first guess was to calculate the difference to 1m for all distance points, then maybe fit the LDR output to that error points.

Data: https://pastebin.com/2fXcSUQB

Column 1: distance in meters, Column 2: voltage from LDR


texts: What is the overall capacity of two cascaded BSC channel?

texts: Is it true:$C = C_1 * C_2$ where $C_1$ and $C_2$ are the two cascaded BSC which makes $C$.


texts:
I want to detect filopodia (red arrow indicates filopodia in a vessel image).

My strategy is that if $I$ is the original and $J$ is the another image after opening operation. Then, $I-J$ is the filopodia. (See Figure 6 in this paper).

To do this, I want to enhance the very thin lines. How to do it? I also have information of front line like this:

Or, do you have any idea for detecting filopodia?


texts: Can't make sense from VHDL butterworth filter implementation

texts:
I am trying to understand the VHDL implementation of the 3rd order Butterworth filter published on opencores. Supposedly the filter implementation gets away with only division and multiplications by 2 which is cheap when implemented in hardware (only shifting the bits left/right).

I've translated the VHDL to MATLAB to simulate a step response.

x = int16([zeros(1,500) 512*ones(1,1000)]); % step-input

a = zeros(1,3); % accumulator registers
w = zeros(1,5); % analog filter state variables

s = 4; % scaling parameter to adjust cutoff frequency

for t=1:length(x)
  in = x(t);

  a1 = a(1) + w(1) - w(3);
  a2 = a(2) + w(2) - w(4);
  a3 = a(3) + w(3) - w(5);

  w1 = in - w(2);
  w2 = a(1)/power(2,s);
  w3 = a(2)/power(2,s+1);
  w4 = a(3)/power(2,s);
  w5 = w(4);

  % update registers for next cycle
  a = [a1 a2 a3];
  w = [w1 w2 w3 w4 w5];
  y(t) = 2*w(5);
end

The step response looks quite promising; dotted line is the step input, solid red, orange and blue are the filter response for s=3, 4 and 5 respectively.

However, the implementation really does not make much sense to me. Can anyone enlighten me? What are the coefficients of the filter? Which form does the implementation follow?


texts: Bode plot of discrete-time transfer function $H(z)$

texts:
$H(z)$ is the transfer function of a biquad filter as described here.

I would like to plot the Bode plot of the magnitude response of $H(z)$.

Scipy has a bode method (scipy.signal.bode) for continuous-time transfer functions. Is there an option I didn't see for discrete-time functions?

I found in this nice file formula 18. I tried plotting it with gnuplot directly, but that didn't look correct.

set samples 100000, 100000
set logscale x

b0= 0.2514
b1= 0.5028
b2= 0.2514
a1=-0.1712
a2= 0.1768

H(x) = sqrt((b0**2 + b1**2 + b2**2 + 2*(b0*b1+b1*b2)*cos(x) + 2*b0*b2*cos(2*x))/(1 + a1**2 + a2**2 + 2*(a1+a1*a2)*cos(x) + 2*a2*cos(2*x)))

plot [1:22050] H(x)

So, how do I do this? Is  it even possible to Bode plot the magnitude response of a discrete-time transfer function?


texts:
I have implemented a pre-emphasis filter by the following (pseudo) code:

a = 0.5;
s1[0] = s[0];

for (n = 1; n < N; n++) {
  s1[n] = (s[n-1] * a) + s[n]
}

The problem I am finding is that due to summing, there is clipping happening.  I thought the solution would be to find the difference between the max value in s1[n] and 1.0 and then subtract s1[n] by that amount...

However, that results in my audio turning into complete garbage, which I am confused as to why........  Subtracting a constant amount from an entire signal should simply lower it's amplitude by that amount, correct?


texts: What does it mean when you zero pad data and why does one do it?

texts: I am struggling to understand what effect zero padding data would have on EMG data. How does it improve the signals? Any help in simple terms that's not too mathematical would be greatly appreciated!


texts: Are Convolution and Deconvolution Kernels the Same?

texts:
I need to clarify this and rather confused by this.

Lets say:

  $x = h * g$

  $x$ - measured data

  $h$ - raw data

  $g$ - instrument response function (convolution kernel)

So lets say I have $x$ and $h$. Can I say if I were to deconvolve $x$ with $h$ I will be able to obtain $g$?

Or if I have $x$ and $g$, deconvolving $x$ with $g$, I will be able to obtain $h$. So till here, is my understanding correct?

So can I say $g$ is both convolution and deconvolution kernel? Is this correct?


texts: Envelope behavior dIfference between QPSK, OQPSK and $\pi/4$-QPSK

texts:
I am trying to analyze whats the major difference between these three variants of 4PSK.
I know that
QPSK -  signal points are at 0,90,180,270
OQPSK - The bit period of Quadrature component is shifted by one bit period
pi/4-QPSK -  Signal points are at 45,135,225,315

Some website says that its to avoid abrupt phase transitions through zero.. I cannot understand and figure it out the actual reason behind it.. Can anyone please point it to me


texts: I somehow "proved" that given any LTI system, its transfer function has to be constant. What am I missing?

texts:
A transfer function is defined as the Laplace transform of the ratio of output to input.

Also, every LTI system has an eigenfunction. Given such eigenfunction as an input, the ratio of the output to input is a constant, which means the transfer function is constant.


texts: Why does the Y component in YCbCr represent luma/intensity?

texts:
As far as I know, the RGB to YCbCr conversion is just a linear transformation, yet the Y channel successfully captures the intensity value of each pixel. How?

In other words, where do the following magic numbers come from?

$$
\begin{align}
E'_Y &= 0.299E'_R + 0.587E'_G + 0.114E'_B \\
E'_R - E'_Y &= 0.701E'_R - 0.587E'_G - 0.114E'_B \\
E'_B - E'_Y &= -0.299E'_R - 0.587E'_G + 0.886E'_B
\end{align}
$$


texts: advantage of using impulse response function for LTI systems?

texts:
i have a discrete-time LTI system $L$ that takes input signal $x[n]$ and gives the output signal $y[n]$. since $L$ is LTI, $y[n]$ can be derived as a sum of shifted and scaled impulse responses of $L$. but why is that a useful fact and can you show a practical example where this saves us work?

this is my attempted answer and i would like to know if it's correct: suppose $L$ is an accumulator,

$y[n] = \sum_{k=-\infty}^{n}x[n]$

then by LTI properties, this can be rewritten as a sum of shifted and scaled impulse responses:

$y[n] = \sum_{k=-\infty}{\infty}h(n-k)$

where $h(n-k)$ is the impulse response function.

what does this gain us, intuitively? it seems that this means that if we could only measure the response of $y(n)$ to a single (unit) impulse, we could derive through shifting, scaling and summation $L$'s output to any arbitrary input signal $x[n]$. so we wouldn't have to measure $L$'s output to a more complicated signal, a single pulse would be enough. is this right?

if so how would it work in this example? suppose the input signal $x_0[n]$ was:

$x_0[0] = 1, x_0[k] = 0, \forall k \neq n$

then clearly $L[x_0[n]] = 1$. how can we use this to derive the response to a new input signal $x'[n]$:

$x'[0] = 1, x'[2] = 1, x'[3] = 2$

the result should be 4. how do you formally derive this answer using only the system response to $x_0$?


texts:
Location of poles are given in unity circle and 4 zeros are given at origin for every plot,How to check which plot shows which filters?

For the Option (C) i am doing like this

I can write the Tranfer function from the plot like this

$H(z)=\frac{z^4}{Z^4-1}$

Now at z=1,$H(z)=\infty$

at z=j at mid frequecny,its its also coming $H(z)=\infty$

at z=-1,$H(z)=\infty$

How can i conclude that which filters it should be ?

What should be the general procedure to do this?


texts: Relation between Fourier Series & Fourier transform

texts:
So i was just revising some basic DSP concepts. Just wanted to verify this fact.

Fourier series represents a periodic signal $\hat{x}(t)$ with period P as a countably infinite sum of sinusoids of frequency $0$, $\frac{1}{P},\frac{2}{P},\frac{3}{P}...$. This converges to the signal in the interval, $-\frac{P}{2} < t < +\frac{P}{2}$, and if the time domain signal is periodic, then over the whole time domain.

Fourier Transform is sorta like a limit of the fourier series where P goes to $\infty$.

So i know that the fourier transform of $\operatorname{rect}(t)$ is $\operatorname{sinc}(f)$ ( ignoring the scaling factors ) . And that the fourier series of a $\operatorname{rect}()$ is given by http://mathworld.wolfram.com/FourierSeriesSquareWave.html ( which is also a $\operatorname{sinc}()$ in the frequency domain ) .

I just wanted to confirm the following

If I sample the $\operatorname{sinc}()$ i obtain from the fourier transform of a $\operatorname{rect}()$, and use those values to reconstruct a fourier series, will i end up getting a square wave ?


texts: Ground plane estimation from a single image (was: confusion about homography and calibration camera parameters)

texts:
I'm new to this multiple view geometry issues and I'm a little bit lost, and some concept are not clear at all.

I have a calibrated camera ( = I have distortion coefficents, field of view, and with those I can find what is called intrinsic matrix).

With this camera I'd like to make an application that can estimate the ground plane from a picture. I have those constrains:

In the picture there is no chessboard (or other kind of known pattern used in computer vision for calibrating cameras)
I can manually define some coplanar points (like points that I know are the projection of a square) or some

So I need something that take the 2D coordinate of a trapezoid in the picture, knowing that in the real world it is a unit-square on the ground plane in the 3D real world coordinates and give me the transformation.

From what I have understand this is procedure is "computation of extrinsic parameter given given intrinsic parameters, a few 3D points, and their projections".

Now, my questions:

Is it possible?
Are 4 points enough?
This is not related with homography, right? (because homography is not related with 3d points but only with 2d).. but how can homography help me? (I can find the homography that transform the trapezoid in a square)
Is there an alternative way to do it (like with triangles and normals)?

For the sake of curiosity I'm using openCV for the implementation.


texts: Given a downscaled image, produce Gaussian blur of the original image

texts:
Let's say I a perform convolution of an image with a Gaussian kernel of some size.

Is it possible to produce the exact same result from a downscaled version of the image? If so, what's its minimal size?

I'm not looking for an exact answer, some directions where I should start looking are enough. (I have some basic DSP and image processing knowledge but not much beyond that.)

Edit: A closely related question is - does Gaussian blur lose information, i.e. is it impossible to reconstruct the original from the blurred image?


texts: evaluating individual terms of 2D DFT (Goertzel?)

texts: I'm looking for a fast way to do a 2D discrete Fourier transform of an image at many arbitrary frequencies. I know the Goertzel algorithm works for 1D, but is it possible to generalize it in 2D? Or any other useful algorithms? Two solutions that I know, but which I'm afraid might be too slow, are directly using the DFT equation, and zero padding and doing an FFT. Alternately, the frequencies might not need to be completely arbitrary but could be sub-arrays of the array of frequencies given by an FFT with high frequency resolution.


texts: Is there an algorithm for second order filtering of data where the filter frequency is not known a priori?

texts:
It is well known that a moving average algorithm done in the time domain is equivalent to a filter with frequency response $\mathrm{sinc}(\omega\tau)$ where $\tau$ is the averaging time. (see this related answer)

This has the following beneficial property: you are streaming a time series of data ${x_n}$, and the average at any point ($a_n$) is just:
$$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$

Thus you may apply the above recursive algorithm for an arbitrary amount of time ($\tau$), and when you stop, the value you have is filtered by $\mathrm{sinc}(\omega\tau)$, and has a correspondingly reduced variance. Now the $\mathrm{sinc}$ function is a first order low pass, modulated by a $\sin$ envelope. So in effect you have done a first order low pass where the characteristic low pass time constant $\tau$ is equal to the length of the data stream, and $\tau$ was not necessarily known before you started.

My question is: is there some analogous procedure which allows for an (approximate) second order low pass where the time constant is not known a priori?

A possibility is to "average the averages" but that requires keeping all the averages in memory. Is there some law preventing such a procedure with small memory requirements?


texts: How to see the amplitude of each frequency?

texts:
do you know software that exports spectrograms to images of 500x22000 pixels?
all the spectrograms that I know show diffused frequencies and it does not look very good
I just want to see the amplitude of each frequency

each vertical line is one second of the sound
each horizontal line is a frequency

the entire audible spectrum enters the height of the image

the amplitude of each frequency is the brightness of the pixel


texts: Phase rotation without Fourier Transform

texts:
Is there any known math equation or other method to perform phase rotation on signal without decomposing it to frequency domain?
In frequency domain it is obvious. You just need to perform phase shift for each frequency.
But decomposing signal from time domain to freq domain, then perform phase shift on each frequency, then again composing it to time domain it is very computationally demanding for frequency-rich signal.
That's why I wonder if I can achieve it without decomposing signal to frequency domain?
To be clear I want to perform on the signal something like that:

But I don't know the frequencies in the signal. I just have time domain signal. Is that possible?


texts: How to take into account phase differences when doing Fourier Transform?

texts:
I apologize if the title seems very vague, I don't really know how to put it a different way.
From what I understand, the Fourier transform is essentially just doing a dot product between multiple sine waves of different frequencies with the original signal. When the "chosen" sine wave has the same frequency as a component in the signal, it will be in phase and produce peaks. If the "chosen" wave isn't in the original signal, then it will either be cancelled by the original signal, or produce smaller peaks.
Let's say that there's a 100Hz sine wave component in a signal with a phase of $\frac{\pi}{2}$. If I were to take the dot product with a pure sine wave of 100Hz with no phase, it would appear as if two signals are present, one at 95Hz and one at 105Hz:

Is there a way to negate the effects of this? I was thinking of taking the maximum of pure cosine and sine wave in order to take into account phase differences, but I get much more noise around the frequency component:

Note: I'm not using any windowing function (Hanning, etc.) I can post my Python code if that will help.
EDIT: I tried using the formula as per Marcus' comment, and I'm getting weird results:

def do_fft(y, rate, f):
    N = len(y)
    summ = 0

    for n in range(N):
        summ += y[n] * np.exp((0 - 1j) * 2 * np.pi * n * (f/rate))
    return float(summ.real)
    #return float(abs(summ))

def selective_fft(y, rate, frequencies):
    ar = []
    for f in frequencies:
        ar.append(do_fft(y, rate, f))
    return ar

"frequencies" is an array of frequencies that I want to "test" for.
EDIT: Here is how I'm doing the pseudo-fourier transform (originally):
t = np.linspace(0, 0.1, 4410)

ar = np.arange(0, 500, 1).reshape((-1,1)) #The frequencies to test

m = np.tile(ar, (1, t.shape[0])) #create matrix for frequencies
time_matrix = np.tile(t, (ar.shape[0], 1)) #create matrix for t vector

a = 100 #frequency components
b = 200
c = 300

y = np.sin(2 * np.pi * a * t + np.pi/2) + np.sin(2 * np.pi * b * t) + np.sin(2 * np.pi * c * t)

yc = y.copy() #copy signal

y = y.reshape((-1,1)).T #Add column and transpose to get one row and 4410 columns

y_stacked = np.tile(y, (ar.shape[0], 1)) #stack the function based on how many frequencies to test (500 from before)

sub = np.sin(2 * np.pi * np.multiply(time_matrix, m))
sub2 = np.cos(2 * np.pi * np.multiply(time_matrix, m)) #use cosine

dotted= np.abs(np.sum(np.multiply(y_stacked, sub), axis = 1).reshape((-1,1))) #basically the dot product
dotted2 = np.abs(np.sum(np.multiply(y_stacked, sub2), axis = 1).reshape((-1,1)))


texts: Phase-vocoder normalization for different speed factors

texts:
Question in short: How can I normalize audio from a block-based phase-vocoder system without knowing the global minimum and maximum during runtime?
Longer version:
Hey! I am very new to DSP and need some help with figuring out an issue that I have (or understanding the math behind it). I am using the RubberBand audio time-stretching/pitch-shifting C++ library and stream audio to and from it. But, it results in differently-normalized audio.
What I mean with this, is that I stream sample-blocks that are normalized within [-1, 1]. However, it returns audio within some range [x, y] where x and y are seemingly quite random, except for a speed-factor of 100%, where it stays perfectly normalized. My question is: is there some way I can know x and y, without needing to stream all my audio beforehand (for real-time usage)? A graph of the minimum and maximum of the audio is given below, so you can see what I mean:

I hope to find out what the mathematical source of this is. I have read online and believed that phase vocoders are supposed to keep audio normalized as long as the input was normalized?
Thanks in advance.


texts: Potential issues arising from too stable discretization

texts:
When numerically simulating a system, usually some kind of discretization is necessary, obtained by some kind of z-transform, such as, for instance, the bilinear transform $s\mapsto \frac{2}{\triangle t}\frac{z-1}{z+1}$, which is kind of the same as using the Trapezoidal rule to approximate integration.
It also has the very nice property that it preserves BIBO-(in)stability from the s into the z-domain. However, when it is implicitly used as a differentiator, there is a very real possibility of unwanted numerical oscillations.

Now take, for example, the Euler backwards method, corresponding to the substitution $s\mapsto \frac{z-1}{\triangle t}$. This method reduces numerical oscillations quite effectively, which would (ignoring for now the slower convergence) seem to be quite a big advantage. However, if we do look at the BIBO-stability areas before and after the z-transformation, this method leads to an effectively increased stability region, i.e. there are some transfer functions that are unstable in the s-, but stable in the z-domain.

1) That the latter method is able to dampen numerical oscillations seems intuitively connected to a greater "numerical stability". However, as I understand, not all systems where numerical oscillations occur after using the bilinear transform are inherently unstable. A simple RL-circuit with transfer function $\frac{U}{I}=\frac{1}{R+Ls}$ will a current source attached will, for instance, show numerical oscillations when the source is suddenly switched on or off. So are these phenomena - increased stability domain after the transformation, and damping of numerical oscillations - connected?

2) In literature, I do find a lot about choosing a suitable z-transformation such that the resulting discrete system is stable. However, there was no regard given to possible cases where the original system might have been unstable. To me, this seems like something of an oversight - when simulating an unstable system, I would expect the discrete system to also be unstable. Because, well, the output signals for stable and unstable systems would possibly differ a lot. Have I misunderstood something there? Or is there a reason why one can safely assume the systems in question (namely, electrical power systems) to be stable?

I do know that many algorithms (say, eg, EMTP and all related) usea mix of the above, but while this decreases the difference in stability regions, it is still there.

Edit:
By numerical oscillations using the Trapezoidal rule I mean phenomena such as those outlined in http://www.ece.uidaho.edu/ee/power/ECE524/spring14/Lectures/L39/numerosc.pdf , or similar (found using google, there seem to be different texts about the same).


texts: Can I combine an IIR (an EWMA) with a maximum?

texts:
I'd like to have a exponentially weighted moving average (EWMA) which is raised to the value of the input signal whenever this is higher than the filter output.  This taking of the maximum shall take place during the iteration, not afterwards (which would be simple).

So, the algorithm should look like this:

o = 0  # or some arbitrary initial value
for i in input:
    o = (o * 99 + i) / 100
    if i > o:
      o = i
    print o

So effectively, my output shall fall slowly after a quickly falling edge in the input without rising quickly before a rising edge.

In fact, any other algorithm which just copies the input except for the places where it has a rapidly falling edge (there it should play parachute and let the output fall slower) would do fine.

The only restriction I have is that this needs to be implemented in the usual Python libraries (numpy, scipy, pandas); iterating myself will not be fast enough.

I tried applying a rolling maximum, followed by a rolling mean with a Gaussian window of the same width (using pandas.Series(arr).rolling()); this gives me nice smoothing of the fall after a falling edge; unfortunately it also gives me an equal smoothing before the rising edges (which I don't want).

If there was a way to apply an asymmetric Gaussian window (e. g. rising slowly and falling quickly) that would be a solution I guess.  But I haven't found a way to achieve this yet.


texts: PSOLA how to calculate synthesized Marks

texts:
I have successfully extracts all Pitch Markers from waveforms, I used autocorrelation to find my pitch periods, The markers can be observed in the figure below:

I want change my signal (Pitch or Speed) using PSOLA, I need find the synthesized Marks, for example if I choose change my Pitch one octave UP(0.5), for it I need find the new markers to synthesize my signal, you can see in this figure below one example.

How can I calculate the new marks to be synthesized (as Synthesized markers)?


texts: Matched filter for recieved RF digital packet by LPF

texts:
Im transmitting by cc1350 launchpad microprocessor digital RF Packets, transmitted modulation FSK. bit rate of transmitting is  50Kbouad, frequency of transmitting is 868MHZ . Im sampling the signal(RF PACKETS) that Im transmitting by a dongle SDR-RTL and it samples on frequency 2048000Hz , I want to decode the RF packets that Im transmitting, so I build demodulator and I see the data of my RF packet which it's: (I plot it in matlab plot)

So if you see here on the plot it shows data of my sampled RF packet, and in order to decode it to 0/1 (binary), I want to do LPF in order to smooth/smear much more my received Rf packet "to look like smoothly sinusoidal"
any help how can I do that?
I mean by demodulator is that I've done two matched filters (Mf1, MF2) as I shown in the photo below:


texts: Use of the Dirac delta as a sampling operator

texts:
As we are all taught in signal processing classes, the Dirac delta is not a normal function but a generalized function. It is sometimes defined as either the limit of some normal function whose support goes to zero as some parameter goes to zero, while keeping unit area. This apparently leads to the usual representation of the Dirac delta function as being defined and only having meaning under the so-called sifting integral,
$$
\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-t\right)d\tau=x\left(t\right)
$$

Yet, pick up any book or paper on signal sampling and you will see the sampling process defined as a modulation process with the delta standing alone as a normal function, multiplying the analog signal to be sampled, like this:
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)
$$
Why is this legal? Why isn't the sampling process modeled as
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau
$$
or
$$
x\left(nT\right)=\intop_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau\text{?}
$$

Similarly, in the frequency domain, we know that a linear system's response to a sinusoidal input is a common measurement. Using the common definitions for $X$, $Y$, and $H$, we know that
$$
Y\left(f\right)=X\left(f\right)H\left(f\right).
$$
If
$$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
is a sinusoid, then
$$
X\left(f\right)=\delta\left(f-f_{0}\right)
$$
and
$$
Y\left(f\right)=H\left(f_{0}\right)
$$
where we have once again used the stand-alone delta as a sampling operator.

Can someone please clearly explain my confusion on this matter? Why do we seem to be so cavalier about how we use the Dirac delta when modeling sampling?

(While we're in the neighborhood, maybe someone could explain why we're even allowed to think about the Fourier transform of a sinusoid since it it is neither absolutely integrable nor square integrable.)

By the way, I find the discussion at Why is dirac delta used in continuous signal sampling? to be unsatisfying. Can someone do better?


texts:
I'm a complete beginner to object recognition.

I need to identify and count how many of a certain type of Lego, in a picture full of different pieces of Lego, using MATLAB.

The pieces of Lego I need to identify, are of a certain shape, colour and size.

I need suggestions for algorithms that can perform this task, and links to where I can find more information about them.

From the research I've done, so many are focused around facial recognition, and provide no help for me!

I followed a tutorial that used SURF points to identify points of interest, in the scene pictures and the picture of the Lego picture to find. But it wasn't able to pick up on all the aspects of the Lego piece, and therefore didn't have a chance of finding it in the scene.

If there aren't any algorithms that can perform this, how can I at least separate the pieces based on their colours?

Any suggestions would be really handy! Thanks in advance.

Edit:

Sample picture where I will need to identify certain pieces:


texts: Quantify the difference between two audio spectrograms

texts:
I have two speech signals of the same length. They contain the exact same speech content but are recorded from two different microphones. I want to compare the differences between them.

1) Assuming that the two signals are perfectly aligned, I can compute the spectrograms for both of them and calculate a distance/coherence metric (L2 distance or magnitude-squared coherence) between the two spectrograms. Are there any other ways of quantifying spectrogram differences?

2) If the signals are not perfectly aligned - even a misalignment of a few milliseconds - can mess up the L2 distance between spectrograms. What could be some strategies in this scenario?

Thanks for your inputs.


texts: Squaring signal for amplitude demodulation

texts:
I read the following here. Can someone explain/elaborate on this?

  Squaring the signal effectively demodulates the input by using itself
  as its own carrier wave. This means that half the energy of the signal
  is pushed up to higher frequencies and half is shifted down toward DC

I didn't find any good source giving more information or examples of this.


texts: How to Calculate the Color Temperature / Tint of the Colors in an Image?

texts: How can I calculate the (black body equivalent) color temperature of the light source illuminating a given image? Below are screenshots of Adobe Lightroom manipulating color temperature and a shift in the RGB Histograms. Given RGB components of the image, how do I go about calculating it? I should expect a single value - The black body equivalent temperature of the illumination source, right?


texts:
I am trying to understand the real DFT and the DFT and why the distinction exists.

From what I know so far the DFT uses $e^{i2\pi kn/N}$ for basis vectors and gives the representation $$x[n]=\sum_{k=0}^{N-1}X[k]e^{i2\pi kn/N}$$ The sum is written from $k=0$ to $N-1$ for historical reasons I think instead of writing it in a way analogous to the Fourier series with the sum going from $k=-N/2$ to $N/2-1$: $$x[n]=\sum_{k=-N/2}^{N/2-1}X[k]e^{i2\pi kn/N}$$
This relying on a peculiar anomoly of the DFT where high frequencies are the same as negative frequencies: $e^{i2\pi kn/N}=e^{i2\pi (k-N)n/N}$.

Continuing the analogy with Fourier Series the real DFT gives the representation
$$x[n]=\sum_{k=0}^{N/2}\left(X_R[k]\cos\left(\frac{2\pi kn}{N}\right)-X_I[k]\sin\left(\frac{2\pi kn}{N}\right)\right)$$
This can be viewed as pairing $e^{i2\pi kn/N}$ with $e^{-i2\pi kn/N}$ in the DFT representation where the sum ranges from $k=-N/2$ to $N/2-1$. This is very much like the pairing $c_n e^{in\theta}+c_{-n}e^{-in\theta}=a_n \cos n\theta + b_n \sin n\theta$ which connects the two representations of a Fourier Series:$$\sum_{-\infty}^\infty c_n e^{in\theta}= \frac{a_0}{2} + \sum_1^{\infty}(a_n \cos n\theta + b_n \sin n\theta)$$

My question then is why is the DFT so much more prevalent than the real DFT? One would expect that since the real DFT is using real valued sines and cosines as the basis and is thus representing the geometric picture better that people would like it more. I can see why the DFT and the continuous Fourier Transform would be preferred in a theoretical sense as the algebra of exponentials is simpler. But ignoring the simpler algebra, from a practical computational applied viewpoint why would the DFT be more useful? Why would representing your signal with complex exponentials be more useful in various physics, speech, image, etc. applications than decomposing your signal into sines and cosines. Also if there is anything subtle I'm missing in my above exposition I would like to know: I'm puzzled that the DFT is seen as more connected to the continuous Fourier Transform than to the Fourier Series.


texts: Choosing minimum sampling frequency of a filtered signal

texts: I have a time series of acceleration sampled at 1000Hz to capture random impacts on a structure. The ringing of the structure (not of interest) produces noise which a 50Hz filter seems to eliminate while preserving the impulses' amplitude and width. I'd like to downsample the filtered data to reduce data volume. What is the minimum sample rate I could choose? Based on this article, I should choose a sample frequency higher than 2*50Hz (Nyquist) since 50Hz is just the filter cutoff and there is still energy at higher frequencies that I don't want to alias. Is 10*50Hz = 500Hz reasonable?


texts: How can you perform echo cancellation on features of audio signals?

texts: I am working on a project that takes some features of audio signals and classify them with a neural network. However, the audio signals are added with their own echo with the environment. What would be some of the processing techniques that can be applied on the features, including chroma stft and Mfcc.