texts: Looking for a analytical formula to compute the central frequency of a signal analyzed by a discrete wavelet at a given scale

texts:
I am reading this paper by Han et al. (2014). In this article, the authors extract detailed information from geomagnetic data sampled at 1Hz using a daubechies 5 (db5) wavelet: they reconstruct the signal using the coefficients from the 6th level only. The authors explain that the central frequency of the resulting reconstructed signal is 0.01Hz.

I have reprocessed similar data following the same steps and have convinced myself that the central frequency of the signal reconstructed using only coefficients from the 6th level is indeed around 0.01Hz. Below are the PSD plots for both the original and reconstructed signals, generated using Python:

But is there a way to analytically determine the central frequency of the reconstructed signal?

On one hand, I have found explanations stating that, considering a SR of 1Hz, the central frequency of the signals reconstructed using the 6th level should lie within [0.5, 1]/2**6, i.e. [0.0078125, 0.015625]Hz, which is the case. Is this the correct answer? If this was the case, it would mean that the central frequency of the analyzed signal would be "wavelet-independent"...
On the other hand, I know about the formula from Abry, 1997, stating that the pseudo-frequency of a wavelet at a given scale Fa is the central frequency of the wavelet Fc divided by the scale a: Fa=Fc/a. But my undestanding is that this formula yields the pseudo frequency of the wavelet, not the central frequency of the analyzed signal!

Any clarification would be appreciated. Thanks!


texts: We know that (7,4) hamming code has 4 message bits and 7 code bits and can detect any single bit error.  The (12,4) repetition code has 4 message bits and 12 code bits and can detect any single bit error.  I tend to think that since in hamming code we are using less number of parity bits, there should be some disadvantage vis a vis repetition code.  Is it true and if yes, what is the downside of hamming code? (One known disadvantage is that of encoding and decoding complexity.)


texts: To find the impulse repsonse using the difference equation

texts:
A causal linear time-invariant filter has transfer function

a) Denote the input signal by x[n] and the output signal by y[n]. Find the difference
equations for the filter.

f) Find the impulse response of the system. (Eliminate j (the imaginary unit) from
your answer.)

I used the partial fractions method and obtained the following solution.

 However it was time consuming. Is there a way that I can simplify this difference equation to get the impulse response of the IIR filter to get the same solution. Any help will be appreciated(as I am not aware of solving by difference equations)


texts: simulating multipath via tap delays in FPGA

texts:
I have a Xilinx FPGA running at 200MHz with an ADC and DAC on board.  I am taking a 2462MHz sine wave (pretending it is 802.11b chan 11), mixing it down to 40MHz putting it through my FPGA, mixing it back up and looking at it on a Spectrum Analyzer.

I have 4 taps (255 clock delays deep) to play with in my FPGA and am currently only using 2 of them.  I worked out the phase delay on each tap index by saying that the phase for a particular index (i) is phase=i*(2*pi)(2.462e9/200e6).  So index 1 has a phase of 0.13(2*pi) and index 250 has a phase of pi.

Armed with that, I put my signal through with no delay on the first tap, and then delayed the signal on my second tap by 250 clocks (so it should be 180 degrees out of phase).  Yet the signal doesn't disappear on the specAn.  I haven't found a tap that does make the signal cancel out, so that was where I was curious what I was doing wrong.

I was hoping there was something simple I am missing, as opposed to finding out it can't be done without implementing some sort of fractional delay filters (or something to that effect).


texts: Variance of zero-mean signal after convolution for SIR computation

texts:
my goal is to scale desired, interfering signal at the receiver in order to achieve desired SIR (signal to interference ratio) for beamforming (source separation) application.

Let be:

$s(t)$ a known speech signal with zero mean $\mu_s = 0$ and known variance $\sigma^2_s$.
$q(t)$ a known speech signal with variance $\sigma^2_q$.
$h_s(t)$ and $h_q(t)$ two known deterministic acoustic impulse responses.
$y(t) = (h_s \ast s)(t) + (h_q \ast q)(t)$ the reverberant speech signal.

Suppose we have access only to the premix and the mix, that means to
$x_s(t) = (h_s \ast s)(t)$, $x_q(t) = (h_q \ast q)(t)$ and $y(t)$.

To achieve desired SIR at the receiver, I can simply make $x_s$ of unit variance, and scale $x_q$ and then mix again the two quantity accordingly.
However now, how do I scale $s(t)$ and $q(t)$ accordingly.

I think a more general approach is then, given the deterministic LTI filter $h(t)$, what is the variance of $y(t) = (h \ast s)(t)$?

$$ \sigma^2_y = \mathbb{E}[ (y - \mathbb{E}[y])^2 ] =
\mathbb{E}[(h \ast s - \mathbb{E}[h \ast s])^2]$$

Since $s$ is zero-mean periodic speech signal and $h$ is a deterministic filter, then $\mathbb{E}[h \ast s] = 0$.

It follows that $$ \sigma^2_y = \mathbb{E}[ (h \ast s)^2 ]$$
Using the Parceval Theorem and the convolution theorem, I can write
$$ \sigma^2_y = \mathbb{E}[ (H X)^2 ]$$

However if I write everything directly in the frequency domain, as
$$\mathcal{P}(Y) = \mathbb{E}[ H^2 X^2 ] = | H |^2 \mathbb{E}[X^2] = | H |^2 \mathcal{P}(X)$$
where $H$ and $X$ are the DFT of $h$ and $s$ respectively, while $\mathcal{P}(\cdot)$ is the PSD of operator.

And here I am not sure how to continue, since the PSD is defined over frequencies, while the variance of the signal is a scalar.

Thanks


texts: Which equation should I use to compute the Extended Kalman Filter?

texts:
Compute the Extended kalman filter can be done in several ways.

The first one compute the convariance matrix $P(t)$ from the Riccati Equation:

$$\dot{P} = F(t)P(t) + P(t)F^T(t) - P(t)H^T(t) R^{-1}H(t)P + Q$$

The second one compute the covariance matrix $P$ from this equation:
$$ \dot{P} = F(t)P(t) + P(t)F^T(t) - K(t)H(t)P(t) + Q(t)$$
https://en.wikipedia.org/wiki/Extended_Kalman_filter#Continuous-time_extended_Kalman_filter

So which one is it? It dosen't matter right? I still can find $P$ by using:

$$P(t) = P(t) + \dot{P} dt $$

Where:

$dt$ is the sampling time

Or I can use $fsolve$ in MATLAB to solve $P(t)$ from this equation:
$$0= F(t)P(t) + P(t)F^T(t) - K(t)H(t)P(t) + Q(t)$$

What is best?


texts: Demodulating FM from IQ Samples using Simulink

texts:
I am using simulink and have got IQ samples from an RTL-2832 / E4000 usb dongle. I would like to learn how to decode an FM station to audio and attempt this project on simulink. How do I proceed with the IQ Samples towards decoding?

Treat this a  question from a DSP novice.


texts: This is a general question, I am just trying to understand the concept. Why should we use higher-order filters,  other than a 1st order (of any type for that matter, but we can keep the discussion to digital Butterworth filters). I understand that the phase shift reduces dominantly as the order increases, but is it the only advantage of increasing the order of the filter or is there any other?


texts:
On Wikipedia, is written that

  MUltiple SIgnal Classification (MUSIC) is an algorithm used for frequency estimation.

One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

I premise that I'm not an expert on the signal theory, therefore I do apologize if this question it will not be much precise.

Since MUSIC is an algorithm used for frequency estimation, and since purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities,
my question is the following. Is MUSIC a method to detect the number of peaks in a power spectrum? Is it used for this purpose?

References and answers are welcome. Thanks in advance!


texts: Number of antenna to perform MUSIC algorithm

texts: When the MUSIC algorithm is used for direction of arrival (DOA) estimation it can identify up M-1 signals, where M is the number of array elements. MUSIC can also be used for spectral frequency estimation. Does MUSIC require more than one receiver (or more than one receive antenna) for frequency estimation of a single signal?


texts: How to interpret the sign of a cepstral coefficient?

texts: I feel like I've got a decent intuition for what the magnitude of a cepstral coefficient represents (i.e. how much of the signal's energy is contained in harmonics of a given wavelength), but I don't really understand what the sign of the cepstrum represents. I feel like it's important to know since some speech processing applications take the max of the real part of the cepstrum (see here), and I'm wondering if there's a reason why the max is in any way more informative than the min, or if I should instead be taking the max of the absolute value of the cepstrum?


texts: I have been looking around on Google, and I have not been able to find an explanation of how to do baseline placement in chromatography. what are some algorithms for placing a baseline on a chromatogram, how do those algorithms work. (Referral to Python implementations are welcome).


texts: What does the term "railed" mean in signal processing?

texts:
I'm having trouble finding a definition of "railed" that relates to signal processing.

Am I correct in my guess that this term is in fact from this field?

My signal data comes from an EEG device. The lightly documented open source software I'm using doesn't define it, but it shows that term when there is no signal data being displayed.

Is that all it means (no data)? Or does it mean something like the signal being read is too great to be displayed or correctly measured?


texts: Please explain the following Linear Prediction estimate graph

texts:
I am reading this doc https://librosa.org/doc/latest/generated/librosa.lpc.html#librosa.lpc?
I can't understand the graph for linear prediction coefficients. Please help explain it as I a new to this field.
Thanks


texts:
Normally a system transfer function is represented by H(z)

I want to know wether all those system transfer functions are only representing a digital filter?

Or their any other thing/entity in z domain that can be represented by H(z)


texts: Does speaker recognition system need the same sample rate audio?

texts:
I'm a new man in speaker recognition project. I've train a model with some .wav files, its sample rate is 16000.

I want to test this model with my record, but my record's sample rate is 8000

Doed it need the same sample rate for training and testing speaker recognition system?

Hope to have your answers,

Thanks


texts: How do I play and record audio at the same time?

texts:
So this question is more related to Matlab itself.

Some background info: I'm trying out the Exponential Sine Sweep method to obtain the Room Impulse Response with my laptop speakers and microphone for testing.

I can generate the Sine sweep perfectly but I have trouble playing and recording it without a lot of lag. I'm now trying to use "audioPlayerRecorder" from the Audio System Toolbox but when I do:

f1 = 80;       % Start frequency (Hz)
f2 = 22000;    % Stop frequency (Hz)
T = 3;         % Sweep duration in seconds
fs = 44100;    % Sampling frequency

[sweep, inverse] = ESS_generator(f1, f2, T, fs);

playRec = audioPlayerRecorder(fs);
recording = playRec(sweep);

I get an error on the last line saying:
Error using coder.internal.assert (line 33)
No full-duplex audio device detected

But I know that my laptop is able to use the microphone and speakers simultaneously so I'm thinking maybe the error gets thrown because of a missing configuration?

I'm running the code on Linux with the ALSA audio drivers. I'm using external speakers through the AUX audio port, and the default laptop microphone.


texts: How to get discrete values from a nosiy data set that have a trend?

texts:
I have a sensor that produces values over time (appr. 50/sec). I visualized these values (X axis is the time), this can be seen on the picture, measurement A and B - these are two independent measurements.

I'm a total beginner, I don't even know how to phrase the question. I just need to get... I don't know... "discrete" (?) values from this noisy data set.

I need something like this (this is happening say under 10 seconds):

Up... Up... Up... Up - I'm only interested in the "ups".

I was thinking of summarizing the values and when it's positive that's an up and when it's negative that's a down. But I don't know where to start and until when do I have to summarize the values...

or if is there (I'm sure there is) any other methodology that is used for this kind of problems.

UPDATE:

After talking with Dan Boschen and applying what he recommended (moving average), I got lot smoother result:

Still got difficulties to ignore the false movements, but this is a one-step-forward. I'll continue this topic in a different question.


texts:
I would like to determine the SINAD of a distorted signal from a measurement and have troubles to get consistent data. The problem is (at least to my understanding) that I only have a very limited number of signal periods and I'm not sure, whether I have enough data to provide this measurement at all.
In order to determine the SINAD I need the power of my signal as well as the power of "all the rest", i.e. the power of noise and distortion. To obtain these powers, I took my oscilloscope data and calculated the spectrum via an FFT. I didn't use any window function (i.e. I used a rectangular window) to calculate the FFT since I know the fundamental frequency I expect from my signal and can choose an integer multiple of the fundamental period as my window width.
From the frequency bin at my desired frequency I determined the fundamental amplitude. The noise + distortion power was then calculated by subtracting the signal's RMS value from the total RMS value according to

From these two RMS values I then calculated the SINAD according to:

The problem is now that my results vary depending on the number of fundamental periods (I took between 1 and 33 periods, I don't have more) from the measurement I take to calculate the FFT, where increasing the number of fundamental periods degrades the SINAD.
This makes sense to me to a certain extent since if I extend the window size of my input data, the frequency resolution of my spectrum increases. This means that parts of the noise + distortion in the vicinity of the fundamental bin which were combined into this bin for a course frequency resolution are now contained in the neighboring bins of the fundamental and treated as noise/distortion and not as part of the signal anymore.
I couldn't really find a lot on how to eliminate this issue. Googling on how to measure SINAD mostly leads to the definitions themselves. Is there a way to determine the SINAD for the limited amount of signal periods? How would I have to calculate it? I also tried the matlab function sinad(), however, it has similar issues (although with this function, the dependence on the number of signal periods is less compared to my own calculation).


texts: Why the inverse discrete fourier transform of the Ricker pulse isn't the same as the Ricker pulse in time domain?

texts:
Question
I'm trying to use Python's scipy library to compute the IDFT of the Ricker wavelet function and compare it with the analytical time-domain version of the same function. When I compare the results, they do not match. I want to understand why.

Details
I'm using Python for the experiment. The Ricker pulse function in time domain is defined as follows:
def ricker_pulse_time(w0, t):
    return (1. - .5 * w0**2 * t**2) * np.e ** (-.25 * w0**2 * t**2)

def ricker_pulse_freq(w0, w):
    return (2. * w**2) / (np.sqrt(np.pi) * w0**3) * \
        np.e ** (-w**2 / w0**2)

So, for 2Hz, for example:
hz_to_rad = lambda x: 2. * np.pi * x

w0 = hz_to_rad(2.0) # [Rad/s]
time_step = 1e-6 # [s]
t = np.arange(-2.0, 2.0, time_step)
y = ricker_pulse_time(w0, t)
N = len(t)

omega = np.arange(0.0, 100.0, 1e-4)
yy = ricker_pulse_freq(w0, omega)

plt.plot(t, y, 'r-');
plt.show();
plt.plot(omega, yy, 'b-');

But, when I try to make the ifft of the ricker pulse in the frequency domain, I don't seem to get the expected curve in time domain, as above...
omega = np.arange(-300.0, 300.0, 1e-4)
yy = ricker_pulse_freq(w0, omega)
plt.plot(omega, yy);
plt.show();

y_ta = fft.fftshift(ifft(yy))
tt = np.arange(-2.0, 2.0, 4. / len(y_ta))
plt.plot(tt, y_ta.real, 'b-');
plt.xlim(-.005, +.005)
plt.show();

I don't understand the reasoning behind this, and what should I do to get the same results. I understand that, perhaps, I won't get the same amplitude of signal, but I want to at least get the same shape...


texts: I want to obtain the energy of this signal $$x_i(t)= e^{j2\pi(2i-1)f_0t}$$ where $f_0 = 1$ Hz, $i= 1,...,4$ and $j$ is a complex. I think I need to use the Fourier transform or the Parseval theorem but I'm not sure how to apply it.


texts: Regularization for Inverse Problems using the Singular Value Decomposition (SVD)

texts:
I am reading these lecture notes on Optimisation and Inverse Problems in Imaging, and I have difficulties understanding how figures on page 20 (Figure 3.2) or page 21 (Figure 3.3).
Precisely, I don't understand what numbers on horizontal and vertical axes mean. I would appreciate if you could explain me this. Here is the code for Figure (3.3).


texts:
I have started analyzing a signal from a MEA experiment that consists of reading from electrodes inserted in a brain slice. Anyway, the signal is sampled at 2 KHz but the readings we have are very strange. It looks like it is very noisy and we cannot understand why it is so much noisy so one of our doubts was if the sampling frequency was to high that it incorporated many artifacts or folded down frequencies into our first half of the Fourier transform frequency spectrum.
I will upload a figure so that you can see it and opinions are appreciated!

The time domain plot of the signal is the shown below, as you can see, it is very noisy but we have the important part, the peaks.


texts: In this article about pre-emphasis on here, its use is defined as "to compensate for the average spectral shape". My question is, how is the average spectral shape calculated? This is the graph included in the article:


texts: Why multi-path channel has linear phase within the coherence bandwidth?

texts:
I read some sentence that said "It is usually assumed that the multi-path channel has a constant gain and linear phase within the coherence bandwidth".
I can't understand what it mean " linear phase ".
To my best of knowledge, difference path delay would caused difference phase shift in the receiving signal, but it's would be rarely that different path delay is linear, or it should be said that there is almost no such situation. Therefore, I am very curious about what it mean about "linear phase".
Can someone please explain it?
Thanks in advanced.


texts: DFT for something similar to convolution

texts:
I have the following problem:
Let $x$, $y$ be finite real valued sequences defined on $0, \ldots, N-1$ and let $g$ be a non negative integer.

Define:

$$z_n=\sum_{k=0}^n x_{n-k}y_k, \quad \text{also on $0, \ldots, N-1$.}$$

Note that the sum is from $0$ to $n$, not $N-1$. In addition, the DFT of $y$ is known in closed form.

Is there a way to write $z$ as some cyclic convolution, so that with the help of the convolution theorem $z$ can be calculated in $N \log N$ instead of $N^2$?

I tried following the convolution theorem proof but I get stuck:

\begin{align}
z_n&=\sum_{k=0}^nx_{n-k}y_k= \sum_{k=0}^{N-1}x_{n-k}\cdot \theta (n-k)\cdot y_k\\
\Rightarrow \tilde{z_q}&=\sum_{n=0}^{N-1}\sum_{k=0}^{N-1} x_{n-k}\cdot \theta (n-k)\cdot y_k e^{-2\pi iqn/N}\\
&=\sum_{n=0}^{N-1}\sum_{k=0}^{N-1} x_{n-k}\cdot \theta (n-k)e^{-2\pi iq(n-k)/N}\cdot y_k e^{-2\pi iqk/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n=0}^{N-1} x_{n-k}\cdot \theta (n-k)e^{-2\pi iq(n-k)/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n'=-k}^{N-1-k} x_{n'}\cdot \theta (n')e^{-2\pi iqn'/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n'=0}^{N-1-k} x_{n'}e^{-2\pi iqn'/N}=?????\\
\end{align}

The problem is that the second sum depends on $k$ so the double sum doesn't factor to the product of DFTs.

What am I missing?


texts: OpenCV: Is NormalBayesClassifier::train function running ok?

texts:
I am have the code from here. I have changed the images with some of my own and I have seen that SVM is not working. That is not a problem because I have found on another question the answer. Now I have changed the evaluation to NormalBayesClassifier. I am waiting for more than half an hour for the classifier to train. Is the NormalBayesClassifier::train function so slow? I have put 7 images in one class and 7 in the other (the size is about 640x480). (I have just 2 classes). I have changed the size of the images to much more smaller: around 180x100 and I have cropped the area where was too much background, but I am still waiting for training. What could be wrong?

Another thing: Is trainer working? I can do a lot of things meanwhile the training and my comp is not running slow...


texts: MATLAB Image-Processing Toolbox: asking for a roadmap!

texts:
I plant trees in the north of BC for a living and have next to
    none sp experience, please be gentle :^]

TLDR:
The problem is essentially this:
 tracking and calculating the "height" of three differently-colored regions whose proportions change relative to each other over a number of images using k-means segmentation and color-histogram vector.

(see fig. below)

It's a little chemistry experiment which kicks off with a transparent beaker full of light-blue-aquamarine foam, then a transparent oil is added and over the next few minutes the light blue foams condenses into a distinctively-darker blue homogenous liquid at the bottom of the beaker and there is no light-blue foam left whatsoever. The boundary between the foam and the emulsifier(oil) is.. well.. foamy at all times as i tried to indicate in the drawing, but i think a blur with an appropriate kernel should solves that pretty well.
 We have the full control over the setup(any lighting, any angle and number of photos, etc.)

They want me to grab the MATLAB image-processing toolbox which I, unfortunately, happened to have and given a set of oh let's say 360 snapshots of the above experiment plot the "height" of the boundary lines between the oil, the light-blue foam and the dark blue liquid at the bottom against time. Without sticking a ruler next to the vessel! Bummer!

How would a signal-processing guru like yourself have gone about this? Terms like "k-means segmentation" and "dominant vectors of a histogram" prevailed in the conversation with the engineer who has imparted this task to me, but knows less MATLAB than I do. He has also made it abundantly clear that he wants it done in MATLAB. Woulda been a simple problem in opencv, but oh well!

Could you possibly throw like a bullet-list of functions/ sequence of steps one would undertake in solving something like this? So i know what to research at least.
 It has to be k-means, histogram and matlab unfortunately.
The difference between going this road and any other for me is like between getting laid and getting shot, folks. I count on your understanding and am indebted in advance!

artie


texts: Length of the FIR filter in LMS algorithm

texts: I'm working on an interference cancellation problem where interfering signals are two linear FM chirps of different slopes and the data signal is a spread spectrum signal. I wanted to know how does the length of the filter affect the convergence ?


texts: How does one describe the “frequency response” of a mixer?

texts:
Suppose one uses a vector network analyzer to measure the $S_{21}$ parameter of an RF mixer over a frequency range. Would it be appropriate to call this the “frequency response” of the mixer? I find this confusing because mixers are not LTI systems, and it is therefore meaningless to speak of the frequency response of a mixer.

But, from what I understand, experimentally measuring the $S_{21}$ parameter through a LTI system such as a lowpass filter, amplifier, or cable precisely yields the frequency response of the system. In that case, what does the $S_{21}$ parameter through a mixer (not LTI) correspond to?


texts:
In communication systems, the transmitted signal is often passed through a root-raised cosine filter to avoid, in textbook language, intersymbol interference (ISI). An identical filter at the receiver end processes the signal so that the combined response of the two filters is a raised cosine filter.

The raised cosine filter has two properties: roll-off factor and oversampling rate.

I have a bipolar sequence ($\pm 1$ values) that I am supposed to process through a raised cosine filter. It is not clear to me what exactly is intersymbol interference here, and how can I visualize it. Here the symbol has only one bit. How should I choose the roll-off factor and the oversampling rate?

-ryan


texts: Find branches in skeleton image using matlab

texts:
I have a skeleton image. I am able to detect branchpoints, endpoints and able to label image using BWMorph, BWLabel. Now I want to differntiate between skeleton center line and branches. In short I want to detect branches in skeleton. I want to perform operations on skeleton branches and I do not want to consider center line in these operations. Can someone help? Following is skeleton image.
http://www.flickr.com/photos/92388309@N03/8397759970/in/photostream/
Arrows in following image show center line of skeleton.
http://www.flickr.com/photos/92388309@N03/8404495186/in/photostream/

Following is code I have used.

mn=bwmorph(y,'branchpoints');
[row column] = find(mn);
branchPts    = [row column];
endImg    = bwmorph(y, 'endpoints');
[row column] = find(endImg);
endPts       = [row column];
figure;imshow(y);
hold on ;
plot(branchPts(:,2),branchPts(:,1),'rx');
hold on; plot(endPts(:,2),endPts(:,1),'*');
% Labeling the Branches
branches = imsubtract(y,mn); % set branch points to zero
figure; imshow(branches);title('no bp')
branchesLabeled = bwlabel( branches); % label connected components
vislabels(branchesLabeled);
% Calculation of Length of Branches and Circular Neighbourhood Method for
% for detection of normal and abnormal branches.
sts = regionprops( branchesLabeled,'Area', 'Perimeter','MajorAxisLength','Centroid' );
figure
imshow(branchesLabeled)
hold on
% Loop for circles
for i=1:size(sts)
r= sts(i).MajorAxisLength/2 ; %desired radius
centerx = sts(i).Centroid(1);
centery = sts(i).Centroid(2);
th = 0:pi/50:2*pi;
xunit = r * cos(th) + centerx;
yunit = r * sin(th) + centery;
plot(xunit, yunit);
end


texts:
I have a signal of the form $s(t)=A(t) \sum \cos(\omega_i(t)t +\phi_0) + n(t)$, where $n$ is gaussian noise.

I can only read the signal+noise and thus can not separate them.

I want to phase shift the signal to $A(t) \sum\cos(\omega_i(t)t)+n'(t)$ and I am at a loss on how to do this. During lectures / courses I've always done phase shifts simply by multiplying with $e^{i\phi_0}$.

My signal is of the form $s(t)=A(t) \sum \cos(\omega_i(t)t +\phi_0) + n(t)$, as opposed to $s(t)=A(t) \sum e^{\omega_i(t)t +\phi_0} + n(t)$, so I can not simply multiply with $e^{i\phi_0}$

Is there any way to do this? I'm asking because I am interested in cross-correlating two signals.


texts: My X ray images have black background, which turn gray after low-pass filter. How to recover blackness?

texts:
I am applying a low-pass filter to X ray images, which changes the pixel distribution such that the background goes from black to gray. In the attached figure, this amounts to the mode increasing from 0 to 50. I want to postprocess the filtered images to have a black background, while minimizing additional information loss. One idea I have is to compute the mode of the filtered image, then reduce the brightness of the image by this amount. Is there a less ad-hoc way of achieving my goal?

Here is my procedure for low-pass filtering:

Apply DFT to the input image
Shift DC component to the center
Apply a radial filter, setting all frequencies outside of radius $r$ from the center to zero
Inverse shift
Inverse DFT
Normalize pixel values to be in $[0, 1]$


texts: Modelling the amplitude and phase through an amplifier

texts:
I am writing a software simulator to simulate a real physical acoustic interferometer. I now want to write the module to model the microphone amplifier(s). I'd like to be able to specify the frequency response of the amplifier both in terms of amplitude (gain) and phase. So I will feed in a sampled waveform and get a sampled waveform out that is modified according to the properties of the amplifier circuit.

I'm no expert in signal processing but happy to read the background maths. But to save me time going down dead end ideas could someone please tell me the standard method of doing this and what maths I have to read up on to be able to understand it?


texts: How do i select my frequency range to plot when fft is made on time domain

texts:
I have a vibration signal that i need to convert from time domain to frequency domain using fft in python.

from scipy.fftpack import fft
yf =  fft(df["x"])
plt.plot(df["x"])

And i would like to plot it without DC value at 0Hz. So i neglected yf[0] and took N/2 frequencies to plot as per Nyquist theorem.

But I am having little trouble how to calculate x axis frequencies. how do i select my range of frequencies to plot? Also the amplitude is plotted from 0 to 25. Does this amplitude (Frequency domain) corresponds to sum of amplitude of a particular signal in time domain? How do i know which amplitude from time domain corresponds to amplitude in frequency domain?


texts: wide sense stationary of dynamic process

texts:
I am trying to understand the definition of wide sense stationary on my own and probably have some silly questions. Wikipedia says, wide sense stationary is a process with constant mean and  autocorrelation function over time. Now, lets assume a process,
$$X(n) = 0.9X(n-1) + V(n) + U(n),$$
where $V$ is white noise and $U$ is an external input.

Lets assume $U(n) = 0\quad \forall n$, i.e. $X_1(n) = 0.9X_1(n-1)+V(n)$. In this case can $X_1$ be called WSS? What does mean being constant actually implies? Can I assume $0.9X(n-1)$ is mean since it is deterministic and say mean is decreasing over time, and this is not WSS?
Lets assume $U(n) = C \text{ const.}$, i.e. $X_2(n) = 0.9X_2(n-1) + V(n) + C$
If $X_1$ was WSS, I believe this should be WSS as well. Please confirm.
If $U(n)$ is an external signal changing over time, i.e. $X_3=X$ then would still be WSS, if $X_p= 0.9X_p(n-1)+V(n)$ were originally? In this case $X$ is depended on a signal which is not part of the original process.


texts: PLL for Phase Demodulation and Carrier Tracking

texts:
I've been playing around with using a PLL for phase demodulation and it's working pretty well. The PLL locks to the carrier and the error term is the demodulated message signal. However, depending on what the underlying message signal is, it is sometimes necessary to configure the PLL to have a very small loop bandwidth - the effect is to prevent the PLL's internal state to respond too quickly to deviations caused by the message signal.

On the otherhand, it's very useful for the PLL to automatically handle locking to the carrier at an offset frequency and/or track doppler - in this situation, a larger loop bandwidth is essential.

With this in mind - is it a valid approach to have a PLL with a default acquisition mode with a large loop bandwidth, and then once lock is detected, transition to a smaller loop bandwidth to give good phase demodulation performance? One idea I had for this is computing the frequency estimates via the loop's phase (e.g. $f_{loop}[n] = \phi_{loop}[n] - \phi_{loop}[n-1]$) vs the detected input sample phase ($f_{i}[n] = \phi_{in}[n] - \phi_{in}[n-1] $), if these are about the same (after averaging/lpf) the loop would be considered locked and the bandwidth reduced.

Are there any well known/used approaches to take in this problem - (solutions beyond the single PLL are welcomed too).

thanks.


texts:
I have a GMSK modulator in MATLAB and wish to generate AWGN with a specific SNR, and I'm having some issues figuring out just how.

I have found this step-by-step process to generate zero-mean Gaussian noise with a specified power spectral density W0:

Form Gaussian distributed random variable: w = randn(1,N)
Map top zero mean: w = w - sum(w)/N
Compute average power: Pw = sum(w.^2)/N
Form w = w.*sqrt(W0*fs/Pw)

Note: As pointed out by Deve, steps 2 and 3 should not be necessary, since randn() generates a random variable with zero-mean and power 1. However, it is possible that the mean and power may be slightly off, since the vector is of finite length. For my application, neglecting these steps give a slightly less accurate SNR.

So, how do I determine which power spectral density I need to achieve a specified SNR? I understand that SNR is defined as $SNR=P_s/P_n$. Measuring the signal power using $P_s=\frac{1}{N}\sum_{n=0}^{N-1} s[n]^2$ and then isolating $P_n=P_s/SNR$ seems trivial, but how do I relate $P_n$ to $W_0$?

Thanks in advance!

Edit: Using Deves answer, I have written the following code:

s           = signal;
L           = length(s);

% Convert SNR from dB
SNRdB       = 10;
SNR         = 10^(SNRdB/10);

% Measure average power of signal
Ps          = sum(s.^2)/L;

% Calculate wanted noise power
Pn          = Ps/SNR;

% Generate random vector and ensure 0-mean
w           = randn(1,L);
w           = w-sum(w)/L;

% Scale to wanted power
Pw          = 1/L*sum(w.^2);
w           = sqrt(Pn/Pw).*w;

% Measure resulting SNR
Pw_meas     = 1/L*sum(w.^2);
SNRdB_meas  = 10*log10( Ps/Pw_meas ); % Gives 10.0000

And here's a plot in time domain for anyone interested:

It doesn't seem like a lot of noise - but I guess the calculation of SNRdB_meas is quite foolproof.

[PSD plot removed]


texts: Adding AWGN noise with a correct noise power to the signal

texts:
I would like to add AWGN noise to my BPSK and QPSK signal. So assume the BPSK signal is x1 and QPSK signal is x2. The EbN0 is $10\textrm{ dB}$. Therefore I will first convert the EbN0 from $\textrm{dB}$ to linear scale. Assume Eb=1 and Es=1.

I would like to know can anyone explain to me, should the noise power be N0 or N0/2?

I am really confused since according to this link
N0 is used to add noise to one dimensional noise and N0/2 is used to add noise to two dimensional signal. Therefore if Eb=1 and EbN0=10 $\textrm{dB}$

Eb=1;
N=1e6;

EbN0=10;
EbN0_lin=10^(EBN0/10);
N0=Eb/EBN0_lin;
noise=sqrt(N0)*randn(1,length(x1))

or should it be

noise=sqrt(N0/2)*randn(1,length(x1))

The same thing for QPSK, is it

Es=1;
EsN0=10;
EsN0_lin=10^(EsN0/10);
N0=Es/EsN0_lin;
noise=sqrt(N0)*randn(1,length(x2))

or

noise=sqrt(N0/2)*randn(1,length(x2))

I was thinking that according to convention the N0 is sigma^2/2 however the link above is saying different things.

Can anyone explain the correct way of calibrating the noise?
And explain if the above reference is incorrect or not?


texts: Iterative Hard Thresholding always thresholds same indices

texts:
I am confused by the fact that the thresholded indices in IHT do not change during the recovery. I used the code from this question and also added the condition that $$\vert \vert \Phi \vert \vert_2 < 1$$ by

# Generate problem instance
A = np.random.randn(n,N)
# normalize to columns to have unit norm
A = A/(np.linalg.norm(A,2)+0.1)

Is it correct that always the same indices are kept/thresholded? For example:
$$x_0 = [ 0.1923, 0, 0, -0.2234, 0 ]$$
$$x_1 = [-0.1234, 0, 0,  2.6289, 0 ]$$
$$x_2 = [ 3.9432, 0, 0, -0.0012, 0 ]$$
$$x_3 = [-1.2345, 0, 0,  0.1298, 0 ]$$
$$ ... $$
for any $ A \in \mathbb{R}^{4 \times 5}$. If IHT will always keep the same indices being not zero, I could find the best approximation by solving a simple equation system, therefore I guess there is an error in the code, isn't it?


texts: Grayscale conversion before applying sobel gradient operator

texts:
I am currently replacing a hand-made image processing library by OpenCV. The programmer who worked on this library had implemented the sobel gradient function. This implementation is right except that he first converts his image to grayscale with the formula thereafter (the arithmetic mean):

Y = 0.333 * R + 0.333 * G + 0.333 * B

In OpenCV, I see that they use the formula:

Y = 0.299 * R + 0.587 * G + 0.114 * B

To convert an image to grayscale.

Is the hand-made implementation wrong or can it be motivated by some purpose. Should I always use the OpenCV conversion version?


texts: Using Spreading spectral with Maximal likelihood

texts:
As asked here, HERE  if we have the signal

x =
[0.7 + 0.7i;
0.7 - 0.7i;
-0.7 + 0.7i;
-0.7 -0.7i];

Which was spread over code c and transmitted over channel H whose dimension is [4x4], so the convolution of signal after spreading will become:

r = reshape(H*reshape(x,4,[]),[],1);

comparing it with signal without spreading, it was

 r = H*x;

That was explained well in the above link.

My question, suppose I am using Maximum likelihood estimation, in case if we didn't spread the signal, we will check the likelhood compared with the channel H, but what's about after using spreading ? how will become the channel ? It supposed to be a vector of [16x1], is that right ? but how will it be ?

thank you!


texts: How to apply "Substitution Method" in Difference Equations?

texts:
How to get the value of $y_p(n)$ using substitution method?
What exactly do we need to substitute to get the value?


texts: Expectation maximization of moving average with binary source input

texts:
I am trying to do blind system identification of a univariate linear FIR model:

I am unsure if the approach is correct or not and any help to further proceed with the maximization will be great. Thank you.


texts:
I have a signal $X(t)=\sum_{n=-\infty}^{\infty} Z_n \delta(t-n\tau)$, $Z_n$ is a random variable with equal possibility of $\pm 1$ and I know the power spectrum of this signal is $\frac {1}{\tau}$ from the solution.

I have try to calculate the power spectrum by the auto-correlation and then calulated by its Fourier transform. As following

\begin{align}
R_X (t+\tau, t) &=E\left[X(t+\tau) X(t)\right]\\
&=E\left[\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} Z_n Z_m \delta(t+\tau-n \tau)\delta(t-m\tau)\right]\\
&=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t+\tau-n \tau)\delta(t-m\tau)\\
\end{align}

Then, I have no idea how to calulate its Fourier transfrom to get the solution as simple as $\frac {1}{\tau}$. I don't know which part is gone wrong, appreciate any help.


texts: How to Combine / Cascade two `3 x 3` Filters into One `5 x 5` Filter

texts:
If you have a uniform 3x3 box filter T which is:

1 1 1

1 1 1

1 1 1

And an 3x3 laplacian filter W which is:

1 -2 1

-2 4 -2

1 -2 1

Can these 2 filters be combined into one 5x5 filter, the order being first T and then W ? (The combination should be possible) And what is the math behind this combination. I can make a 5x5 filter of all one's and after that stride the W filter over that 5x5 and add everything up, is it that simple?


texts: OpenCV Motion Analysis: Optical Flow vs Motion History

texts:
OpenCV Motion Analysis Documentation lists multiple techniques for motion tracking:

Sparse Optical Flow | calcOpticalFlowPyrLK,
Dense Optical Flow | calcOpticalFlowFarneback, and
Motion History Image | calcMotionGradient, segmentMotion, calcGlobalOrientation.

I need help in understanding their pros and cons. Which ones should be preferred under what circumstances?

Context:
I am trying to track lane departure/lane changes. Any abstract/generic response would be good enough. Thanks!


texts:
I need to transform this Laplace function to the z-domain:

From the answer I received:

$s=(1-z^{−1})/T$

Then substitution into my Laplace function would give:

$t(z) = 2R/(m*(1-z^{−1})/T + 2R)$

Is that it?

Then $z$ is the input signal's current sample? Is $z^{-1}$ is the input signal delayed by one sample?

ie. Is it:

output_of_function = 2*R / (m* (1 - input_signal_1)/T + 2*R)

What basic mechanism might this function expect to perform? Eg. a LPF? Or what?

Does the fact that you're subtracting $z^{-1}$ from 1 imply that the input must be between -1 and 1?

Thanks. (I edited this in reply to the answer given from Matt L which helped.)


texts: A few questions about the output value of dsp.SpectrumAnalyzer

texts:
I am testing the output values of dsp.SpectrumAnalyzer for a sine wave mixed with white noise. The $\sigma^2$ of the white noise is 1e-4, so I use 0.01*randn(1024,1) to generate it. The sine wave is generated by dsp.SineWave with default amplitude (=1). The sampling rate is 44100 Hz. If the 'SpectrumType' property of dsp.SpectrumAnalyzer is set by default ('Power'), the scope output is as follows:

According to the analysis in "Convert the Power in Watts to dBW and dBm" section of Estimate the Power Spectral Density in MATLAB, The peak is at 26.9897 dBm, which is what the above screenshot shows us. But other values in the figure confused me. I will ask them one by one.

(1) The power spectrum of the white noise is $\sigma^2$=1e-4 (right?), which is 10*log10(1e-4/1e-3)=-10dBm according to the definition in the above link. But the figure shows that the power spectrum of the white noise is at -40dBm. How is this value obtained? PS: if we calculate the power spectrum in dBW, 10*log10(1e-4)=-40 which is close to what the figure shows, but here I am using dBm unit!

(2) If I change the 'SpectrumType' property of dsp.SpectrumAnalyzer to 'Power density', the output of the System object is as follows:

According to the explanation of this property in documentation dsp.SpectrumAnalyzer System object, now the peak of the sine wave should be "the magnitude squared of the spectrum normalized to a bandwidth of 1 hertz.". The spectrum is 0.5, so its square is 0.25, what I expected to get in dBm should be 10*log10(0.25/1e-3)=23.9794dBm, but what I can see in the figure is 12.6948 dBm. How is this value obtained?

(3) Let's look at the power spectrum density of the white noise again. According to the definition mentioned in the above question, the value should be 10*log10(1e-4^2/1e-3)=-50dBm. I'm not sure if the calculation is correct, but it is close to the figure, although there seems to be some distance. The question is, the section "Convert the Power in Watts to dBW and dBm" of Estimate the Power Spectral Density in MATLAB says the $P_{whitenoise}$=1e-4*(22050/21.53)=0.1024, but if we calculate the dBm, we get 10*log10(0.1024^2/1e-3)=10.2060. To my surprise, the webpage then says "the total power of the white noise is -39.87 dBm." How is this -39.87 dBm obtained?

I'm totally confused about these values. I will greatly appreciate it if you can help me clarify them. Thanks a lot!

The following are some reference code for your experiment. I am using the latest Matlab R2016a.

Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024, 'PhaseOffset',10,...
    'SampleRate',Fs,'Frequency',5000);
SA = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectrumType','Power',...
    'PlotAsTwoSidedSpectrum',false,'Window','Hann',...
    'ChannelNames',{'PSD of the input'},'ShowLegend',true);
for Iter = 1:5
    Sinewave1 = step(Sineobject1);
    Input = Sinewave1;
    NoisyInput = Input + 0.01*randn(1024,1);
    step(SA,NoisyInput);
end


texts: Given h(t)~H($\omega$) Fourier transform pair, what is the counterpart of H(-$\omega$)

texts: Thanks in advance. Any help will be much appreciated.


texts: What components that I should have in using PCA to denoise noisy ECG signal?

texts: I know that I have to have more than one component to use PCA method to process some signals. But I only have one noisy ECG signal that has to be denoised. What other components in PCA beside the noisy signal itself that I should have in order to denoise the signal? Is it an artificial noise or a clean ECG signal?


texts: I am working for my class project on making animal repellent , I need to know how I could go about taking a normal audio file and shift its frequency/Pitch very high to make it an ultrasonic sound.


texts: The Instantaneous Frequencies from IMF Aren't Consistent with Its Spectrum


texts: FM synthesis algorithm is not very efficient

texts:
I have very little knowledge of DSP, I'm trying to implement an FM synthesis algorithm in code / software.

So far I've successfully implemented a FM algorithm using the following:

sine (hz * (integrate (s 0 t)))

where hz is the frequency, s is the 'signal' and t is the time. And integrate is the function I've implemented that does integration using the "Trapezoidal rule".

However if I then try to generate a bunch of samples (amplitude over time), it has to run the integrate function n times for each sample - which adds up to a lot of operations I assume and hence is rather slow.

Are there efficient algorithms I could look into? Or alternative implementations worth considering?

One approach I've considered is instead of integrating the modulating frequency on each sample, is being able to calculate these discrete and allow them to be reused between samples.

So instead of integrating:

sample 1: 0s -> 1s
sample 2: 0s -> 2s
sample 3: 0s -> 3s

Do:

sample 1: 0s -> 1s
sample 2: integrated(sample1) + (1s -> 2s)
sample 2: integrated(sample2) + (2s -> 3s)


texts: How is Wi-Fi transmitter sends complex number to receiver?

texts:
I am recently studying physical layer the Wi-Fi (802.11), after doing the FFT, the signal bits are transformed from '0' and '1' to complex numbers, but the complex number is not exist in real life, so my question is:

How is transmitter send those complex numbers?
What is the transmitter actually send?


texts: Develop a simple image classification system by using Laplacian of Gaussian (LoG)

texts:
I am tring to write a simple CAD system to classify some images within two groups by using Laplacian of Gaussian (LoG). I am using scikit-image for this tasks and I want to use DNN in keras to train a classifer.

The problem is that I am not sure how to create feature list by using LOG.

My images could have different sizes so LOG could have different feature set.

I am using below codes

from math import sqrt
from skimage import data
from skimage.feature import blob_dog, blob_log, blob_doh
from skimage.color import rgb2gray

img1T = Image.open( "C:\Users\Public\h1T.jpg")
img2T=Image.open( "C:\Users\Public\h2T.jpg")
data_img1T = np.asarray( img1T , dtype='uint8' )
data_img2T= np.asarray( img2T, dtype='uint8' )
data_gray_img1T = rgb2gray(data_img1T)
data_gray_img2T = rgb2gray(data_img2T)

blobs_log_img1T = blob_log(data_gray_img1T, max_sigma=30, num_sigma=10, threshold=.1)
blobs_log_img2T = blob_log(data_gray_img2T, max_sigma=30, num_sigma=10, threshold=.1)

when getting the LOG feature for tow images:

>>> blobs_log_img1T.shape
(73L, 3L)
>>> blobs_log_img2T.shape
(51L, 3L)

now these two images have different shapes, but I need to find a way to use identical number of features for all images. Do I have to resize images to be of identical sizes or I have to take another approach?

Thanks,


texts: Advantage of applying a Window Function in Analysis and Synthesis of STFT?

texts:
I'm studying Short-time Fourier Transform and I learned that if you apply a windowing other than rectangular, you need to overlap the frames---for example, 50% overlap for Hamming family---to reconstruct the signal properly (like add them up to 1). However, I only need to apply this window function during the Analysis part and I can reconstruct the signal properly during synthesis.

But in some literature (like Audacity's Noise Reduction) I saw some methods that apply a second window function during the synthesis part and employ a 75% overlap-add. What's the advantage of doing the window function twice (analysis/synthesis) with 75% overlap compared to a single window function (analysis) with 50% overlap? I can't seem to find any explanation justifying the advantage between these two windowing techniques.

Edit: I also tried to implement the 75% overlap with Hann window applied on both analysis and synthesis part and I was able to yield a x1.5 of the original signal's amplitude. Did I do it correctly?


texts: Haar Wavelets vs STFT for graceful audio quality degradation (newbie question)

texts:
I'm playing around with FFT (using Naudio library) and with Wavelets (using TarsosDSP library) trying to get better understanding of both of them.

I've noticed that when I zero out FFT/Wavelet coefficients (for the same window size - and then perform inverse transform and listen to the result) starting from low to high energies and compare the loss of audio quality Wavelets seems to suffer less. They have less noise for the same number of zeroed out coefficients (or require less non-zero coefficients to have roughly the same amount of noise) and require less non-zero coefficients for audio still to be recognizable.

Is it inherent property of Haar Wavelets or have I messed up somewhere and FFT should perform equally well in this regard?


texts: scipy.signal.deconvolve not working, Blind deconvolution in python

texts:
I am trying to deconvolve my original signal with a given wavelet (Ricker wavelet). But scipy.signal.deconvolve is not giving any output to the signal.
Impulse signal can be generated by below code.

dt=1e-8
fs=1/dt
f0=2.5e6
t0=1/f0
nCyc=2
fmax=1/dt
fmin=0;
fnyq=fmax/2
rick=np.zeros((nt))

for it in range(nt):
  temp=math.pi*f0*(it*dt-t0)
  temp=temp*temp
  rick[it]=(1.0-2.0*temp)*math.exp(-temp)

plt.plot(rick)

Original Signal data is attached in this link
https://drive.google.com/file/d/1vHbv0WUY6jRFLFy_J_JJuL-TzdYVBgpe/view?usp=sharing

first column is time and second column is amplitude.

Kindly help how to deconvolve the signal. Is there any other way to do that.

is it possible to perform blind deconvolution of this signal?


texts: Determining the Nyquist rate of signals added and multiplied together

texts:
Suppose I have a signal in the frequency domain $x(\omega)$.

Suppose I know that its Nyquist rate is $\omega_0$ (and let $\omega_m$ be the its maximum frequency).

I need to get the Nyquist rate of

$$x(\omega)+x(\omega)e^{-j\omega}$$

The answers say it's still $\omega_0$.

Why is it not

$$\omega_m+\omega_me^{-j\omega}$$

and if we substitute $2\omega_m = \omega_0$

$$\frac{1}{2}\omega_0+\frac{1}{2}\omega_0e^{-j\omega}$$


texts: Generating data with Gaussian distribution

texts:
We want to generate a Gaussian vector, (call it $\Delta$) - with an arbitrary size-, with zero mean and variance of $\alpha$. If $\alpha$ is chosen such that $\parallel\Delta\parallel_2 \leq 0.5$ with probability of e.g. $90\%$.

How can we do that in Matlab?


texts:
I have an OFDM signal affected by AWGN, if I know the power of the gaussian noise, can I recover the OFDM signal? Can I perform deconvolution to recover the original signal?
thanks


texts: What is the proper way to select the time-bandwidth product when doing multi-taper spectral estimation?

texts:
Frequently programs for calculating the spectrum of a signal using the multi-taper method require the user to specify the time-bandwidth product.

What considerations must be made when selecting the time-bandwidth product?
Would would be the effect on the resulting spectrum if the selected time-bandwidth product is too large or too small?


texts: Algorithms for computing FFT in parallel

texts:
I am trying to parallelize the computation of an FFT on terabyte-sized signal files. Right now such an FFT using an open-source library takes many hours, even running through CUDA on the fastest GPU I have. The framework I am trying to adapt to this process is Hadoop. In very basic terms, Hadoop distributes a problem over any number of server nodes in the following manner:

  • You split your input file into (key, value) pairs.
  •   These pairs
  are fed into a “Map” algorithm, which transforms your (key, value)
  pairs into some other (key, value) pairs based on what you put inside
  the Map.
  •   The framework then collects all the (key, value)
  outputs from the Maps and sorts them by key, as well as aggregating
  values with the same key to a single pair,  so you end up with (key,
  list(value1, value2, ..)) pairs
   •  These  pairs are then
  fed into a “Reduce” algorithm, which in turn outputs more (key, value)
  pairs as your final result (written to a file).

There are many applications for this model in practical stuff like processing server logs, but I am having a hard time applying the framework to chopping up an FFT into “map” and “reduce” tasks, especially since I am not really familiar with DSP.

I will not bother you with the programming mumbo jumbo, as this is a DSP Q&A. I am, however,  confused on what algorithms exist for computing FFTs in parallel; Map and Reduce tasks can’t (technically) talk to each other, so the FFT must be split into independent problems from which the results can somehow be recombined at the end.

I have programmed a simple implementation of Cooley-Tukey Radix 2 DIT that works on small examples, but using it for recursively calculating odd/even index DFTs for a billion bytes will not work.  I have spent a few weeks on reading many papers, including one on a MapReduce FFT algorithm (written by Tsz-Wo Sze as part of his paper on SSA multiplication, I cannot link more than 2 hyperlinks) and the “four-step FFT” (here and here), which seem similar to each other and to what I am trying to accomplish. However, I am hopelessly bad at mathematics, and applying any of those methods by hand to a simple set of something like {1,2, 3, 4, 5, 6, 7, 8} (with all imaginary components being 0) gives me wildly incorrect results. Can anyone explain an efficient parallel FFT algorithm to me in plain English (one that  I linked or any other) so that I may try and program it?

Edit: Jim Clay and anyone else who may be confused by my explanation, I am trying to do a single FFT of the terabyte file. But I want to able to do it concurrently on multiple servers in order to speed up the process.


texts:
[![enter image description here][1]][1]

For this problem:
Actually,i think that the question asks us to get the IDFT from the DFT,so I did it..but I am not sure,till now, whether this problems asks us to do transformations or?


texts:
Is there a way to make an inverse notch filter block in matlab simulink?
I have found the peak-notch block, but I need to amplify the signal instead of attenuate it.

Thank you.


texts:
I'm using MATLAB's fft2 to generate the fourier power spectrum of an image. Naturally the spectrum is two dimensional as well, showing the direction of the direction and the frequency of the spectral components. I'd like to get rid of the directional component and just look at the frequency, i.e. reduce the dimension of the spectrum by one.

The way I see it, one can either:

sum over a ring containing all pixels with
$$
r_1 \leq \sqrt{x^2 + y^2} \leq r_2
$$
(which would make it tricky to select the correct boundary radii to really count the pixels that correspond to a certain frequency)
or one can convert the image to a polar coordinate system and sum over the axis corresponding to the radial component of the matrix.

I wonder which method would be preferrable and for what reasons and what normalization I would have to make (of course, for the first method, there would be much less pixels contained in the ring element for lower frequencies than for higher frequencies and I guess some similar scaling problem will occur for the second method)


texts: Regarding BIQUAD IIR implementation in MATLAB

texts:
I want to implement a high pass filter (order 4) in MATLAB using cascade form.

[b,a]=butter(4,0.0083,'high');% I am using butterworth filter
[sos,gain]=tf2sos(b,a);

Should I use [sos,gain] for implementing the same in cascade form or just [sos]. I read a post related to same and it says just [sos] should be used, however, I am not sure if that is the correct way.
P.S. I want to create a MATLAB function for IIR filter(don't want to use the standard filter function in MATLAB) which would get the coefficients b,a & input.


texts: How to compute fundamental frequency from a list of overtones?

texts: Given a list of overtones (F1, F2, F3, etc), how do I compute the fundamental frequency? Can I do something like F2/F1=F1/F0? Is it the correct method to use?


texts:
I'm new to signal processing and working my way through a textbook. There is an exercise where a causal LTI system is given that responds to a rectangular pulse.
I have an exercise where a causal LTI system responds for a rectangular pulse
$$ s(t) = \Pi(t) = \operatorname{rect}(t) $$
with a triangular pulse.
$$ g(t) = \Lambda(2t) = \operatorname{tri}(2t)$$
One task is to find out what the impulse response for this system is. The solution is
$$h(t) = \frac{d}{dt} h_{\epsilon}(t) $$
because $$\delta(t) = \frac{d}{dt} \epsilon(t)$$
where $\epsilon(t)$ is the unit step function, $\delta(t)$ the unit impulse function.
and looks like this

I don't really understand the solution and how it turned to an endlessly oscillating rect pulse. Unfortunately, the textbook doesn't explain its approach. So, why does the output signal looks like this?
EDIT:

Translation:
A causal LTI system responds to a rectangular pulse  $s(t) = \operatorname{rect}(t)$ with $g(t) = \Lambda(2t)$
a) What is the response for $s_1(t) = \operatorname{rect}(\frac{t-1}{2})$?
b) What is the unit step response?
c) What is the impulse response?
And the solutions for a) - c) from the textbook


texts: FIR Impulse response & transfer function

texts:
I am working though a dsp past paper and I have come across the questions "find the impulse response in the time domain" & "find the transfer function in the time domain" I know its really simple but the "time domain" part of the question catches me out, I know how to find the transfer function in the Z domain, its not for homework, its a past paper I am looking through.

How would I find the impulse response in the time domain?

AND

How would I find the transfer function in the time domain?

I have already tried YouTube and google and a few books but haven't found anything that is specific to those questions. The question is attached in the photo.


texts: Developing Lead Compensator to Decrease Settle Time

texts:
This is a follow up to this question.
I've approximated the transfer function for a system to be H = zpk([0.012 -1.05 18],[-0.22 -0.22 -45 -1000],10000);. My goal is to try to reduce the settle time to under 3 seconds. Thanks to answers on my previous question, I found the issue was that I have a slow, unstable zero at 0.012 - meaning I have a non-minimum phase system.
I've been looking into how to fix this issue, and found compensators could be used to counteract the slow zero. At first I implemented a parallel compensator, which did indeed help. I was able to get the settle time down to 20 seconds when previously it was above 50. From my research it seems that a lead compensator can help in reducing the rise time and settle time, but I've been having trouble with the implementation, as the math to determine the zeros and poles for the compensator are bit complex.
I believe it will be of the form $$C(s) = K_c\frac{s-z_0}{s-p_0}$$
where the zero, pole, and gain are determined, in part, using the root locus plot. If anyone could help me in determining an effective compensator, that would be greatly appreciated.


texts: Why bother impluse invariance, bilinear transform?

texts:
Now that we can operate in the frequency domain by way of FFT. We can implement ideal lowpass, highpass and bandpass digital filters. then why there many textbooks talk about designing digital Butterworth filter, digital Chebyshev filter etc, by the so called impluse invariance method, bilinear transform method?

I cannot warp my mind around this.

Thanks a lot!


texts: Use Fourier Transformation to compare 2 short audio files

texts:
I want to use Fourier Transformation to compute the similarity of 2 audios. I am dealing with small audios having the same size: 1 second length, mono channel.
My question is how can I use the Fourier Transformation of these 2 audio files to check wether they have the same sound or not.
My focus is about high intensity sounds like hand clap for instance. The 1st audio file is the reference one. It represents a hand clap sound, for example. And I am checking if the 2nd audio includes or not this hand clap sound.
Would be possible to describe how can I perform this check (I am programming in java) ?
Thanks in advance.
Kind regards.


texts: OFDM symbol acquisition block in GNU Radio gives rotating constellation even if there is no carrier offset. How is this possible?

texts:
Why does the constellation rotate when using OFDM Symbol Acquisition block in this flow graph below ?

The only explanation I could think of would be carrier frequency offset but in this flowgraph there is no carrier offset, so how can the constellation rotate if there is no carrier offset ?

I think the OFDM Symbol Acquisition block only finds where the OFDM symbol is (by finding the cyclic prefix using the autocorrelation function - I believe) and does not do any further magic.

I am asking this because I am trying to demodulate a recorded OFDM signal for which I do not know the Sync Word but I can do carrier synchronization "by hand", so the OFDM Symbol Acquisition block would be useful in this case.

Here is the GRC file : http://pastebin.com/2WvRG2Ly

EDIT:
Inspired by Marcus' answer I made the following simpler flow graph:

Which gives this result:

This is not spinning anymore :)

Now I wonder why is the constellation so noisy if there is no noisy channel?

It should not be noisy right ?

EDIT2:

It seems to be noisy because the source is noisy, strange. Why is that ?

EDIT3: The source noisiness comes from the RRC filter.


texts:
I working on FPGA block which requires an Sinc3 Filter for Sigma-Delta ADC coming in.
I read about it and Sinc filter, is filter with frequency response of sinc, that means its just a moving average filter right?
Frequency repsonse:

Time Domain:

so if I need Sinc3 filter implemented, which is 3rd order of Sinc filter, that means it's the first order Sinc filter cascades in serial 3 times, right?
So I Know there is an easy implementation of moving average filter using FIR filter, so to create Sinc3 filter all I need is to take the first FIR filter and cascade it 3 times in a row to create Snc3 filter, am I get it right?
Thanks


texts:
In their textbook(Example 1.1), Oppenheim and Willsky provide the following method for obtaining $g(t)=x(-t+1)$ from $x(t)$ : Replace $t$ with $-t$ in $x(t+1)$. In other words, advance the signal by 1 unit and flip the result.

But there could be another interpretation as well. Based on the fact that $-t+1 = -(t-1)$, we can have the desired signal as $h(t)=x(-t+1) = x(-(t-1))$,i.e. Delay the signal by 1 unit and flip the result.

But the two signals $g(t)$ and $h(t)$ are not identical. Why is the latter interpretation incorrect ?


texts: removal of blackman harris window results in signal higher than expected

texts:
We are capturing data using a PX14400 adc.  The samples are 16bit unsigned integer, 200M samples/sec, I data only.  100MHz bandwith.  We have a cw signal at 51MHz, -11.4dBm.  The -11.4dBm is measured by a power meter so that includes all losses right before being digitized.

Question 1:  when I apply a blackman harris 4 term window, the fft data calculates to -14dBm which sounds reasonable since bh loss is about 3dB, however when I dont apply the window the fft data calculates to -5.3dBm.  I would have expected a 3dB jump, not almost 9.  Why the large jump?

Question 2:  Does anything jump out in the code below that might be off in my calculations.  For brevity not all the code is shown (threads, board init, etc).

// Initialzied elsewhere in the class
const int fftSize = 8192;         // size of the fft
complex<double> array[fftSize];   // array of complex (r+i)...used post fft
double dData[fftSize];            // array of doubles, used pre and post fft
fftPlan = fftw_plan_dft_r2c_1d(fftSize, dData, reinterpret_cast<fftw_complex*>(array), FFTW_ESTIMATE);

void processNewBuff(uint16_t *rawdata) {

        // convert the samples to V.  This formula is provided by the vendor
        // dData is an array containing the converted to V data (double)
        // rawdata is an array containing the initial adc counts
        // since we only fft the first 8192 values, we effectively throw out the rest of the data
        // we skip the rest of the 10ms sample (200MS/s * .01 = 2MB sample, 8192 is processed the rest is discarded
        double dInputRange = .220;
        for (int x=0; x<fftSize; x++) {
            dData[x] = (-dInputRange)/2 + ((rawdata[x]/65532) * dInputRange);
        }

        // apply blackmanHarris
//        double bhtmp;
//        for (int x=0; x<fftSize; x++) {
//            bhtmp = (double)(2 * 3.1415 * x / fftSize);
//            dData[x] *=(double)(0.35875 - (0.48829*cos(bhtmp)) + (0.14128 * cos(2*bhtmp)) - (0.01168*cos(3*bhtmp)));
//        }

        // do the fft.  fftPlan is a "fftw_plan_dft_r2c_1d", discreet fourier transform, real 2 complex
        // see: http://www.fftw.org/doc/Real_002ddata-DFTs.html
        // after the fft, we throw away the upper half of the fft since we only started with 'I' data
        fftw_execute(fftPlan);

        // normalize the fft for power
        for (unsigned int x = 0; x <(fftSize/2); x++)
            array[x] = array[x] / complex<double>(fftSize, 0);

        // convert the fft to power spectrum
        // I have this separated into 3 distinct lines for readability, final code will be all in a single line
        for (unsigned int x = 0; x <(fftSize/2); x++) {
                dData[x] = sqrt(real( array[x] * conj(array[x]) ));
                dData[x] = ((dData[x]) / 50) * 1000;     // convert to W, then mW
                dData[x] = (20 * log10(dData[x]));
        }
}


texts: Given are two Signals $x_1[n]$ and $x_2[n]$. $x_1[n]$ is in Intervall $[0,2]$ different than null and $x_2[n]$ is in Intervall $[0,3]$. The convolution product is Null outside intervals: $[0,6]$ and $[-2,10]$. Can someone explain me why? I understand for 0,6 it's simply multiplication but why $[-2,10]$?


texts: Fast Fourier Inversion: Functions of a Complex Argument $f:\mathbb{C} \rightarrow \mathbb{R}$

texts:
I originally posted this question on math stack exchange, but I think it may be better suited for this community.
I'm interested in functions $f: \mathbb{C} \rightarrow \mathbb{R}$ with associated Fourier decompositions $$ f(a + ib) = \int_{-\infty}^{\infty} F(\lambda) \ e^{i \lambda (a + ib)} \ d\lambda.$$
We can assume these functions are band-limited in that there exists some $\mu \in \mathbb{R}$ such that $F(\lambda)$ is finite for $\lambda \leq \lvert \mu \rvert$ and $F(\lambda) = 0$ for $\lambda > \lvert \mu \rvert$.
I'd like to be able to express the above integral (or its discretization as a sum) in such a way that would allow me to perform an FFT (or other fast integral transform), though this has escaped me thus far. From what I've read, it seems like a Fourier-Mellin transform would be a good candidate (or maybe a clever application of the residue theorem), however I haven't been able to find this case in the literature.
Is it possible to evaluate this integral in a numerically efficient manner, and if so, how would I do so?
EDIT
Per the suggestions of Andy and Marcus, the problem I'd like to solve is really a discretized version of the above. Specifically, I have a discrete set of Fourier coefficients $$\{ {\bf F}_\lambda \}, \ \lambda \in \mathbb{Z}, \ -\frac{N}{2} \leq \lambda \leq \frac{N}{2}.$$ I know that these coefficients correspond to a function $f(z)$ defined on a sub-interval of the complex line such that $$ f(z) = \sum_{\lambda = -\frac{N}{2}}^{\frac{N}{2}} {\bf F}_{\lambda} \ e^{i \lambda z},$$ where $$ \Re(z) \in [0, \ 2\pi), \ \ \Im(z) \in \mathbb{R}.$$
My aim is to recover the values of $f$ on a discrete sampling of points in its domain.
One way to do this is as follows. Writing $z = a + ib$, the above sum can be written as $$ f(z) = \sum_{\lambda = -\frac{N}{2}}^{\frac{N}{2}} {\bf F}_n \ e^{i \lambda (a + ib)} = \sum_{\lambda = -\frac{N}{2}}^{\frac{N}{2}} \widetilde{\bf F}_{\lambda}(b) \ e^{i \lambda a},$$where $\widetilde{\bf F}_{\lambda}(b) = {\bf F}_{\lambda} \ e^{-\lambda b}.$ From here it's clear that for any $b \in \mathbb{R}$, we can recover the values of $f$ at the set of points $\left\{ \frac{\pi n}{N} + ib \in \mathbb{C} \ \vert \ 0 \leq n < 2N \right \}$, by running a 1D inverse FFT. Assuming we want to sample $b$ at $M$ discrete points, we’ll perform $M$ inverse FFTs, giving a total complexity of $O(MN \log N)$ in evaluating $f$ at the sample points.
I want to know if it's possible to do this more efficiently, though the odds seem long. In my application, almost anything would be fair game - assuming $f$ is periodic in $b$ on some finite sub-interval of $\mathbb{R}$, etc.


texts: Similar to this question: Visually stunning math concepts which are easy to explain, what are some great visualizations of basic DSP concepts such as FFTs, filters, etc?


texts: Extending Goertzel algorithm to 24 kHz, 32 kHz and 48 kHz

texts:
I'm learning to implement Goertzel's algorithm to detect DTMF tones from recorded wave files. I got one implemented in python from here. It supports audio sampled at 8 kHz and 16 kHz. I would like to extend it to support audio files sampled at 24 kHz, 32 kHz and 48 kHz.

From the code I got from the link above, I see that the author has set the following precondition parameters/constants:

self.MAX_BINS = 8
if pfreq == 16000:
    self.GOERTZEL_N = 210
    self.SAMPLING_RATE = 16000
else:
    self.GOERTZEL_N = 92
    self.SAMPLING_RATE = 8000

According to this article, before one can do the actual Goertzel, two of the preliminary calculations are:

Decide on the sampling rate.
Choose the block size, N

So, the author has clearly set block size as 210 for 16k sampled inputs and 92 for 8k sampled inputs. Now, I would like to understand:

how the author has arrived at this block size?
what would be the block size for 24k, 32k and 48k samples?


texts:
I am writing a synthesizer. One part of the synthesizer specification says that different sources can specify an attenuation to be applied to a signal, and that these attenuations should be summed.
The attenuations are specified in cB and are given as positive values. An attenuation of 0cB means no attenuation.
The naive approach, which I've tried, would be:
$$
x_{total} = 20 \log_{10}\left( \sum 10^{\frac{x}{20}}\right)
$$
for each $x$ as the attenuation (converted to dB) being summed.
This has problems, since if I add two 0cB attenuations (i.e. no attenuation, twice), I get a total attenuation of 60.2cB. This is obviously not what I want.
How can I add my attenuations to get a sensible result out the other side? Thank you.


texts:
I want to understand how SINAD is calculated. In MATLAB there is a funciton called SINAD and I would like to see if I can obtain "the same" answer as the inbuilt function does.
From the defintiion, I understand that it is the ratio of RMS of the signal to the noise and it takes all the values including the harmonics of the signal but the DC value.
So given the following code:
Fi = 2500;
Fs = 48e3;
N = 1024;
X = sin(2*pi*Fi/Fs*(1:N)); %Signal
Y= 0.001*randn(1,N); % Noise
x=X+Y;              %Signal + Noise
SINAD = sinad(x,Fs)

Matlab gives  me the value of SINAD for this particular signal. I believe that DC value is the first input to the signal, so if I wanted to calculate it "by hand" I would calculate the RMS of all the samples but the first one like this:
A=X([2:end]);
B=Y([2:end]);
z=sqrt(mean(A.^2));
w=sqrt(mean(B.^2));
sinad=20*(log10(z/w))

The answers are close, usually within 0.5 but are not the same. Is what I am doing correct? If not, how would I reproduce the SINAD function?


texts:
I am looking here at the course of current fluctuations and have calculated the RMS values, but somehow it seems strange to me: with the black line there are small fluctuations, but the RMS value is still comparably large with, for example, the blue...
The RMS values are:
red:0.83mV
blue:0.44
green:0.52
black: 0.31mV
compare with the picute


texts: What type of images are required for voxel classification?

texts:
In an algorithm, voxel classification approach is used to perform image segmentation. The aim was classifying voxels either as cartilage or background in MRI image. The project has to be completed in matlab.

MRI images has .dcm extension, it may have .avi and other types too, I am not sure.  Is it possible to perform voxel classification on .dcm images or it require another type of images?


texts:
For a digital lock-in amplifier (synchronous demodulator) I have a sinusoidal reference signal in the range of 10-15kHz, the amplitude and phase I want to observe changes very slowly <10Hz.

My question is do I gain any benefit by oversampling the input signal, or is twice the reference frequency enough, i.e ~30kHz in my case?

I ask because the A/D converted signal is mixed with the reference oscillator and low-passed filtered to almost DC.


texts: What is the difference between circular and linear convolution? When would I choose one over the other? In image processing where a filter is applied to an image with a mask which type of convolution should I choose?


texts: Mitigate floating-point numerical errors for very-low corner frequency low-pass filter with DSP

texts:
I'm designing a low-pass filter for a digital signal processing application that ideally just passes a very small bandwidth above DC. I'm using an IIR biquad filter for this, where the coefficients are derived using the instructions here. A smaller bandwidth leads to a longer filtering time (larger time constant) but yields a more accurate result whereas a larger bandwidth can be filtered faster but is less accurate. Both of these are valid use cases.
Here's the code I've got
#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz

# calculates filter coefficients using link above
# fc is corner frequency, fs is sample freq
def iir_lp_coeffs(fc, fs):
    w0 = 2 * np.pi * fc / fs
    q = 1 / np.sqrt(2)
    alpha = np.sin(w0) / (2 * q)
    b0 = (1 - np.cos(w0)) / 2
    b1 = 1 - np.cos(w0)
    b2 = b0
    a0 = 1 + alpha
    a1 = -2 * np.cos(w0)
    a2 = 1 - alpha
    b0 /= a0
    b1 /= a0
    b2 /= a0
    a1 /= a0
    a2 /= a0
    a0 /= a0
    return (
        np.array([b0, b1, b2], dtype=np.float64),
        np.array([a0, a1, a2], dtype=np.float64),
    )

fc = 2  # low pass corner frequency (Hz)
fsample = 500e3
b, a = iir_lp_coeffs(fc, fsample)

w, h = freqz(b, a, worN=int(1e6), fs=fsample)
fig, ax = plt.subplots()
ax.plot(w, 20 * np.log10(abs(h)))
ax.set_ylim(-40, 10)
ax.set_xscale("log")
plt.show()
print(w[0:10])
print(abs(h[0:10]))

The current settings use 64-bit floating point with a cutoff frequency of $2\,\text{Hz}$. This all works fine, and I can even decrease the corner frequency substantially as long as I increase the granularity of freqz (with worN=).
For instance here's a plot of the gain response with the above code (note that I've cut the x axis off at the higher frequencies):

However, my actual application requires 32-bit floating point. When I do this (set dtype of iir_lp_coeffs to np.float32), I get non-unity gain in the passband. For instance, here's a gain response with fc=10 using 32-bit:

If I set the corner frequency higher, the gain response looks correct again (e.g., fc=100 looks fine).
Am I running up against the limit of what's possible with 32-bit FP? Or, is there another strategy that would allow me to get away with the lower precision of 32-bit? Have I correctly diagnosed this issue as a floating-point issue?


texts: Gaussian Energy Broadening implementation to meet detector response

texts:
Respected Coding experts,

I am trying to perform Gaussian broadening on an energy spectrum. The final result should convert "specrum.png" to a similar form of "Exp.png". The method is explained in the "snap.jpg". In "snap.jpg" each bin is converted in a Gaussian where the "μ" is the mean, and sigma of the Gaussian is a polynomial function dependent on the energy i.e. x-axis (see code). Here each bin is converted to a Gaussian in a for loop and I want that the second bin gaussian includes a pileup from the Gaussian of the first bin. Similarly, the third Gaussian should include the pileup from second and first. The I am using for this is below:

E=linspace(0,10,1026);    % Energy bins
for i=1:1026
    sigma(i) = (-0.00789) + (0.06769)*E(i) + (0.021159)*(E(i)^2);                          % Sigma polynomial function
    gaussFilter(i) = (1/sqrt(2*pi* (sigma(i)^2)))*exp(-((E(i)-0.5)))/(2*(sigma(i)^2));   % Gaussian function
    gaussFilter(i) = gaussFilter(i) / sum (gaussFilter(i));
    AAA=W_NAIF/sum(W_NAIF);
    newSPECU(i) = conv(AAA(i),gaussFilter(i));
end
figure(1);
stairs(E,newSPECU);
set(gca,'YScale','log');
xlabel('Energy Deposited (MeV)');
ylabel('Counts');
title('Geant4 Gamma energy spectra (NaI NEAR - Tap Water)');
grid minor;
grid on;
ylim([1 1e6]);

The data for this plot can be downloaded at DATA


texts: Matching features points of rotate image using SURF

texts:
I'm trying to match SURF points between 2 images using SURF with opencv. One is the rotation of the other.
The problem is that the few matching it can found are wrong.

image of matching points

I mainly took my code from opencv http://docs.opencv.org/doc/tutorials/features2d/feature_flann_matcher/feature_flann_matcher.html#feature-flann-matcher

Here is the code:

string imageName1="test_right_rotate.jpg";
string imageName2="test_right.jpg";
Mat image1 = imread( imageName1, 1 );
Mat image2 = imread( imageName2, 1 );
int minHessian = 500;
SurfFeatureDetector detector( minHessian );
std::vector<KeyPoint> keypoints_1, keypoints_2;

detector->detect( image1, keypoints_1 );
detector->detect( image2, keypoints_2 );

SurfDescriptorExtractor extractor;
Mat descriptors_1, descriptors_2;

extractor.compute( image1, keypoints_1, descriptors_1 );
extractor.compute( image2, keypoints_2, descriptors_2 );

FlannBasedMatcher matcher;
std::vector<vector<DMatch > > matches;
matcher.knnMatch(descriptors_1,descriptors_2, matches, 2);
std::vector< DMatch > good_matches;
for(int i = 0; i < min(descriptors_1.rows-1,(int) matches.size()); i++) //THIS LOOP IS SENSITIVE TO SEGFAULTS
{
    if((matches[i][0].distance < 0.6*(matches[i][1].distance)) && ((int) matches[i].size()<=2 && (int) matches[i].size()>0))
    {
        good_matches.push_back(matches[i][0]);
    }
}
Mat img_matches;
drawMatches( image1, keypoints_1, image2, keypoints_2, good_matches, img_matches, Scalar::all(-1), Scalar::all(-1), vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS );

It works well for scaled images and translated images but fails for rotated images.

EDIT:

I don't know what is wrong but instead of using SURF, I used ORB as detector and FREAK as extractor and it works a lot better.
You can find the code on this SO : https://stackoverflow.com/questions/12491022/opencv-freak-fast-retina-keypoint-descriptor


texts: How would one design a (quasi) linear phase adaptive notch filter for a single complex tone?

texts:
While IIR notch filters are attractive, I need to retain phase linearity at the filter output.  I imagine that it's possible to use a standard IIR notch  filter:
https://www.researchgate.net/publication/224174145_A_Complex_Adaptive_Notch_Filter
followed by some form of phase equalization.
An alternative would be an FIR filter with conjugate symmetric coefficients.  I am unaware of any FIR designs in the literature explicitly for a single complex tone.

Any thoughts on the best way forward?


texts: How to define the scales in a biorthogonal, dyadic Discrete Wavelet Transform in Python

texts:
I am sorry if this has been discussed before but I could not find any answer.

I am trying to implement an algorithm that I found in a paper which asks for a biorthogonal, dyadic DWT with scales 1 and 2 using 'algorithme a trous'. I use PyWavelets' Stationary Wavelet Transform function to achieve this but I could not find any way to specify scales.

I would really appreciate if someone could help me about this either by using Pywavelets or any other library.

In below link's page 25, you can find the functions provided by PyWavelets.

Link


texts: Efficient computation of dot product of convolutions

texts:
I have two FIR filters $f_{1},f_{2}$ and two short time-slices of signals $s_{1},s_{2}$. I need to compute: $$\left<s_{1}*f_{1},\ s_{2}*f_{2}\right>$$

Which I compute using six FFTs: each $s_{i}*f_{i}$ is computed by pointwise multiplication of the (padded) FFT vector of $s_{i}$ by the corresponding FFT of $f_{i}$, then an IFFT is computed and the result is truncated (this truncation is the reason computing the dot product in the frequency domain does not work).

The $f_{i}$ are known ahead of time, so this takes four FFTs at each iteration of my code.

Is there a way to get away with computing less? If the conjugate transpose of $f_{1}$ could be composed on $f_{2}$ and a FIR filter was the result, it would be posible to manage with just 2 FFTs with something like: $$IFFT\left(FFT(s_{2})\cdot FFT(f_{2}^T\circ f_{1})\right)$$

(where $FFT(f_{2}^T\circ f_{1})$ is computed ahead of time,) but this doesn't quite seem to work out.

A team member of mine asked a similar question (on another stackexchange community) but did not get helpful replies.

edit (details about truncation in my code):

To compute $\left<s*f\right>$ (supposing that $s$ is of length 20000 and $f$ of length 4000, though these numbers are made up)  both vectors are zero-padded to length $20000+4000-1$, FFTs are multiplied pointwise, IFFT is computed on the product, and the samples in indices $\left[3999, 19998\right]$ are taken from the result (this represents a closed range of indices on a vector indexed from 0.)

This gives us all of the "meaningful" result computable from the slice of the signal: all dot products between $f$ and a shift of the signal-slice $s$ such that $f$ is shifted onto a completely-known part of $s$ (this is done over successive time slices of a signal as it is recorded.) It's important to notice the result is of length 16000 (the length of $s$ minus the length of $f$).


texts: Why is the raised cosine pulse for digital communications so popular?

texts: The raised cosine pulse is essentially a windowed sinc function, designed with some roll-off in mind. Why do we use this particular window, as opposed to other windows such as Hann, Hamming, Kaiser, etc? All of these have zero ISI when applied to a sinc function, and have tunable frequency characteristics just like the raised cosine. Is the raised cosine optimal for digital comms in any sense?


texts: Null autocorrelation function and stationary

texts:
I can show that a process $X(t)$ is Wide Sense stationary (WSS) by showing that $E[X(t)]$ is constant and that its autocorrelation function is in function of $\tau=t_1-t_2$, that is, $R_X(t+\tau,t)=R_X(\tau)$.

My question is:
If $E[X(t)]$ is constant and $R_X(t+\tau,t)=0$ can I say the process is WSS? What does $R_X(t+\tau,t)=0$ mean? Can I say $R_X(t+\tau,t)=0=R_X(\tau)$ and therefore is WSS?


texts: How to generate a beam pattern from a beamformer on 3D

texts:
given a 3D microphone array with $M$ microphones, $d(f,\Theta)\in \mathbb{R}^{f\times M}$ is the steering vector for a given direction $\Theta=(\varphi,\theta)$, inclanation and azimuth respectively.

I have generated both delay and sum $h_{DS}=\frac{d(f,\Theta)}{M}$ and minimum variance distortionless response $h_{MVDR}=\frac{\Phi_y^{-1}d(f,\Theta)}{d^H(f,\Theta)\Phi_y^{-1}d(f,\Theta)}$.

I have a function that generates a plot given a beam pattern - 3 column table corresponding to $(\varphi,\theta, A)$, where $A$ is the amplitude in that direction.

How do I generate such a beam pattern from each filter? Is there any value in generating a beam pattern from the steering vector regardless of the filter?


texts:
Most fft tutorials use sinusoidal signals for demonstration, which makes the user already know what the fft is supposed to give as an output, but what about real time signals where we don't have any previous knowledge except the sampled data with the following problems:

No previous knowledge
Variable sampling time
Unknown Noise

I have a signal which has the previous characteristics, so I used the following Matlab code to interpolate and select a proper period (visually):

clear all
load('ExperimentData.mat')
close all
% Because the sampling period is not constant, I have first to interpolate the data
Ts=min(diff(t));Fs=1/Ts;
t_new=t(1):Ts:t(end);
x_new=interp1(t,x,t_new);
%% Extracting one period and Applying ZERO PADDING
tStart=0.7086;[~,idxStart]=min(abs(t_new-tStart));
tEnd=0.7686;[~,idxEnd]=min(abs(t_new-tEnd));
x1p=x_new(idxStart:idxEnd);
t1p=t_new(idxStart:idxEnd);
F=figure;plot(t,x,'b');hold on
plot(t_new,x_new,'r:')
plot(t1p,x1p,'.-g');
ylabel('x');xlabel('time[sec]');legend('Original signal','Interpolated Signal','One Period')

The outcome is the following:

I then applied the FFT once on the complet signal an once on one period with ZERO padding using the following code:

X = fft(x_new-mean(x_new));
N = size(x_new,2);
f=Fs/2 * linspace(0,1,N/2+1);
figure;
AX(1)=subplot(211);semilogx(f,abs(X(1:N/2+1)));ylabel('Magnitude');xlabel('Frequency(Hz)');title('Entire signal Without Zero Padding')
AX(2)=subplot(212);semilogx(f,angle(X(1:N/2+1)));ylabel('Phase');xlabel('Frequency(Hz)')
linkaxes(AX,'x');axis  tight;
%Zero Padding
x_1PZP=[x1p-mean(x1p) zeros(1,100000)];
X_1PZP = fft(x_1PZP);
figure;
AX(1)=subplot(211);semilogx(f,abs(X_1PZP(1:N/2+1)));ylabel('Magnitude');xlabel('Frequency(Hz)');title('One Period with Zero Padding')
AX(2)=subplot(212);semilogx(f,angle(X_1PZP(1:N/2+1)));ylabel('Phase');xlabel('Frequency(Hz)')
linkaxes(AX,'x');axis  tight

The results are:

My wonderings are:

Does the interpolation of a signal alter the FFT ?
From the Time Signal, the FFT on the complet signal is closer to reality. but it still doesn't show all frequencies. Why?
I almost can't extract any information fro the phase plot. What am I missing?

The data can be obtained from this link.


texts: How does complex data of an FFT and IQ modulation differ?

texts:
If I have an arbitrary time domain signal (e.g. a simple cosine) i can sample it with an ADC and get a time discrete representation.
I can perform an FFT on the samples and get a complex output.
The phase of the signal is:

$$\phi=\arctan\frac\Im\Re$$

Using an IQ down modulation of the same signal gives complex data, too, with the phase being:

$$\phi=\arctan\frac QI$$

Is there a difference between these two and when is one preferred over the other?

EDIT

Thanks for the input!

@Dilip Sarwate: Here is what I did: set $\phi=60°$ (see your example)

The results are exactly what I expected. I can see the frequency of 1 Hz in the spectrum (plot 3) and the initial phase of $60°$ in the phase spectrum (plot 4). I get the same result when I use a non-complex function, e.g. $x(t)=cos(2\pi t+\phi)$.

I read everywhere that phase information is lost if I don't use IQ down-modulation, but that does not seem to be the case, as I can see the phase information in the plots.

Or am I missing something here?


texts: How to append two bandlimited signals and make the result bandlimited without modifying the first signal

texts:
I need to append (not add) a bandlimited signal F(b) to a bandlimited signal F(a) and keep the result bandlimited without modifying the part that corresponds to F(a).

Both are bandlimited to the same frequency (for example to 1/4 the Nyquist frequency). If F(a) is 100 samples long and F(b) is 100 samples long, the result F(c) is 200 samples in length and the first 100 samples of F(c) are identical to F(a).

An obvious "solution" is to filter the boundary so that the discontinuity can be made bandlimited, but in this case, I am restricted from modifying the samples in F(a).

How can I filter (or manipulate) F(b) such that the result of appending the two signals is also bandlimited to the same frequency and preserves the integrity of F(b) as much as possible?

Edit: For context, here's the problem I'm trying to solve. I am performing procedural terrain generation using 2D heightmaps and need to add a tile of terrain to a pre-existing landscape that may or may not have been procedurally generated. I'm looking for an algorithm that will "smooth" the borders of the tiles so that the tile boundaries are not noticeable.


texts:
I have been trying to restore an image that was blurred with a known Point Spreading Function and corrupted with noise using a kalman filter. I have looked at theory and have a basic understanding of how kalman filters work. But I can't find any good material on how I can use it to restore an image. I have tried a method where I considered that the image was the state to be predicted with no control signal and a observation model equal to the blur function with a known measurement noise i.e,

$$X(k+1) = X(k)$$
$$Z(k) = H*X(k) + R$$
where,

X(k) is the 2D image
Z(k) is the observed noisy blurred image
H is the PSF
R is the noise

here is the code is used,

I = im2double(imread('E:\Wallpapers\2.jpg'));
I = I(:,:,1);

LEN = 2;
THETA = 5;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');

SNR = 20;
noisy = awgn1(blurred,SNR);

noise = noisy - blurred;
R = cov(noise);
H = fft2(PSF,size(I,1),size(I,2));
K = 1e5*eye(size(I));
P = 1e5*eye(size(I));
X = zeros(size(I));
S = zeros(size(I));

for i = 1 : 1000
    S = R + H*P*H';
    K = P*H'/S;
    X = X + K*(noisy - H*X);
    P = (I - K*H)*P*(I - K*H)' + K*R*K';
end

MSE_noise = sum(sum((noisy-I).^2));
MSE1_restored = sum(sum((X-I).^2));

The image is filled with NaNs which is because of K = P*H'/S being badly scaled.
Where am I going wrong? Is there a problem with the code or am I supposed to change the model to restore the image?


texts: Orders of complimentary butterworth filters for IMU

texts:
I'm currently trying to design a complimentary filter for deployment on a drone, as part of an assessed semester project. Even though processing power is not an issue (we have GHz, lots of pipelining and all the good stuff), I am still electing to use an IIR filter in the final implementation, as this is easier to represent mathematically in my opinion.

As I have read Here in this answer, certain orders of filters do not produce monotonic unity gain in output. As such, at this stage I am mostly interested in proving, in a round-about way, which orders of filters produce monotonic unity gain. Phase difference is almost inconsequential at this stage, as I am reasonably confident in not using such a high-order filter that we hit any timing constraints.

So, I have written some initial matlab code to test my application on a first order Butterworth pair, as seen below. It is worth noting that our sampling frequency is 1 Khz, meaning the nyquist frequency is 500 Hz. As such the normalized cutoff frequency of 0.25 equates to 125 Hz, although this is arbitrary at this stage. Possible problems include the summing of the 2 filters in the last line, or the need to add some integrator term to the high-pass. We also currently have some stock complimentary filters that came with the drone (all came as a big package from Quanser, their autonomous vehicle research suite thingy), but it seems lackluster to just use them and not know what they do. I will also be taking a closer look at the filters for the IMU and magnetometer in a short while after this is complete, as I think I can do a little better.

clc;
clear;
close all;
w_n=0.25; %cutoff frequency

[zL1,pL1,kL1] = butter(1,w_n,'low');
[numL1,denL1] = zp2tf(zL1,pL1,kL1);
sysL1=tf(numL1,denL1);

[zH1,pH1,kH1] = butter(1,w_n,'high');
[numH1,denH1] = zp2tf(zH1,pH1,kH1);
sysH1=tf(numH1,denH1);

sys1=sysL1+sysH1;

Making the bode plot of this transfer function gives me something completely different from what I was expecting, as also seen below.

Any help and guidance is appreciated, many thanks!


texts: Choosing the best measure for machine learning on a Spectrogram

texts:
I am working on an application to classify Talas(a piece of music with onsets). For that, I have the spectrogram of a Tala, which I want to compare with the test Tala, or compare 2 Talas.

For this I intend to find the onsets in the spectrogram and then form a classifier based on the frequencies that I get at the onset points in the Tala.

Now as we already know that a spectrogram contains multiple frequencies for a time t , I wish to quantify this as a single value to be used in machine learning, so that I have n values for the n onsets(say), and not multiple frequency values for the same time t.

So here is my question - Given all the frequencies at time t, how can I get a value that represents all these frequencies(like mean, median, etc). What should be the measure I should go ahead with?


texts: How to Generate Band Limited Gaussian White Noise in MATLAB?

texts:
In some papers, I read that the additive noise is band limited Gaussian white.
How can I simulate this type of noise use MATLAB?


texts:
I am currently looking over filters and understand the frequency response of a system, e.g. a low-pass filter. However I am confused by the 'Power/Amplitude response' of a system as detailed in my notes below:

  Square of the amplitude response, or power response:

  $$|H(\omega)|^2=H(\omega)H^*(\omega)=H(\omega)H(-\omega) $$

I created a basic matlab scenario for the a low-pass filter and plotted the frequency response and the power response:

From the top graph I understand that as the frequency increases along the x-axis you get an attenuation for difference frequency components in db.
However I do not understand what the significance or intuitive value is gained by squaring the attenuation (or gain). I can see that the gain becomes more negative as frequency increases and the amplitude response increases.
So for example my 3dB cut-off point is at 1rad/s, 10^(0). This has a corresponding amplitude response value of 9. Does this mean the amplitude has decreased by a ninth? Because surely if my gain were positive my amplitude response value would also be a positive number and so the opposite cannot be true?

Any help would be greatly appreciated!

Regards

Pat


texts: Filter response to sinusoidal input signal

texts:
This question is related to Does 46 dB gain of a filter for a frequency imply 200 times more amplitude? In  the accepted answer it is mentioned that "For a sine wave input the amplitude of the steady state output is simply the amplitude of the input multiplied with the magnitude of the transfer function at that frequency"

Above relation is derived as an integral; as a convolution between filter and the input. But in practice we have all discrete signal, therefore can we get 200 times more in practice?


texts: FIR filter design for response inside unit circle

texts:
I would like to design an FIR filter such that its Z-transform has a certain profile in certain regions.

For example, if I'd like to have an FIR filter that nulls the decaying exponential $\left(\frac{1}{2} \right)^n$, I can simply put a zero at $z = \frac{1}{2}$.

My question is, is there a more systematic way of designing such filter for arbitrary Z-transform response? Something similar to the Remez algorithm or maybe using convex optimization? So far all the papers I've looked at online looks at responses on the unit circle only.


texts: Selecting the Nyquist rate for a combination of signals

texts:
If several sine waves were combined but each had a separate frequency, for example 10 Hz, 20 Hz and and 30 Hz, what would the required Nyquist rate be for analysis?
From my research, it suggests that Nyquist rate should be double the maximum frequency and for this case it should be 60 Hz. Is this correct or does the fact that it is a combined signal change this concept?


texts: What tools do you use to create signal processing plots?

texts:
I am in process of creating some notes in Latex on signal processing.

such as

Does anybody know if there are softwares or tools that are useful for quickly generating these diagrams?


texts: IIR low pass filter: floating point and fixed point

texts:
I have implemented an IIR low pass filter in direct form 2, the input to the filter is a unit impulse $\delta[n]$. I have calculated the output of the filter for the mentioned input using the difference equation of direct form 2, the equations are given below:

$$\begin{align}w[n] &= x[n] - a_1 w[n-1] - a_2 w[n-2] \\
y[n] &= b_0 w[n] + b_1 w[n-1] + b_2 w[n-2] \end{align}$$

How to plot the frequency response of the output with out using inbuilt commands in MATLAB? Also I want to implement the same in fixed point using 32 bit and plot the frequency reponse of it.
For plotting should I use freqz command or plot(abs(fft(my ouput))?. Both should give me same frequency response, but I am getting different responses,


texts: Expand and shrink waveform visualization depending on tempo changes

texts:
I am not quite sure if this is the right state exchange site or not. Pls let me know.
So, I noticed that Logic and Pro Tools, etc. have this waveform shrinking and expanding effect when dragging the region through tempo changes.
Here's a short video: https://imgur.com/a/ol8JB6J
I was wondering how one would achieve that?
Given that you have a point cloud representing your peaks and a window size in which the waveform should be rendered, how could one shrink and expand the waveform in real time? Because the window size obviously changes depending on the tempi given during the window.
I guess you'd have to know the viewport pixelToSampleRateFrame ratio for every vertical pixel line, but still.
It seems super expensive to redraw it with every movement the user does.
I can't really wrap my head around how that would work in a efficient way and I thought that maybe someone here has an idea.