texts: Pesq numbers and clean speech comparison to processed speech

texts:
I am trying to evaluate my processed speech for my DSP project. My project is an implementation of an MMSE algorithm proposed by Y.Ephraim and D. Malah, "Speech Enhancement Using a Minimum Mean Square Error Short-Time Spectral Amplitude Estimator," in the IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-32, No. 6, Dec. 1984.
As per project requirements, I ran PESQ, LLR, and segment SNR. I got numbers, but I am not sure what they mean. Is high number is good or bad? Also, what does it mean if I get a negative number?

Please help and thanks in advance


texts: If noise is your signal, what is your noise?

texts:
Consider the following contrived situation. Imagine a Gaussian white noise process $x[t]$, with bandwidth $Δf$, with PSD equal to some quantity $A$ which you would like to measure.

So the way to measure this seems to measure the variance of the process $x[t]$, which by Parseval's theorem will be $AΔf$.

So you measure points with some frequency $f_m$, probably $2Δf$. And at a lower rate, say $f_v$, you compute the variance of the preceding block of $f_m/f_v$ points and take that to be a measurement. What will the noise/variance be in this measurement of the variance of $x[t]$? How can I approach this question?


texts: Compare feature vectors of an image in MATLAB

texts:
I am working on OCR. I have extracted features of 26 alphabets, and I want to compare them with an experimental character feature.

I have tried using Euclidean distance and correlation functions, but I am not getting good efficiency. What are some alternate methods for the same?


texts: Are the SNR, Data Rate and Delay values correct in the below calculation?

texts:
I have these values:
power_received_from_one_subcarrier = -80 dBm
bandwidth_for_one_subcarrier = 180 KHz
noise = -109 dBm

When using {/, *}, the dBm values are converted to Watts using the formula:
value_in_watts = 10 ** ((value_in_dBm - 30) / 10)

When using {+, -}, the Watts values are converted to dBm using the formula:
value_in_dBm = 10 * math.log10(value_in_watts) + 30

The SNR formula is:
snr = power_received_from_one_subcarrier / (noise * bandwidth_for_one_subcarrier)

Which gives the value:
snr = 0.00387

The data rate's formula is:
data_rate = bandwidth_for_one_subcarrier * log_base_2(1 + snr) = 180e3 * log(1 + 0.00387)

Which gives the value:
data_rate = 1002.26 bps

I do not know if my calculations are correct, hoping for some guidance.


texts:
I am currently trying to train a CNN-HMM acoustic model for speech recognition.

The CNN model is able to detect a center monophone given a context window of x (limits has not been tested yet - but works with 50) frames from a spectogram. The CNN provides me with a posterior probabilities of all possible monophones for the center monophone, but I am not sure how i am supposed to combine it with a transistion model such as HMM (Hidden markov model), given the CNN (Convolutional neural network) already provides posterior probability of all possible phones.

How am I supposed to train a HMM? Since i am classifying monophone, i seem to have a hard time understanding where the use of HMM would be appropriate here, as each HMM would just have one state, and encode one monophone, which the posterior probability in itself provides.

Is it possible to combine the CNN and HMM in a meaning full way, I am using kaldi, and the dataset consist of utterances of yesno from one speaker. (Simple case).


texts: How do I analyze a music file ( only contains 5 seconds of note A) using the FFT plot to determine the note being played?

texts:
To perform this task why can't I just import the file, and directly perform the fft and use the stem plot to determine the note? From the index k of the peak I can do w = 2pi*k/N where N is total number of points.

Fs, x = wavfile.read('piano.wav')
X = np.fft.fft(x)
plt.plot(abs(X))
plt.show()

However I'm not sure about the correct way to map the fourier coefficients from the plot to the frequency to determine the note.

Please advise.
Thanks.


texts: How to add Rayleigh Fading Coefficient to a BPSK signal with AWGN?

texts:
I have a problem with my code. I saw many programs that add the coefficients of rayleigh like this:
h=1/sqrt(2)*(randn(1,N)+1i*randn(1,N)), where N=number of bits.
For 1 user I think this code work great, but I´m simulating a BER for BPSK in Rayleigh Channel for 4 users, and I tried with the above code but It doesn't work. When you compare gaussian noise with coefficiente Rayleigh both are practically the same. So I tried:
SNR = 0:40;
P=1;
a1 = 0.75; a2 = 0.18; a3=0.05; a4=0.02;
ber1=zeros(1,length(SNR));
ber2=zeros(1,length(SNR));
ber3=zeros(1,length(SNR));
ber4=zeros(1,length(SNR));
for u = 1:length(SNR)
    %Generate random binary data
    data1 = randi([0 1],1,N);  %Data bits of user 1
    data2 = randi([0 1],1,N);  %Data bits of user 2
    data3 = randi([0 1],1,N);  %Data bits of user 3
    data4 = randi([0 1],1,N);  %Data bits of user 4
    %Do BPSK modulation of data
    x1 = 2*data1 - 1;
    x2 = 2*data2 - 1;
    x3 = 2*data3 - 1;
    x4 = 2*data4 - 1;
    %Rayleigh coefficent |h1|^2<|h2|^2...<|h4|^2
    h=zeros(2,4);
    h(1,:)=raylrnd(1,1,4);
    h(2,:)=h(1,:).^2;
    h=h';
    H = sortrows(h,2);
    H=H';
    %Superposition
    x = sqrt(P)*(sqrt(a1)*x1 + sqrt(a2)*x2 + sqrt(a3)*x3 + sqrt(a4)*x4);
    N0 = 1/10^(SNR(u)/10);
    %Received signals
    y1 = H(1,1).*x + sqrt(N0/2)*(randn(1,N));
    y2 = H(1,2).*x + sqrt(N0/2)*(randn(1,N));
    y3 = H(1,3).*x + sqrt(N0/2)*(randn(1,N));
    y4 = H(1,4).*x + sqrt(N0/2)*(randn(1,N));
    %Equalize
    eq1 = y1.*H(1,1);
    eq2 = y2.*H(1,2);
    eq3 = y3.*H(1,3);
    eq4 = y4.*H(1,4);
    %AT USER 1--------------------
    %Direct decoding of x1 from y1
    x1_hat = zeros(1,N);
    x1_hat(eq1>0) = 1;
    %AT USER 2-------------------------
    %Direct decoding of x1 from y2
    x12_hat = ones(1,N);
    x12_hat(eq2<0) = -1; %here
    y22_dash = eq2 - sqrt(a1*P)*x12_hat;
    x2_hat = zeros(1,N);
    x2_hat(y22_dash>0) = 1;
    %AT USER 3----------------------------
    x13_hat = ones(1,N);
    x13_hat(eq3 < 0) = -1; %Estimate user 1's signal first real
    y23_dash = eq3 - sqrt(a1*P)*x13_hat; %Extract useless signal 1 from y3
    x23_hat = ones(1,N);
    x23_hat(y23_dash<0) = -1;%here
    y33_dash=y23_dash - sqrt(a2*P)*x23_hat; %Extract useless signal 2 from y3
    %Decode x3 from rem
    x3_hat = zeros(1,N);
    x3_hat(y33_dash>0) = 1;%Final bits for user 3
    %AT USER 4-----------------------------------------
    x14_hat = ones(1,N);
    x14_hat(eq4 < 0) = -1; %Estimate user 1's signal first
    y24_dash = eq4 - sqrt(a1*P)*x14_hat; %Extract useless signal 1 from y4
    x24_hat = ones(1,N);
    x24_hat(y24_dash<0) = -1;%here
    y34_dash=y24_dash - sqrt(a2*P)*x24_hat;%Extract useless signal 2 from y4
    x34_hat = ones(1,N);
    x34_hat(y34_dash<0) = -1;%here
    y44_dash = y34_dash - sqrt(a3*P)*x34_hat;%Extract useless signal 3 from y4
    x4_hat = zeros(1,N);
    x4_hat(y44_dash>0) = 1;%Final bits for user 4
    ber1(u) = biterr(data1,x1_hat)/N;
    ber2(u) = biterr(data2,x2_hat)/N;
    ber3(u) = biterr(data3,x3_hat)/N;
    ber4(u) = biterr(data4,x4_hat)/N;

As you can see, I'm trying use raylrnd (Rayleigh random numbers) but it doesn't work as I expected. I hope you can help me, I've been a long time with this doubt/problem and I can't find a solution.


texts: Negative exponential signal's energy and power

texts: Earlier I have dealt with exponential functions multiplied with unit step function. But, energy and power of exponential function alone comes out to be infinite when I put limits of the integrals to from $-\infty$ to $+\infty$. How can I find its energy or power if the signal is not multiplied with unit step function?


texts:
In processing analog signals using discrete-time systems, it is generally desirable to
minimize the sampling rate. This is because the amount of arithmetic processing required
to implement the system is proportional to the number of samples to be processed.
If the input is not bandlimited or if the Nyquist frequency of the input is too high,
prefiltering maybe necessary. An example of such a situation occurs in processing speech
signals, where often only the low-frequency band up to about 3 to 4 kHz is required for
intelligibility, even though the speech signal may have significant frequency content in
the 4 kHz to 20 kHz range. Also, even if the signal is naturally bandlimited, wideband
additive noise may fill in the higher frequency range, and as a result of sampling, these
noise components would be aliased into the low-frequency band.
Could someone explain please the point below and why it is a problem?
Ie, why is it a problem if noise components are aliased into the low-frequency band?

"Also, even if the signal is naturally bandlimited, wideband
additive noise may fill in the higher frequency range, and as a result of sampling, these
noise components would be aliased into the low-frequency band."


texts: Relationship between FFT amplitude and sample size

texts:
Would anyone know if there is any relationship between the FFT amplitude vs. data length?

For instance, in Matlab,
if signal_30k had 30,000 data points,
sampled at 1000 Hz,
and I do this:

fft_amplitudes = abs( fft( signal_30k ) );

then dominant frequency's resulting amplitude is some number.

but the same signal with the same frequency components,
signal_500k had 500,000 data points,
also sampled at 1000 Hz,
and had the same operation done on it:

fft_amplitudes = abs( fft( signal_500k ) );

the same dominant frequency's resulting amplitude would be another number but much larger.

Is there any relationship between the number of data points and the resulting amplitude that comes out of the FFT?
The goal is to ensure that regardless of the number of samples, the amplitudes have some means of being compared with one another.

Any help would be greatly appreciated.


texts: Conceptual question on FFT and chirp signal

texts:
If I take the FFT of a sinusoid I will get a plot whit all the energy of the signal concentrated at the sinusoid frequency. But what happens if I have a signal in which the frequency keeps changing?(like the chirp: linear and exponential)

Is it possible to compute the FFT of a chirp signal? And what should I expect as result?

I tried to do it in Matlab creating a chirp using the matlab function and then passing it to the FFT function but I do not know if this makes any sense and how to 'read ' the result I get.

Could you also suggest me some readings (papers, books) on the topic?


texts: Compare two words from different speakers

texts:
I have a set of audio files which contain a word each one, and I'd want to be able to compare a new word (coming from a smartphone mic) to check if it resembles to a specific file and how much it resembles to that audio file. I don't need to know which word it is, or doing a STT system, actually, just checking how similar two words are. The words may come from different speakers, male, female, old people, kids, even non-speakers for a specific language.

Now, the question is, which should be the processes to do something quite simple? I've been reading about MFCC, DTW and I can ensure I understand little to nothing about how I should proceed. I'd need to create an algorithm to match those words, but I have not a single clue on how to proceed after the "Conception steps".

I have, however, tried an Android library called Musicg, which should compare sounds based on:

A frequency range for the whole sound
An intensity range for the whole sound
A standard deviation range (yes, for the whole sound)
A High pass and low pass limits
Range of times the sound crosses the "zero" line

But I am pretty sure this might not be enough to make sure two sounds are similar or different. The examples they handle are only onomatopoeii, (for those who code, here's an example to detect Claps).

Now, the question. Is there a magic algorithm that does all the heavy work of comparing two audios from different speakers that is able to give a similarity score? If so, can you put me back to the tracks with a simple (well, as simple as it can be, of course) example, please?

Bonus track: Is any of the algorithms able to do the processing with precomputed values from the set of audio files?


texts: Write a 440 Hz Sine Wave to .wav-File using Python and SciPy

texts:
This is my first post on DSP.SE! I'm a student in Computer Sciences and am just getting started with Signals and Systems, thus getting in touch with discrete signals and analyzing them.

I wanted to explore this world a little further and so my first goal is to create a Sine Wave that has a Frequency of 440 Hz and write it to a .wav-File. Of course i have all the calculus in mind, but i can't seem to figure out the problem myself!

So far, i have this:

#!/usr/bin/env python3

import numpy as np
from scipy.io import wavfile

sampleRate = 44100

t = np.linspace(0, 20, sampleRate * 5)  #  Produces a 5 second Audio-File
y = np.sin(440 * t)  #  Should have frequency of 440Hz

wavfile.write('Sine.wav', sampleRate, y)

I understand why this code always creates a File, that is 5 seconds long, but i can't figure out, how to get the frequency right!

Have i even gotten the rest of the code right, Dsp-Guy-wise?


texts:
Its true that HTC one M8 uses dual camera to achieve the 3D Dimension plus effect, but can anybody tell what is it called, if its a series of image manipulation?

Demonstrated here,
Video


texts: How to calculate the total signal strength and phase shift of multiple signals in a reciever

texts:
I have a simulation to make where I have an array of transmitters that transmit the same signal. At a random point, which I have to consider as a receiver I have to measure the phase shift of the signals and also measure the total signal strength. The requirement is that the transmitter emits the same signal but with different phase offsets.

As far I understood because the transmitters are separated apart, the received signal is already offset by some phase because they different travel distances. Is my assumption right?
I know the formula to find the phase offset given the distance travelled of two signals - how do I calculate it for multiple signals? Or am I understanding the question wrong?
How do I calculate the total signal strength for all signals combined? I know for one signal assuming a free space, line of sight connection I can use friis equation. But how do I use it for multiple signals?

Any help is very much appreciated. I am Software Engineer who have recently taken up Wireless communications course, and I'm not sure if I understand the motivation and meaning behind the simulation right. Thanks!


texts: What is the mean power of a complex random variable?

texts:
Say $\alpha$ is a complex random variable, then which one of the following expressions is correct?

$\mathbb{E}[\alpha^2]$ or
$\mathbb{E}[\alpha \alpha^*]$?


texts:
I have a basic question about stochastic processes:

When some informations such as wss, uncorraleted sampled, white about random signal (say x[n]) are given, what do we exactly have?

For example white information:
I know if x[n] is white, its power density function $ P_{x}(e^{j\omega}) $ is constant value independent of $ \omega $. I suppose it means autocorrelation function of x[n] consists of dirac? (I am not sure here, If I am wrong, please correct me).
Uncorraleted samples: I know if the random signal is uncorrelated, its covariance equals to zero which means its autocorrelation equals to square of its mean. But uncorraleted samples means this, or something else?
What about wss?: I know if the signal is wss, mean (average) of x[n] is constant which is useful information when $ x_{c}(t) $ (constinuous time form of x[n]) and its means are given and we can say directly mean of x[n] is the same as mean of $ x_{c}(t) $. Moreover, we have $ \phi_{x}[m]  =  \phi_{x_{c}}[mT] $ where T is the sampling period and $ \mu_{x_{c}}^{2} = \mu_{x}^{2} $ where $ \mu $ is the power of the signals (second moment).

Any correction, addition, comment would be appreciated.

Pass through an example may be more convenient:

Suppose we have a system like this:

where, x[n]: wss, $ \eta_{x} = 0 $ (mean value), uncorrelated samples, $ \sigma_{x}^{2} $ and h[n] are given.

How to calculate $ \eta_{y}[m] $ the mean, $ \phi_{y}[m] $ the autocorrelation and $ \mu_{y}[m] $ the average power of y[n]?
How to calculate $ \eta_{v}[m] $ the mean, $ \phi_{v}[m] $ the autocorrelation and $ \mu_{v}[m] $ the average power of v[n]?


texts: Simulation and theory of first-order low pass filter

texts:
I am trying to teach myself some signal processing and control theory. I'm currently working through THIS tutorial. I've started an ipython notebook to keep my notes in.

Unfortunately, something is wrong with my very first simulation of a simple first-order low pass filter. The differential equation is
$$\dot{x} = \omega_0(u-x)$$
where $u(t)$ is a sinusoidal input. My understanding is that this system is described by the transfer function:
$$G(s) = \frac{x(s)}{u(s)} = \frac{1}{1+s/\omega_0}$$
I am interested in using this transfer function to calculate the amplitude and phase lag of $x(t)$. To do this, we evaluate $G$ at $s=i\omega$:
$$|G(i\omega)| = \frac{|x(t)|}{|u(t)|} = \frac{1}{\sqrt{1+\omega^2/\omega_0^2}}$$
$$\arg G(i\omega) = \text{phase lag} = -\tan^{-1}\frac{\omega}{\omega_0}$$
Are the above equations correct? Is my interpretation of them correct? Because when I numerically simulate the above differential equation in python, I get results that do not agree. Particularly, the blue dots and red dots should line up with the black lines in my plots. In other words, I have been unable to reproduce Fig. 3 in this tutorial

The ode is so simple, I don't see how I could have messed up the simulation. You can see my code in the links above. Thanks in advance for any comments/suggestions. (P.S. I hope I picked the right stackexchange site for this question, let me know if I got it wrong).

Update: I have tried decreasing the time step for the numerical integration routine, but this does not improve the fit. The code is still available at the same url: [Click here for code]


texts: Detecting similar behaviour in audio files

texts:
Can anyone help, I'm trying to find a solution that would allow me the ability to look at a batch of .WAV files and detect ones which share a similar behaviour in the recording like the image attached.

Recording -> Large Gap -> Recording

Thanks!


texts: Why every musical instrument sounds different than any other while listening?

texts:
This is soft question.I am not expert in music. We know that English language has 26 letters (which are again grouped into vowels and consonants) , while language of music is same all over the world and it has only 7 letters (all are vowels). There are number of musical instruments. one can blindly recognize which instrument is being played just by listening to it.

Each musical instrument sounds (i am not finding proper word here in case of "sounds" ) different than any other musical instrument.Every instrument seems to be unique while listening.My basic question is why is it so?

Also,is it possible to replace one musical instrument by combining two or more other instruments playing together or like that ? (For example is it possible to obtain  same Guitar sound (i.e.sound/tones our ear listen while Guitar is being played) by  combination of other musical instruments say violin and harmonium etc played together ?)


texts:
$B(z)+B(-z) = 2c$, explain the structure of $b[n]$ and find the constraint of its length given that $c$ cannot be $0$.

This is a homework problem. "Explain the structure" means that $b[n]$ is zero for certain values of $n$, and has a certain shape.

I'm trying to take the inverse $\mathcal Z$-transform of $B(z)$ and $B(-z)$, but I'm not sure how the inverse $\mathcal Z$-transform of $B(-z)$ is related to $b[n]$, so I'm stuck...can anyone give me some advice?


texts: Auto correlation correspondence to original time domain window

texts:
I'm trying to figure out how exactly an auto correlation corresponds to a time domain signal.

Now I'm trying to find pitch periods in an audio file and 99% of the time I'm getting it spot on.  Unfortunately however that 1% of the time is throwing me some annoying issues.

Now at present I scan through my window for the highest peak (ie the highest abs( sample value )).  This seems to match up really nicely in most cases.

Now when I auto correlate I would expect the central peak to be this peak that I have identified (It will always have a positive number if i read it right).  Is this correct?

Secondly, A peak that is n samples away from the central peak of the auto correlation will be n samples away from the peak I've identified in time domain, ok?

Sooo, if this is correct then I'm having a strange problem.  I'm finding a peak in the auto correlation 91 samples from the center of my auto correlation.  Unfortunately when I look at the signal in the time domain (inside an audio editing app) then my central peak corresponds perfectly to a negative peak in the audio.  However, there are no peaks 91 samples away from that peak.

If I look forward the next real peak is 142 samples away and the previous peak to it is 75 samples away.  I do notice a very tiny peak 91 samples or so in front of my peak but the sample never even goes negative.  So why does this end up showing up as by far the strongest peak in the auto correlation?

Any help with this would be hugely appreciated!


texts:
i got similar question as in FFT scaling in GNU Radio
and quite not clear about the answer given so far. I understood that the FFT plot is a relative value in which we have to use measurement instrument to get the actual value.

But, without the measurement instrument, in general, what value will be shown on the plot if the given input is for example complex cosine waveform with amplitude 1?

I've tried the following flowgraph: signal_source - throttle - qt_frequency_sink

and vary the cosine waveform amplitude from 0 to 100. other than amplitude, all blocks use default setting. I expect to get 0dB when the amplitude is 1, but the plot shows ~ -10dB instead,and ~32dB when amplitude is 100 (i'm not sure if this value make sense). SO, how actually the frequency sink calculates these value?


texts: Deriving the DFT magnitude of $A\cos(2\pi nk/N)$

texts: Given that $$x(n) = A\cos(2\pi nk/N),$$ the $N$-point DFT of $x(n)$ can be expressed as follows—the derivation can be found in here: $$X(m) = \color{red}{\frac{A}{2}\sum_{n=1}^{N-1}e^{-j2\pi n(m-k)/N}}+\frac{A}{2}\sum_{n=1}^{N-1}e^{-j2\pi n(m+k)/N}.$$ It is said that the left summation above represents a positive-frequency spectral component. Does it mean to say that, when $m$ is negative, the summation goes down to zero?


texts:
I have heard anecdotaly that sampling complex signals need not follow Nyquist sampling rates but can actually be gotten away with half Nyquist sampling rates. I am wondering if there is any truth to this?

From Nyquist, we know that to unambiguously sample a signal, we need to sample at least higher than double the bandwidth of that signal. (I am defining bandwidth here as they do in the wiki link, aka, the occupancy of the positive frequency). In other words, if my signal exists from -B to B, I need to sample at least > 2*B to satisfy nyquist. If I mixed this signal up to fc, and wished to do bandpass sampling, I would need to sample at least > 4*B.

This is all great for real signals.

My question is, is there any truth to the statement that a complex baseband signal (aka, one that only exists on one side of the frequency spectrum) need not be sampled at a rate of at least > 2*B, but can in fact be adequately sampled at a rate of at least > B?

(I tend to think that if this is the case this is simply semantics, because you still have to take two samples (one real and one imaginary) per sample time in order to completely represent the rotating phasor, thereby strictly still following Nyquist...)

What are your thoughts?


texts:
Im trying to generate a sweep sine wave that increases its frequency till an upper limit and then starts decreasing its frequency to the lower limit at the same rate. Im trying to simulate in matlab.

The frequency of the wave will change based on the following formula: f = fmax-a*abs(T - t); a in the formula equals (fmax-fmin)/T
So the frequency should first increase and then decrease after it reaches time T at the same rate.

The matlab code is:

freq1 = 20;
freq2 = 200;
fs = 44100;
endTime = 0.1;
t = 1/fs:1/fs:(endTime*2);
a = (freq2-freq1)/endTime;
f = freq2 -   a*abs(endTime - t);
sinwave = sin(2*pi*f.*t);
plot(t,sinwave)

It seems to display a very different result when we plot the value obtained. After T the plot suddenly reverses, changes to original value of frequency and again increases in frequency I am not able to understand this behavior as the frequency seems to give values that increase and then decreases. Please help me understand why the plot for waves shows a different value and how to correct this. (Sorry not able to include plot due to web site restrictions)


texts:
It would be very helpful if someone could provide the intuitive reasoning and mathematical proof for:

$$F_a(t):=\frac{1}{a}\left[u\left(\frac{t}{a}+\frac{1}{2}\right)-u\left(\frac{t}{a}-\frac{1}{2}\right)\right]$$

$$  \lim_{a \to 0} F_a(t) = \delta(t) $$

where $u(t)$ is the unit step function and $\delta(t)$ is the unit impulse function.


texts:
I have an aerial image, like this

I want to detect the outline of the ship without shadow. How can I do it? I tried threshold in both graylevel, RGB and HSV. It doesn't work perfectly.


texts: Using a Kalman filter to clean acceleration data

texts: I am collecting acceleration, gyroscope and magnetometer data in all 3 directions. I would like to filter noise from the accelerometer using a Kalman filter. I am not trying to esimate orientatoin or position etc., I just really want to filter the acceleration signals. How should I design this filter?


texts: Particle filter, resampling step implementation

texts:
I am not very good in statistic neither in probability theory, and I do not completely understand the mathematical background of Monte Carlo simulation or SIR sampling, but I'd like computer vision and I'm experimenting with it hoping to get some nice results with already known techniques (I don't want to reinvent the wheel neither to find some new method).

I have done an approximate implementation of particle filter with Bhattacharyya distance as liklihood score function. What I mean with approximate refers to the resampling part.

Now, my rough algorithm does something like that, in python pseudo-code:

for each particles:
    p.position = move with a random walk motion model

    p_histogram = calculate histogram that is under the particle
    p.weight = bhattacharyya_distance( p_histogram, model_histogram )

sort(particles, less_weight_comparator)
# now particles[0] has less weight --> less distance
# and particles[nParticles-1] has more weight --> more distance

# find best k particles
# i.e., the particles that has weight less then a threshold 0.25
k = -1
for each particles:
    if p.weight < 0.25:
         k++;
    else:
         break

# resample all the other particles near those k best ones
if k!= -1:
    for i in range(k, nParticles):
        # choose a random number between 0 and k
        new_index = rand(0,k)
        # resample the particle near the k_th
        # with a gaussian random with 0 mean and 0.6 variance
        particles[i].position = particles[k].position * random_gauss(0, 0.6)
else:
    print "target is lost"

predicted_position = particles[0].position

this work quite well, but not every time, so I'm trying to improving it following some famous implementation, like the one in the Arulampalam_etal_2002 paper.

But.. I'm confused and I don't understand what is happening in the resampling stage.

With my rough and superficial algorithm I'm always aware of what is happening: particle with the weight less then my threshold are too far from the histogram model, and if k == -1 the target is lost. It is clear, in my mind.

In the paper mentioned (and in all the other particle filter papers) I'm reading about normalizing weight as all weights sums up to 1, and calculate a degeneracy value $N_eff$ (Equation 51) and.. I'm lost.  My n_eff has values from 249.444 (at the first frame after the particle filter initialization so that some particles has to be near to the target) to 245.452 (more frame later, when the target is lost) and here, again I cannot find a threshold. If I normalize the weights, after the sorting, I have particles weight that goes from values like 0.00418279 to 0.00240669, and I cannot easily find a threashold value that can help me in understanding which are near the target or not, and it depends on how many particles I have so if I add more particles I'm more lost.

More over, in paper like Nummiaro_2002 the Bhattacharyya is calculated inside a gaussian (Equation 10) and I'm not sure on what this means, and I lost an easy way to understand which particle has an histogram near to the model or not.

I'm totally lost in what including the Bhattacharyya into a gaussian means, and why normalizing and how to find the threshold for understand how many particles has to be resamples and which are out which are in.

I hope someone can explain me a little bit, with simple and easy words.


texts: $\tt awgn$ and $\tt wgn$ functions in MATLAB always output the same signal, no matter what I input as SNR

texts:
I'm very new to all of this.

Right now I'm just trying to add a bunch of different noises to an input signal and see how they affect my signal.

This is the current output of my code:

 signal,Fs = audioread('~/Documents/Methlab/sp03-0.wav');
 %generate white noise from input signal
 [m,n] = size(signal);
 snr = input('SNR: ');
 power = input('Power[DbW] ');
 whitenoise = wgn(m,n,power)

noiseSignal = addnoise(signal,noise,snr);
 %FILTER THE NOISE
%filteredSignal =
 VSLMS(filter_size,ss_uper,ss_lower,input_file,iterations)

subplot(2,2,1);
plot(signal)
title('Input signal')

subplot(2,2,2);
plot(noiseSignal)
title('White noise and signal sum')

If it helps anyone, I'm using this addnoise function to add the noise.

From what I understand if I change power and SNR the white noise should look different based on what I input right? What am I missing ?

I've tried a bunch of powers and SNR combinations. SNR being 0.0004 or 0.9 or power being 30 or 0.0003 had no effect on the signal output.


texts: Need help implementing DCT Type 4 algorithm

texts:
I'm trying to implement Algorithm 2 from this paper:

  Xuancheng Shao and Steven G. Johnson: Type-IV DCT, DST, and MDCT algorithms
  with reduced numbers of arithmetic operations, https://arxiv.org/pdf/0708.4399.pdf

It outlines an algorithm for computing an even-sized DCT type 4 via a DCT type 3 and a DST type 3 of half size.

I've made a few attempts at implementing this, but either the paper has a mistake or I have a serious misunderstanding about what's happening here.

https://pastebin.com/uFiHsMKX Here's my C++ code, which uses FFTW to compute the inner transforms, as well as to compute an "expected output". Aside from FFTW all the code necessary to compile is in this one file. The function compute_dct4 is what does the actual processing -- I did my best to name variables the same as in the paper, so it should be clear what in the code corresponds to what in the paper.

I'd appreciate some extra eyes on what's going wrong. Thanks!


texts: how to apply Short-time fourier transform to code in matlab

texts:
Please Help, I have a white noise signal that I created and need help to apply the Short-time fourier transform to my code so that it can do the bandpass filtering. I am trying to put it in the code so that I don't have to use the FDA tool. Also be able to plot the graphs from the output of the STFT being applied.

L = 5000; %Sample length for the random signal
Pause = 10000; %Sample Pause Gap
mu = 0;
sigma = 2;

%Need to see left signal is not displaying
Left_signal = sigma*randn(L,1) + mu;
Right_signal = sigma*randn(L,1) + mu;

Long_signal = [Left_signal zeros(L,1); zeros(Pause,2); zeros(L,1) Right_signal];

%Plots subplots in graph
figure
subplot(211);
plot(Left_signal, 'b'); grid on;
subplot(212);
plot(Right_signal, 'r'); grid on;


texts: I am working in w mbus and came across this term kchip/s in datarate like 37.768 kchip/s what does the unit specify?


texts: I have 2 color images generated from same X-ray machines, but one of the image is obtained using VGA cable. Therefore the RGB intensities differ in the 2 images. I am using Object detection to detect threat objects in the x-ray images. Is there any way i can remove the noise added due to VGA so that the images have similar RGB intensities?


texts: Image feature extraction for discrimination

texts:
I am using linear/quadratic discrimination to guess who is the artist of a painting. All the images are RGB format and I use matlab. I have used both original image and images that are resize to be square.

I tried to use features like the average of color channels, average skewness an so on. Then I plot them against each other and so far it is not very easy to use as feature for a linear/qudratic discrimination. Here are examples of the mean and norm of saturation and value channel (hsv in matlab) plotted against eachother for two artist:

As you can see there definitely is a difference in distribution but the error rate is to high.

Do you have any suggestions on what kind of features I can look at that would differentiate the artists even more?

When I say average I mean that I use mean2 to get the average of a channel or tranformation.


texts: Pitch fluctuates unexpectedly in simple sine function

texts:
I wrote a very basic C program to play a sine of a user-specified frequency. In the interest of portability, I have it spit values directly to stdout, so hopefully you can reproduce my problem on your own machine. This is the code I wrote:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

#define RATE 48000
#define TWOPI (3.14159 * 2.0)

int main(int argc, char** argv) {
        float frequency = atof(argv[1]);
        float phase = 0.0;
        float increment = TWOPI/RATE;
        int16_t sample;
        while (1) {
                sample = sin(phase) * 32767;
                putchar(sample & 255);
                putchar(sample >> 8);
                phase += increment * frequency;
        }
        return 0;
}

It can be compiled with:

cc sine.c -lm

and played with

./a.out 440.0 | sox -r 48000 -c 1 -t s16 - -d

That is: 16 bit mono audio at a sampling rate of 48,000 Hz.

I've tested this on both OpenBSD and MacOS and have noticed the same behavior: at around 6 seconds, the pitch of the wave jumps slightly. This happens again around 12 seconds. Further changes will take place as the wave continues to play. I have also replaced the sin() while loop with a prebuilt wavetable of a sine accessed through linear interpolation, only to notice the same problems.

Any ideas what might be off here?


texts:
Consider the following analog signal

$$x(t) = 2 \sin(100\pi t)$$

The signal $x(t)$ is sampled  with a sampling rate $F_s=50\textrm{Hz}$. Determine the discrete time signal. Plot the discrete time signal.

Also determine the total number of samples.

I don't understand how to approach this question.


texts: How to implement a moving average in C without a buffer?

texts:
Is it possible to implement a moving average in C without the need for a window of samples?

I've found that I can optimize a bit, by choosing a window size that's a power of two to allow for bit-shifting instead of dividing, but not needing a buffer would be nice.  Is there a way to express a new moving average result only as a function of the old result and the new sample?

Define an example moving average, across a window of 4 samples to be:

ma <= (a + b + c + d) / 4

Add new sample e:

ma_new <= (a + b + c + d) / 4 - (a / 4) + (e / 4)
ma_new(ma, oldest_sample, new_sample) <= ma - (a / 4) + (e / 4)


texts:
Are there any command line programs for Windows, preferably free and stand-alone, which can report the peak/strongest frequency within a given range of frequencies?
I need something like this to automate finding the frequency of a calibration signal which slowly drifts.


texts: As far as I understand, you can have two different continuous-time signals with the same discrete-time frequency spectrum after they are sampled and it may be possible these shifted replicas in the discrete frequency domain do not overlap with each other. In textbooks aliasing is associated with overlapping of shifted replicas. is the case I described a case of aliasing?


texts:
I need to calculate the output power spectral density of the following digital filter

My calculations are as follows:
$y\left(n\right)\:=\:x\left(n-1\right)+d\left[x\left(n-1\right)+x\left(n\right)\right]$
$Y\left(z\right)\:=\:X\left(z\right)z^{-1}+d\left[X\left(z\right)z^{-1}+X\left(z\right)\right]$
$Y\left(z\right)\:=\:X\left(z\right)\left[\left(1+d\right)z^{-1}+d\right]$
$\frac{Y\left(z\right)}{X\left(z\right)}=\:H\left(z\right)=\left(1+d\right)z^{-1}+d$
Did I get this equation right? would somebody help me please because this diagram is kind of confusing for me.


texts: Obtaining a horizontal disparity image from a depth image

texts:
In the paper Learning Rich Features from RGB-D Images for Object Detection and Segmentation available here, the authors encode depth images using three features computed from the depth image, called HHA images.
The first feature they compute is the horizontal disparity feature. The dataset contains the intrinsic and extrinsic parameters of a depth camera, in this case I believe since the paper tests their method on the NYU-depth v2 data set, obtained from a Kinect v2.
My question is what is a horizontal disparity image and how is it computed. I seem to find results for two images (stereo vision), but not for a single depth camera image. Also why is it called "disparity"?
Edit: Examining the matlab code provided by the authors, it seems what they do is to simply take the reciprocal of the depth, and then scale back to numbers in the range [0, 255]. This seems to ignore the calibration parameters of the camera, namely the focal length and the distance between camera lens centers. In this post here the mathematics of showing that disparity is inversely proportional to depth here.
I'm not sure what role the scaling factor of focal_length * distance between cameras plays.
I have included some python code here to describe my understanding of how to compute this transformation for a single gray-scale value.
max_value_unscaled_disparity = 1.0
min_value_unscaled_disparity = 1/256

unscaled_disparity_range = max_value_unscaled_disparity-min_value_unscaled_disparity

min_value_depth = 1.0
max_value_depth = 256.0

min_value_scaled_disparity = 0.0
max_value_scaled_disparity = 255.0
scaled_disparity_range = max_value_scaled_disparity-min_value_scaled_disparity

def getScaledDisparity(depth_value):
    depth_value = depth_value + 1.0
    unscaled_disparity_value = 1/depth_value
    scaled_disparity_value = int(((unscaled_disparity_value-min_value_unscaled_disparity)*(scaled_disparity_range)
                                  /(unscaled_disparity_range)) + min_value_scaled_disparity)
    return float(scaled_disparity_value)


texts: Auto or cross correlation for synchronisation in MATLAB

texts:
I have multiples records of  a motor sound, I would like to synchronize them first and then  calculate their spectrogram. My question is how can I use the autocorrelation or the crosscorrelation to synchronize them? I want to make sure that the only difference I see in the spectrogram is related to the frequency, not to the start position of the motor or the records.

I know to calculate the auto and cross correlation in MATLAB but I don't understand how to use them for synchronization purpose.

Thanks for any help!


texts: Z-transform of difference equations and stability of a process

texts:
According to this paper:

  $y(t)$ is stationary if all of the roots (of characteristic equation) lie outside the unit circle

Here, $y(t)$ is causal.

To me it seems the case is exactly the opposite, the roots should be inside the unit circle for causal process to be stable.

Example, consider the difference equation:

$$y_t-5y_{t-1}+6y_{t-2}=0$$

The characteristic equation is then: $$r^2-5r+6=0$$

Solution for characteristic roots: $$r=2, r=3$$

Solution for the difference equation (assuming 1 as initial conditions): $$y(t)=2^t+3^t$$

Clearly the roots are outside the unit circle (as 1<2,3), and the process unstationary (as it grows with t). To me it seems that the 3rd equation in the paper should have negative exponent for zs (in which case the solution would be the same as mine, even if arrived through the z transform).

Yet, the claim appears in many places so it can't be wrong (I doubt the paper is wrong either). However, for LTI system to be stable the poles must be inside the circle, a other claim that is repeated everywhere. Is the statement just plain wrong or are these things somehow not equivalent?


texts: How can the resolution of an image in matlab be reduced by using 2D haar wavelet by downsampling-decimation? Is there any predefined code or tool within wavemenu? thanks in advance


texts: How can I detect clicks/ticks in a wav file?

texts:
I have a long recording of a train (here is a small sample) there are a number of clicks as the wheels go over joins in the tracks. I would like to detect the location of these.

Looking at the waveform:

There appear to be 8-10 clicks. Some are close together.

I can see that they appear as regions of sustained high amplitude. Is there a filter I could cross-correlate this with that would pick them out or is there another better method?


texts:
Suppose I generate chirp signal with the following python code:

 t = np.linspace(0, 10, 512)
 f0 = 2
 f1 = 5
 fs = 2*f1
 w = chirp(t, f0=f0, f1=f1, t1=10, method='linear')
 result  = rfft(w)
 phase = np.zeros(257)
 mag = np.zeros(257)
 for binnum in range(len(result)):
   mag[binnum] = abs(result[binnum])
   phase[binnum] = np.arctan2(result[binnum].imag, result[binnum].real)
 plt.plot(fs * np.linspace(0, 0.5, 257),np.unwrap(phase))
 plt.show()
 plt.plot(fs * np.linspace(0, 0.5, 257),mag)
 plt.show()

And the phase plot is:

Magnitude plot is:

But it is very different from what is described here.

Well, looking at the phase graph I can see parabola from 0 to 0.8. The magnitude is definitly  not a constant.

What I'm doing wrong or misunderstood?

Thanks in advance.


texts: How we can show that 2 correlated videos are similar?

texts:
Can I tell that different versions of the same video are correlated to each other?
and if so, how can I show this correlation? for example,
high-quality and poor-quality versions of the same video.
Can we say various episodes of the same tv series are correlated since they mostly share the same background, characters, and scene?
and if so, how I can extract this similarity?
Maybe it's better to  ask, how can I compare two videos and show that they are statistically similar or they are correlated?


texts:
Consider a filter which tracks $x_i$, clipping $|x_i - x_{i-1}| >$ some limit, say 1:

in:  0 1 2  10 0 0  10 10 10  0 0 0 0
out: 0 1 2   3 2 1   2  3  4  3 2 1 0  -- |delta| <= 1

(in Python with Numpy:  cumsum( clip( diff( x ), - limit, limit ))

Can anyone point me to a description of such filters,
either in theory or in practice ?
Or is there no theory -- for the task of smoothing concrete $x$,
just vary the limit and see what happens ?

(Added: "This tutorial was written for normal engineers, who do not have nonlinear filters for breakfast." -- F.Daum, Nonlinear filters: beyond the Kalman filter. Is this online anywhere ?)


texts: What is the steady state response of a system of two exponentials to sinusoidal excitation?

texts:
What is the two exponential steady state response?

Here is an example solution for sinusoidal excitation of a system having a single exponential response to a impulse excitation:

Example:

Excitation: impulse: $L(t) = \delta(t)$

Impulse response of system: $I(t) = I_o e^{-t/\lambda}$

Excitation:  sinusoid: $L(t) = a + b \sin(\omega t)$

Sinusoid response of system:  $I(t) = I_o \lambda \left(a + \frac{b}{\sqrt{1 + (\omega \lambda)^2}} \cos(\omega t - θ)\right)$

$\tan(\theta) = \omega \lambda$

Question:  What is a similar solution for a system that has a double exponential impulse response?

Excitation: impulse: $L(t) = \delta(t)$

Impulse response of system: $I(t) = I_{o1} e^{-t/\lambda_1} + I_{o2} e^{-t/\lambda_2}$

Excitation:  sinusoid: $L(t) = a + b \sin(\omega t)$

Sinusoid response of system: ????

Thanks!!!!


texts: Problem understanding steps towards defining the convolution sum

texts:
I have a problem understanding the reasoning behind a step that was taken to characterize a LTI system.

So, we were told the following:

For each integer $k$, we have the following function: $$   \delta_k[n] =
\begin{cases}
1,  & \text{for $n=k$} \\
0, & \text{otherwise}
\end{cases}$$

Thus, if we write $x_k = x[k]$, we have that: $$x[n] = \sum_{k\ = \ -\infty}^{+\infty} x[k]\delta_k[n]$$ where x[n] is a discrete signal.

Finally, we can exploit the linearity of the system as follows: $$H\{x[n]\}\  = \ H\left\{ \sum_{k\ = \ -\infty}^{+\infty}x_k \ \delta_k[n] \right\}$$
$$=\sum_{k\ = \ -\infty}^{+\infty}x_k \ H\left\{\delta_k[n]   \right\}$$ where $H$ is the system.

It is this last step that confuses me a lot. How does linearity explain that we have $\sum_{k\ = \ -\infty}^{+\infty}x_k \ H\left\{\delta_k[n]   \right\}$ ? It only makes sense if we consider $x_k$ as a scalar (or a constant if you prefer), but we defined $x_k$ as a signal, so it can't be a constant. So, why is this correct?

PS: Just for clarification, we have that $\delta_k[n] = \delta[n-k]$ where $\delta[n]$ is the Kronecker-delta function.


texts: Error introduced by smoothing filter at the end of a processing block

texts:
I'm running analysis on a multichannel audio signal and due to the size have decided to process it in blocks (my computer doesn't have enough memory to process it in one go). Some of the data produced by the analysis needs to be smoothed.
The smoothing function means that each sample was an average of the samples contained within a window and the window being centered around sample n. At the edges this obviously isn't possible so the at the start of the signal the window remains static until such point that sample n is in the centre of the window, and then the window can slide along the signal - this also means at the end of the signal there will come a point where the window can't move along and will remain static. So essentially there will be a number of identical values at the start, and then at the end of the signal. This happens at the end of each block (as it can't continuously slide in the same way it would if processing the whole signal in one chunk). I know in order avoid discontinuties at the start of each block an overlap add/save method can be used, but I'm not sure how to fix this problem that occurs at the end of a block.
Just to explain further how the smoothing window moves along the signal, if we have a window length of 5 samples, the first 3 averaged values would be identical, at the 3rd sample, our current time step would now be located at the centre of our window, when moving to sample 4 the window would also slide along by one, keeping our current time step at the centre of the window, and now producing a different value. Once we get to the end of our signal the window will remain static for the last 3 samples, again producing identical values. This is fine when processing something in one big block, but produces discontinuties when trying to process in multiple blocks.


texts: How to determine the number of multipath?

texts: In a multipath propagation environment if we take the time index of reach channel tap by the integer multiples of sampling periods, assume the RMS delay spread $ \sigma_{tau} $ and maximum excess delay $ \tau _m $ is known, then how to determine the maximum number of paths?


texts:
I've been Googling "FIR filters" and "FFT convolution" etc. and I ran across something interesting.

Array Index vs. Pointers:
http://www.bores.com/courses/intro/program/7_array.htm
http://www.dspguru.com/dsp/tricks/fir-filtering-in-c

So, here's my question:
Why is almost every "Fast!!" Fourier Transform function on the web using Array Indexing?
(I just thought I would throw it out there.)


texts: sources for learning Image processing with MATLAB?

texts: I have Digital Image processing (DIP)  in my academics and I am looking for some reliable and daily updating sources for learning DIP Filters for Enhancements,color conversions, Filtering purposes through MATLAB , please suggest me some websites or tutorials which can help me. Thank you.


texts:
The images given below are a result of fruit defect detection algorithm i have made myself...

Just wanted to know is there a way of outlining these defects?

Original image 1:

1st Image:

Original image 2:

2nd Image:

What i want to achieve is this :

        =>      2.


texts:
Looking at some examples in the "handbook of MRI Pulse Sequences" by Bernstein, I see the following:

So if I can generate an RF pulse with a bandwidth of 2 kHz and I want a 3 mm slice thickness, the gradient amplitude can be chosen with:

$$
G_z = \frac{2\pi \Delta f}{\gamma \Delta_z}
$$

Plugging the values, this gives a gradient strength of 15.6 mT/m.

Now, my question is that will this gradient select a 3 mm thick slice at the magnet isocenter?

Now say in my MR experiment I am acquiring 128 such slices. Now, if I want to excite the "bottom most slice" than I need to adjust the RF bandwidth and the change in RF BW is given by:

$$
\delta f = \frac{\gamma G_z \delta_z}{2 \pi}
$$

Here now, my $\delta_z = -64 * 0.003 m$. Hence, the new RF BW will need to be

$$
f = \Delta f - 64 \Delta f = -126 khz
$$

I am not sure if this calculation is correct. The RF value seems quite high to be although I am not sure what is the range of modern MR scanner RF generators.


texts: What is the rms value of a signal and how to use it

texts:
In some papers/books I found infos about a signal given as RMS value.

What is the RMS value and what information does it contain?

How can I compute it?

Moreover, if they say that the signal is a white noise with RMS $100$ band limited onto the interval $[0, 2000] Hz$, how can I use the RMS information to produce the signal they are talking about?
And if they say that the noise contributes to 5% of the signal RMS value?


texts:
There is probably something trivial behind this, but I am missing something. I need to create a stationary random time series data v(t) which is the the sum of another time series u(t) and u(t) with a time shift $\tau$, i.e

$v(t) = u(t) + u(t - \tau)$

The series u(t) is itself a stationary colored Gaussian time-series with some PSD $S_u(f)$. No problem! I wrote a small function (using Python) to make u(t) for a given PSD. I calculate the PSD of the the generated time series using Welch's method and it checks out.

If we take the Fourier transform of u(t - \tau) it should just gain an extra phase compared to the Fourier transform of u(t)

$F[u(t - \tau)] = \int dt \, u(t -\tau) \, e^{2 \pi j f t } =  \int dt \, u(t) \, e^{2 \pi j f t }  e^{- 2 \pi j f \tau } = \tilde{u}(f) e^{- 2 \pi j f \tau }$

Since $v(t) = u(t) + u(t - \tau)$,

$\tilde{v}(f) = \tilde{u} (f) \, (1 + e^{- 2 \pi j f \tau })  $

And the PSD of V should then be (using Euler identity)

$S_v (f) = 4 S_u (f) \cos^2 (\pi f \tau )$

In particular $S_v$ should vanish when $f = \frac{2n + 1}{2 \tau}$ with the smallest zero being at $\frac{1}{2 \tau}$. However when I test my this with the u(t) series I generated above the PSD looks different. There are zeros low frequencies which shouldn't exist. For example  here is the PSD of the v(t) series and the spectra I would expect for $\tau = 2$. (The first zero should be at 0.25 Hz)

Moreover the zeros in the data seems to change with the sampling frequency I am using. The v(t) series for the above plot has a sampling frequency of 0.25 Hz. Here is one with 0.5 Hz, keeping everything else the same.

I am kind of baffled by this. Any help would be greatly appreciated. Sorry about the long post, and thanks in advance!!

EDIT: After using a higher sample rate as suggested by Hilmar below, and fixing the error I described in the comment everything looks good. Spectrum for the fixed case below.


texts:
In LTE systems, up- and downlink are duplexed either using FDD or TDD. Why has this not been designed to utilize the same multiplexing technique as used for the downlink (OFDMA)? Resource blocks could be dynamically assigned for up- and downlink then, making everything much mor flexible.
I know that for uplink sc-fdma is utilized, to minimize PAPR, is this maybe the only reason?


texts: extract trend correctly, including most recent values

texts:
I'm looking to extract the trend of a signal.

i've tried two methods for now, polynomial regression, and wavelet denoising

both methods don't respect the computation of the last values (meaning the last values computed will not be the same if we compute a longer buffer containing new values).

is there a way, excluding FIR/IIR filters to extract a trend that stay consistent across all the values in time ?

matlab code is welcome if possible

here, with Matlab denoising, we can see most recent values aren't the same

thanks


texts: inverse discrete FFT in python, multiple times?

texts:
I was wondering what really happens when taking the inverse discrete FFT on some set of numbers, for 3 times? Because looking at it, it looks like we're getting an output that is identically with the input, and I don't get why.

Also I was trying to implement my own inverse discrete Fourier transform in python, following wikipedia, looking at the IFFT provided by octave:

and at the plot provided by my code:

It just doesn't make any sense.
Here is my code:

#!/usr/bin/python
import matplotlib.pyplot as plt
import numpy as np
import math

points = [24+0j, -18.39230+0.73205j, 8-3.46410j,
        -2+4j, 0+0j, 2.39230-2.73205j, -4+0j, 2.39230+2.73205j, 0-0j,
        -2-4j, 8+3.46410j, -18.39230-0.73205j]

def inverse_fourier_transform(points):
    inversed_points = []
    N = len(points)

    for k in range(0, N):
        s = 0.0
        for n in range(0, N - 1):
            s += points[n]*np.exp(1.j*2.0*np.pi*k*n/N)/N
        inversed_points.append(s)

    return inversed_points

inv_points = inverse_fourier_transform(points)
inv_points2 = inverse_fourier_transform(inv_points)
inv_points3 = inverse_fourier_transform(inv_points2)

x = []
y = []

for pt in inv_points3:
    x.append(pt.real)
    y.append(pt.imag)

plt.scatter(x, y, s=0.5)
plt.plot(x, y)
plt.show()


texts: What is a segmentation mask in the paper "Fast Edge Detection Using Structured Forests"

texts:
I've read the paper about edge detection, in this paper they treat edge detection as a learning problem which takes an image patch as input and output a label, a binary edge map or a segmentation mask, but I don't know what is a segmentation mask?

In the paper their image patch is 32*32 but segmentation mask is 16*16, why is that?

During training they transform a multidimensional segmentation mask into a binary label to calculate the information gain, but when making predictions, how can I transform a binary label back to a segmentation mask or edge map? the transformation seems to be irreversible.


texts: Do you recognize this discrete realization of a low pass filter?

texts:
I have a simple 1st order lag filter implementation that performs well, shown in  the block diagram below.  This is legacy 'code' we've been using for many years, the original developers have long since moved on.  I've been asked to doc the requirements (development of the difference equation)  starting from the s domain and assuming $\frac{Y(s)}{U(s)} = \frac{1}{\tau s+1}$.
$k1$ is $0.5+\frac{\tau}{T_s}$, $k2$ is $0.5$. $\tau$ is the time constant, in sec, and $T_s$ is the time between samples, in sec.
My assumption was I would arrive at a bilinear transform of $\frac{1}{s+1}$ but it's just not happening as the topology is wrong.  Any hints?
Thanks,
Chris


texts: Filtering and Fourier Transforming, does the order matter?

texts:
I have a signal $x(t)$. I want to find the Fourier Transform of it, $X(f)$, and then extract a narrow frequency range from $X(f)$ by use of a Band Pass Filter (BPF) in frequency domain.

Can I instead filter $x(t)$ by using a BPF in time domain and then find the Fourier Transform of the filtered signal?

I believe these two are equivalent.


texts: What's the right scale to represent a frequency spectrum?

texts:
To represent a 16 bits temporal spectrum, it's "easy" : my scale is going from -32000 to +32000. But What's the scale for a frequency spectrum ?

I made that try : I calculated a pure single frequency sine signal and calculated it's energy for a period. So, I supposed that any signal generated on such a period - without saturation, of course - could have more energy. So, I considered the ratio

  20.log10 (signal_energy/maximum_theorical_energy)

as my 0 dbFS.

Then, I made few experimentations, generating signals for different kind of frequencies and, when I had a single sine, each time, I reached the 0dbFS.

But, of course, a music have thousands of frequencies. So, when I represent a spectrum of a "music sound" signal, I rarely reach more than -30 db and my signal is more or less around -60 db for each frequency. It seems very low.

So, my question is : is it correct to represent a signal with such a measurement ? If not, how do represent the "values" for each frequencies on the graph ? With which scale ?

Thank you for all your answers.


texts: What Resources Are Recommended for an Introduction to Signal Processing (DSP)?

texts: What are some recommended resources (books, tutorials, lectures, etc.) on digital signal processing, and how to begin working with it on a technical level?


texts:
I have a stereo music signal corrupted by strong sinusoidal noise that varies over time. Here is the spectrogram of Left channel I plotted with Matlab.
As you can see there are 3 or 4 strong harmonics with frequency that varies over time.

As a first try I computed the difference of left and right channels and it seems most of the harmonic noise disappears, so my guess is the noise in left and right are almost identical. Here is the difference spectrogram to prove this.

Which technique could I use now to filter away all the interfering harmonics while preserving the original music spectrum as well as possible?

Notch filters?
Tracking the frequencies over time and adding a sinusoid with reversed phase?

How do I estimate and use the phase information? I think phase is important to have good denoising, but would not know how to proceed.

Thanks a lot! :)


texts: Is the filter described by specific difference equation time independent

texts:
Is the filter described by
$$y[n] + 2 y[n-1] + y[n-2] = x[n] - 2 x[n-1]$$
time invariant?
($x[n]$ - input, $y[n]$ - output)

Theorem of Foures-Segal says no (?)

When I asked about Foures-Segal I meant the following.
F-S theorem says that if transfer function $T$:

$$y=Tu$$

is time invariant then there exists a function $G(z)$
which is analytic and bounded in exterior of unit disk
on complex plane,
such that

$$T = Z^{-1} G(z) Z$$,

where $Z$ is the Z-transform.

In my case the function $G(z) = \frac{z(z-2)}{(z+1)^2}$
having pole at $z=-1$, is not bounded. I would say
it means the system is not time invariant?


texts: FIR filter Gain of 0dB at Passband Frequency

texts:
I have designed a FIR pass band filter. The impulse response looks good and on applying this on my input data, the filter works as desired. Now I have a requirement in my system to have exact $0$ dB gain at passband frequency and hence I normalize the FIR coefficients using the following method:

for (n=0; n < N; n++)
  Gain += FIRCoefficients[n] * cos(2 * PI_DOUBLE * PassbandFreq * n);

for (n=0; n < N; n++)
   FIRCoefficients[n] = FIRCoefficients[n] / Gain;

Even after this normalizing, the Gain at passband frequency is around $23$ dB. What am I doing wrong?


texts: Reproducing paper results about a wavelet transformation using python

texts:
I have been reading this PhD theis about wavelets and I am trying to reproduce some of the results but I don't know the specific code to use to generate similar results.

The original time series looks like this:

Then, the paper says:

  The customer demand is transformed to the wavelet domain by the
  non-decimated wavelet transform performed with the Symmlet 8 filter.
  Figure 47 shows customer demand in the wavelet domain.

Can someone point me out how to reproduce this results using Python (then R, then matlab, in that order of preference)

NOTE: link to document, page 122: A Multiscale Forecasting Methodology for Power Plant Fleet Management A Thesis, Hongmei Chen, Feb. 2005.


texts: Integer Wavelet with Only Positive coefficients

texts:
Is there any Integer Wavelet Transform which produce only positive coefficients?
Currently i'm using Haar wavelet but it produces both positive and negative coefficients.


texts: how can I generate SIMO underwater acoustic channel

texts:
I'm working in underwater acoustic channel, then I need to generate SIMO channel to use it for simulation in MATLAB,  how can I generate SIMO channel?

thank you


texts:
Let $x[n]$ be real signal and $y[n]=\exp(j3\pi n)$ be a complex signal

Would the convolution between those two signals be $$x[n] * \Re(y[n]) + jx[n]*\Im(y[n])$$?


texts: Best DSP algoritms for ultrasonic background noise cancellation

texts:
I am looking into detection of pressurized gas leaks by means of the broadband hissing sound (white noise) generated from pressurized leaking gas. The instrument will have a broadband microphone followed by a high pass filter to filter away all low frequency audible noise below 15-20 kHz, but when broadband noise from leaks occure, they will pass the filter and activate an alarm. The problem is that in noisy industrial environments, machines (compressors), and other equipment will also have some ultrasound components above 15-20 kHz, and it will do false activation. Is it possible to recomend optimal DSP algoritms to be used, after the raw 15-20 kHz filtering to do a efficient suppression of machine made background noise, without affecting the sensitivity for the leak noise? I cannot tolerate varible gain for the leak noise (white noise) because then it will affect the sensitivity to detect leaks in a predefined distance, so the system must automatically suppress unwanted industrial background noise, but not affect the sensitivity to the leak noise (white noise in frequency range 20 kHz to 50 kHz).?

I would say that the main difference between leak noise i want to detect (the broadband white noise) and the machine noise is that the leak noise is just plain white noise, while very often the machine noise will have some components of high frequency reprocicating noise pattern in it from highspeed running compressors which the white noise leak noise do not have, but the high frequency reprocicating background noise can vary in frequenct from machine to machine so I do noy want to use some kind of fixed notch filter, so I am locked to one machine, at one rotation speed. Actually the problem is that the machines main noise frequency will be belos the 15-20 kHz, but it is the harmonics that repeat itself up above 15-20 kHz, and it is those I need to have filtered away. Seen from a more practical point of view, my own ears will be able to hear the difference, so I woundered if that could be realized in a electronic filter too.


texts: What is so special about this filter coef = [1,2,1]?

texts:
Can anyone tell what kind of filter is this? I can see it's smoothing/averaging the input signal, but why specifically with [1,2,1] instead of something else? What is so special about this coef=[1,2,1]? Thanks

# Python code
length = 1000
signal = np.random.rand(length)
output = np.zeros(length-2)
coef= np.array([1,2,1])
for i in range(length - 2):
   output[i] = np.sum(signal[i:i+3] * coef / 4)


texts:
I am trying to replicate the spetrum as given in Audacity and other commercial software... I am getting crazy with the correct scaling needed in order to have dBFS (with reference to a full scale sine wave).

Suppose $x_t$ is the PCM sample at $k$ bit precision, first of all I rescale the PCM samples

$$
s_t =  \frac{x_t}{  \frac{ 2^{k-1} }{ \sqrt{2}  }  }
$$

Suppose now $w_t$ is some window function, I compute the FFT of of size N of the product $s_t w_t$.

Let $P_1, P_2, \ldots, P_N$ be the FFT complex values. Since I am only interested in one-side spectrum, I just retain the fist half of the sequence. Moreover suppose for simplicity that the time series matches the window size $N$.

Taking

$$
p_i = 2\left| \frac{P[i]}{N}\right|^2
$$

for $i=1,2,\ldots, N/2-1$ should give the power spectrum values. Then one would take $10\log_{10}(p_i)$ to get dBFS.

But I think that in order to express dBFS one needs to scale the $p_i$ in some way taking into account the effect of the window. In fact  I produced a pure 1KHz sine wave then I computed my spectrum and compared it with that obtained with Audacity,  Adobe Audition etc. My FFT stays below by few dBFS so I think it's really a problem of scaling.

I hope somebody knows if there is a standard ways of doing that!

Regards
Pierre


texts: Filter Design near Nyquist frequency (complex curve fitting)

texts: I'm fitting IIR filters to a complex transfer function with Matlab's invfreqz.m. My fitting results are very good, but when I come very close to the Nyquist frequency the phase response runs to 0 and I also get deviations in the magnitude response. I did some research but couldn't find anything regarding that phenomenon. Can anybody point me in the right direction, for which keywords I have to search in order to find something on that topic and how I can try to optimize the fitting near the Nyquist frequency?


texts: I am not sure if this question suits this site, but I want to convert humming into notes or relative values to analyse melodies, word of "Signal Processing" looks suit for this. I am planning to use Flash AS3, and Tortoise and Hare for Cycle Detection. I need very simple script implementation. All I want is find notes more easily with my bass. Is that possible? Or is there any frameworks that already makes this?


texts: Why are the real- & imaginary part of complex white gaussian noise independent processes?

texts:
Assume we transmit a bandpass signal over an AWGN channel adding the Gaussian noise contribution $n$

\begin{equation}
Acos(\omega_ct) + n(t)
\end{equation}

Further, the bandpass signal is generated by IQ-modulation. If we perform homodyne IQ demodulation + lowpass filtering at the receiver should the resulting in-phase & quadrature component not be affected by the same lowpass noise contribution?
Instead the equivalent complex lowpass noise is modeled as

\begin{equation}
N_{LP}=X+jY
\end{equation}
where the imaginary & real part are both Gaussian distributed and independent processes. Hence the real & imaginary part of the transmitted symbol are affected by different noise contributions.

Why is that?


texts: How can we separate two frequency component when they get mixed means when signal get mixed with noise how receiver can judge which part is signal and which is noise ?


texts: Two versions of Constant Q Transform (CQT) doesn't match each other?

texts:
To my knowledge, there's two major CQT papers, the one by Brown in 1991, and the one by Schorkhuber in 2010.

The 2010 paper claims to be a more computationally efficient implementation of the 1992 and 1991 paper by Brown. (The 1992 paper is also an efficient version of the 1991 paper.)

I discovered that the mathematics behind the two papers are not exactly the same.

Definition of $Q$ and $N_k$ by Brown in 1991 and 1992

$Q=\frac{f_k}{\Delta f_k}=\frac{f_k}{f_{k+1}-f_k}=\frac{1}{2^{\frac{1}{b}}-1}$ ----(1A),
where $b$ is the number of bins per octave.

$N_k=\frac{f_s}{\Delta f_k}=\frac{f_s}{f_k}Q$ ----(2A), where $f_s$ is the sampling rate.

Definition of $Q$ and $N_k$ by Schorkhuber in 2010

$Q = \frac{f_k}{\Delta f_k} = \frac{N_k f_k}{\Delta \omega f_s}$ ----(1B)

$Q = \frac{q}{\Delta \omega (2^{\frac{1}{b}}-1)}$ ----(2B)

By combining equation (1B) and (2B), you can get the window length $N_k$ as below

$N_k = \frac{qf_s}{f_k(2^{\frac{1}{b}}-1)}$ ----(3B)

where $\Delta \omega$ is the equivalent noise bandwidth, $q$ is the scaling factor.

Here's the problem that I have

If you try to equate (2A) and (3B) together, you will get another expression for Q.

$\frac{f_s}{f_k}Q = \frac{qf_s}{f_k(2^{\frac{1}{b}}-1)}$

After simplifying, you will get the follow expression.

$Q = \frac{q}{2^{\frac{1}{b}}-1}$. Which looks very similar to equation (2B) with a missing $\Delta \omega$ factor.

So, there's 3 different expressions for $Q$ now? To me, it seems a contradiction. Either of the papers must be wrong in their definition on $Q$.

Another mystery is where does the term $\Delta \omega$ come from? I don't see any necessity in adding this term. It seems to me the source of the contradiction.


texts: How to make a signal conjugate symmetric?

texts:
Take the simple frequency-domain band-pass filtering operation below . . .

NFFT = 128;
x = randn(NFFT,1);
H = zeros(NFFT,1);
H(10:20) = 1;
y = ifft(H.*fft(x), 'symmetric');

This gives a real output because I'm using the conjugate symmetric flag to the ifft operation.

I want a function that returns the conjugate symmetric version of H, so I do not have to rely on the builtin symmetric option in Matlab's ifft. NFFT can be any positive integer. This could be called something like this . .

H(10:20) = 1;
H = MakeConjSym(H);


texts: Adding Poisson noise to a time series signal

texts:
I'm working on a radiation detection problem. I have simulated some measurement count data and want to add an amount of Poisson background noise to achieve a certain SNR. Radiation count measurement (due to the background radiation) is given by the Poisson random variable.

After searching the web, I'm not seeing any clear indication of how to go about it doing this. Is there a similar method as in the case of AWGN? (https://www.mathworks.com/help/comm/ref/awgn.html).


texts: Differences between two closed-loop systems

texts:
What could I say about stability of those systems by looking at the scheme?

System 1

System 2

$y_r(s)$ is input, $a(s)$, $b(s)$, $p(s)$, $q(s)$ and $r(s)$ are some polynomials and $y$ is output.

I calculated tranfer $\frac{y}{y_r}$ and it's the same for both systems.

$$\frac{y}{y_r} = \frac{r(s)b(s)}{a(s)p(s)+b(s)q(s)}$$

What makes the difference between stability and behavior of those systems?

Some progress:

I've edited the schemes. There is something, that obviosly differs in those systems. It's $e_{ss}$ - stady state error.

For the first system is:
$$e_{ss1}=y_r(s)\frac{r(s)}{p(s)}-y(s)\frac{q(s)}{p(s)}$$

For the second:
$$e_{ss2}=y_r(s)r(s)-y(s)q(s)$$

Does it tell me something about stability?

I guess this form is more usefull:

$$e_{ss1}=y_r(s)\frac{r(s)a(s)}{a(s)p(s)+b(s)q(s)}$$
$$e_{ss2}=y_r(s)\frac{r(s)a(s)p(s)}{a(s)p(s)+b(s)q(s)}$$

The second system's steady state error has higher order numinator than denominator. Does it mean, the signal's higher frequencies would be infinitely amplified?


texts:
I'm reading this material but I don't understand this example:

"An example of a 94% overlap transform processing could be demonstrated by a 279 sample frame with 262 samples overlapped from each frame into the next. At 12.5 MHz sample rate, each spectrum frame would occupy 23.76 µs. Each spectrum would start 1.36µs after the previous one and could show a 10 MHz span."
Could anybody please help to explain what is 10 MHz span and why why each spectrum could show a 10 MHz span?
The same problem with the following picture (I don't understand the example):


texts: How can an audio segment report an integrated loudness (LUFS) greater than 0

texts:
According to my very rudimentary understanding of the ITU-R BS.1770 standard, I thought that the LUFS (LKFS) value for any audio signals should not ever go positive.
However, there is a 5-second segment from an audio file I have whose integrated loudness level is +1.28 when measured using ffmpeg, and +1.44 when measured using pyloudnorm:
(stereo, source audio resampled at 48000 Hz, 16-bit signed little-endian integer waveform PCM)
[Parsed_loudnorm_0 @ 0x7fc245008680]
{
    "input_i" : "1.28",
    "input_tp" : "0.86",
    "input_lra" : "2.90",
    "input_thresh" : "-10.16",
    ...

How is this possible? Could this be an issue related to the implementations (i.e. ffmpeg and pyloudnorm) and not the algorithm itself?
I concede that I don't fully understand the math in the algorithm, and so if the algorithm does allow a positive value in certain situations, please explain.


texts:
What do ripples in the spectrum of a bandpass filter physically mean? What do they tell us about the filter? Please check the picture:

The filter should look smooth, like a resonator, such as


texts: Why does signaling overhead for time synchronization scale up with the number of transmitter nodes in a multiple access system?

texts:
I was recently reading about random multiple access methods without time synchronisation (of the pure Aloha type). In this paper it is stated that,

  The need for transmitter synchronization is a major drawback for large networks as the signaling overhead scales up with the number of transmitters independently from their traffic activity factor.

I understand that oscillators tend to drift in time and present instabilities depending on the technology, operating conditions, etc., and therefore a reference signal must be given by the central node with a periodicity that depends on the desired accuracy and the transmitter's clock drift.

I would like to know why there is also a dependency on the network size.

EDIT: The only scenario I can think of where signalling scales up with the network size, is a distributed network (but in that context, we usually talk about logical clocks and not about hardware clocks/oscillators). Note also that the network described in the paper has a star topology, so a periodic broadcast message from the central node would do, independently on the number of users in the network.


texts: Use of Tukey window and its effect on a graph

texts:
I am trying to minimise the noise in a spectrum of a variable. I am trying to to divide the time history in equal parts, to do the FFT of this time history parts and to do an average between these fft.
Doing this procedure I have find that the trend, at the high frequencies seems to be a little bit different from the one without doing this procedure. Following an other suggestion of my advisor I have multiplied every part of the time history by a Tukey window, using the command with matlab
A= tukeywin(L)
Using the Tukey window, I have find that the trends seems more similar with the first one. Can someone explain me why? I have tried to search the theory of this window, but I always find the definition of the filter, but I cannot understand which could be the effect of this window.
A other fact that I have noticed is the fact that obviously using the tukey window the energy content is reduced, are there technique to reduce the noise avoiding this reduction?
I just need a suggestion to start, it is difficult to understand which is the correct technique to reduce the noise.


texts: When to use cross spectral density versus cross correlation?

texts:
All my life, I've used cross-correlation to analyze the strength of the relationship between an input signal and the output signal.  I use this to determine which input signal most closely influences the output signal:

Yesterday, a friend told me that he uses cross spectral density (same as coherence, he say).

To my interpretation, they seem to quantify the same thing.

When would someone use one over the other?


texts: Convert normalized frequency to real frequency in AR model

texts:
Let us suppose that we have modeled signal using AR model, and suppose we have following model:

I used spectral estimation function from MATLAB

pyulear

Now frequencies are given in normalized frequencies and I would like to know how to convert them back to real frequencies? From there it is clear that we have four deterministic model and also plus some white noise, actually I want to know approximate frequencies in each deterministic model. I can of course determine this frequencies using  FFT, periodogram and so on, but I am studying application of the AR/ARMA model. So in my case, I have such frequencies on following figures, how can I determine actual frequencies? Thanks in advance.


texts: Noise Shape Filter to obtain a given PSD

texts:
I have experimental PSD data of disturbances in terms of acceleration on a mass. My first problem is that these data are presented as $ (m/s^2)/\sqrt(Hz)$ which is not the PSD unit. Those unit should be adressed to an amplitude spectral density (ASD) which is directly $ ASD = \sqrt{PSD}$. (I am not sure of that either, beacuase in some books the relations takes into account $\sqrt{2}$.
However, my idea is to reconvert into a transfer function (TF) the PSD and then use that TF as a filter. Infact, if the TF is fed by a white noise ($ W(f) = 1$) the output of the filter operation is $ H(f)*1$ i.e. the colored noise in frequency domain $C(f)$.
Thus
$$ PSD =|H(f)|^2 $$
So I can interpolate the given points and use a least square method to obtain the chosen number of coefficients of the TF
$$ min(a_z, b_p) \sum_{i} (|H(2\pi j f)|- \sqrt{PSD}|)^2 $$
( I don't really know how to do this).
My other option is to hit and trial whit some tf
Below, there is an image of the data.

To finish I am leaving here what I have done on Matlab

% TRIAL AND ERROR PROCEDURE
% clear all
% close all

% EL ACT data
% points
x = ([0.06 0.1 0.12 0.2 0.24 0.4 1 2 3 10].*1e-3);
y = ([500 150 100 40 30 20 13 12 12 12].*1e-15).^2;     % now is a PSD?
% plot(x,y)           %trial and error to make the interpol suitable
%
loglog(x,y)
hold on
loglog(x,y,'o')

%%
xlog = log10(x);
ylog = log10(y);

% fitting
N = 1000;
pp = polyfit(xlog,ylog,3);
freq_log = linspace(xlog(1),xlog(end),N);
ELEM_sensing_ASD_log = polyval(pp,freq_log);
ELEM_sensing_ASD = 10.^ELEM_sensing_ASD_log;
freq = 10.^freq_log;
loglog(freq,ELEM_sensing_ASD)  %fitted?

%%%%%%%%%%%%
% Data from Ref
s = tf('s');
H_EL = ((s+2e-5)*(s+3e-5)*(s+3e-4))/((s+2e-6)*(s+5e-6)*(s+1e-5));
Num = cell2mat(H_EL.Numerator);
Den = cell2mat(H_EL.Denominator);
[A,B,C,D] = tf2ss(Num,Den);

%% Simulink
T_filter = 100;
% open('noise_shape_trial_SIM')
sim('noise_shape_trial_SIM')

%%
segment_lenth_sample = 1e4;
overlap = 0;
DFT_points = 4e4;
f_sampling = 1000;
[PSD_color,f] = pwelch(color_time); %,segment_lenth_sample,overlap,DFT_points,f_sampling);
loglog(f,(PSD_color))
grid on


texts: Fourier domain division causing translation

texts:
I am trying to find the correct filter to convolve with an image so they have the same PSF.

I have the final PSFs of both images:

A: PSF of the image with the wider PSF.
B: PSf of the image with the narrower PSF.

(So the FWHM of B is smaller than A).

For the time being I am trying to do the most simple (naive) deconvolution:  by simply dividing the Fourier transform of A with that of B. The reason is that the type of data that I am working with (astronomical images) will not have zeros in the fourier transform. Also I am using very high signal to noise images of stars (point sources) so the effect of noise is minimal. Later on, I will try to implement the Wiener deconvolution. But my problem seems to be more basic than that!

When I take the naive deconvolution steps, it seems the Final PSF I get is translated by half the image width as you can see below. Just a note: I have set the sums of both images to 1 before this process. Also, the images are all inverted so zero is shown with white.

This is A (the wider PSF):

This is B (the narrower PSF, same image size as A):

This is the Fourier transform of A:

This is the Fourier transform of B:

This is the pixel by pixel division of the Fourier transform of A divided by the Fourier transform of B. Following the advice in this stackexchange-Mathematics question, I have set all the pixels with a very small B value to zero (a plus-like figure in the center that is all white).

And this is the inverse Fourier transform of the division:

It seems reasonable, but it is not centered, its as if it has been translated by half the image size! I know I can multiply the Fourier domain results with constants to center it. But I don't understand why this happens. I am affraid I might have made a mistake somewhere! If I understand it, then applying that translation is easy.

I am using this exact algorithm for convolution (where the Fourier transforms are multiplied) and the inverse discrete Fourier transform result is centered. So, I know there is no error in my inverse discrete Fourier transform steps.

As a summary, my question is that why does the inverse Fourier transform of the frequency space appear to give a translated (by half the image size) spatial image in division?


texts: Is ICA appropriate for separating mixed signals when all source signals are NOT detectable at all sensors?

texts:
A generic implementation of ICA for the separation of a mixture of $N$ signals into their $M$ constituent components requires that the signals be assumed to be a linear instantaneous mixture of the sources.  Every description of ICA that I have come across seems to take for granted the fact that all $M$ sources are present to some extent in all $N$ signal mixtures.

My question is, what if the $M$ sources are only present in some but not all of the signal mixtures?

Does this scenario violate the fundamental assumptions necessary for ICA to be able to separate these signals?  (Assume, for the sake of argument, that we are dealing with an overcomplete or complete system ($N>M$ or $N=M$), and that each of the $M$ source signals are in fact statistically independent from each other).

The implementation I am considering using ICA for, in which this situation arises, is the following: I have data from 4 different types of sensors, each with a different number of channels.  Specifically, I have 24 channels of EEG data, 3 channels of electrooculographic (EOG) data, 4 channels of EMG data, and 1 channel of ECG data.  All data is recorded simultaneously.

I would like to identify the contributions of the ECG, EMG, and EOG signals within the EEG data so that I can remove them.  The expectation is that EMG+ECG+EOG signals will be picked up by the EEG sensors, but not vice versa.  Also, EOG and EMG will likely contaminate each other and be contaminated by ECG, but ECG will probably be pretty isolated from all the other signals.  Also, I'm assuming that where mixing does occur, it is linear and instantaneous.

My intuition tells me that, hypothetically, ICA should be smart enough to return mixing filters with very small (close to 0) coefficients to account for the lack of contribution of sources to a mixed signal.  But I am worried that something about the way that ICA demixes the signals inherently enforces the expectation that all sources will be present in all mixtures.  The implementation I am using is FastICA, which is a projection pursuit-based approach.


texts:
I am working on speech synthesis with LPC and I originally implemented a pre-emphasis filter:

+(void)processBuffer:(Buffer *)buffer {
    float alpha = -0.9375f;
    for (int i = 1; i < buffer.size; i++) {
        buffer.samples[i] += buffer.samples[i - 1] * alpha;
    }
}

And found that my output audio was distorting, so I adding some code to scale the values, and so now it looks like this:

+(void)processBuffer:(Buffer *)buffer {
    float preEnergy = [self energyFor:buffer];

    float alpha = -0.9375f;
    for (int i = 1; i < buffer.size; i++) {
        buffer.samples[i] += buffer.samples[i - 1] * alpha;
    }

    [self scaleBuffer:buffer preEnergy:preEnergy postEnergy:[self energyFor:buffer]];
}

+(void)scaleBuffer:(Buffer *)buffer preEnergy:(double)preEnergy postEnergy:(double)postEnergy {
    float scale = sqrt(preEnergy / postEnergy);

    for (int i = 0; i < buffer.size; i++) {
        buffer.samples[i] *= scale;
    }
}

+(float)energyFor:(Buffer *)buffer {
    double sum = 0.0;

    for (int i = 0; i < buffer.size; i++) {
        sum += buffer.samples[i] * buffer.samples[i];
    }
    return sum;
}

I am finding that the speech I am generating still doesn't sound very good, and so I am just wondering if my implementation for this filter is ok?  Does anyone see anything I am doing wrong or does this look ok?


texts: attenuation of the central difference differentiator

texts:
In the following snippet, I am differentiating a sine wave using the central difference equation,
first through misc.derivative and then through convolution with the kernel [1/2, 0, -1/2]
import numpy as np
import scipy.misc as misc
import scipy.signal as sig
ncycles = 5.0
nsamples = 500.0
x = np.linspace(0,ncycles*2*np.pi,nsamples)
dydxmisc = misc.derivative(np.sin,x,dx=0.1)
dydxconv = sig.convolve(np.sin(x), [0.5,0.0,-0.5], mode='same')

As shown in the plots below, misc.derivative produces the expected output when dx is small enough, but when using convolution the derivative is attenuated.

The attenuation changes based on number of cycles and number of samples and is a factor of pi but I cannot figure out what the exact relationship is. Can anyone explain why the attenuation happens and what the relationship is? In other words, how do I have to scale the kernel in order to get a derivative of the correct dimension?


texts: I have a embedded system where I get some type of sensor readings with varying frequency between 0.3 to 1 Hz. Sometimes these readings may have impulse like noise. To filter out this noise, I've implemented a median filter with total 3 data points. Since noisy values does not occur repeatedly, my filter works well except it adds single reading time as delay. The microcontroller that I'm using has limited sources with no Floating-point Unit. I am looking for alternative filters which are closer to real-time and does not consume much resource.


texts: Object-Recognition using SIFT: Switching the order of images creates two different results

texts:
I've developed a large test-system which can use every interest-point-detector and -descriptor. I use it for evaluation of the given algorithms with the standford-image-collection of cd-covers. This post will not go any deeper with any source-code. I only found out that switching to images (no referred as firstImg and secondImg) will create to different results. The Ransac-Threshold is used to compute the homography and epsilon is the reprojection-error. For this example I use the Scale-Invariant-Feature-Transform by David Lowe.

First test with the secondImg as the queue-image creates the following results I've written down:

Inliner: "1672"
Number_of_Keypoints_FST: "368"
Number_of_Keypoints_SND: "1749"
Outlier: "77"
Ransac: "65"
Epsilon: "10"

Now I'm switching the order. The firstImg is the queue.

Inliner: "310"
Number_of_Keypoints_FST: "1749"
Number_of_Keypoints_SND: "368"
Outlier: "58"
Ransac: "65"
Epsilon: "10"

Questions: I'm confused that switching the order of the images affects the results. Should the matched keypoints be identical for the matched ones when switching the order of the images? Should the inliner be identical?


texts: Linearity with Difference Period in Fourier Transform

texts:
I know that a system is linear if it satisfies

$$\mathscr{F}\{ a\,x(t)+b\,y(t) \} = a\,X(\omega)+b\,Y(\omega)$$

for Fourier transform, $X(\omega)\triangleq\mathscr{F}\{x(t)\}$

But what if $x(t)$ and $y(t)$ don't have the same period?

Can we use same linearity condition in Fourier transform?


texts: How to Linearise an Input Image before applying Color Correction?

texts:
I'm doing color correction on Raw images taken from Camera directly. I want to make my Camera reproduce same colors as my Target camera.

Here are the steps I'm following:

Take Target Color values by Capturing a Macbeth Color chart by using Target Camera & Take set of 24 Color Palette's.
Take the Input Color Values by Capturing the same Macbeth Color chart by using Reference Camera & Take set of 24 Color Palette's.

Now if i Compute Color correction directly at this stage, all the White colors are appearing Rose or Not correct.

So I'm applying Gamma correction as follows:

Split the Three channels. For each Color find the Gamma for that color palette by using this formula.

float Gamma_R = log10(Target_R/255.0))/(log10(Input_R/255.0);

Error Cases:1

If the Input_Channel Value == 255 make it as 254

If the Input_Channel Value == 0 make it as 1

If gamma_Channel > 3 or gamma_Channel < 0.2 make gamma_Channel as 1.

Average the Values of Gamma_R for all 24 colors & make this as a
gamma of that channel.
Apply Gamma for Each channel using the Gamma computed for Each channel using the formula.

Corrected_R = 255 * (Input_R/255)^(1/Gamma_R)

My Problem:

How to make sure the Gamma Values that I've computed are correct at this stage? Also correct me if I'm doing it wrong.
What to do when the Following Error case: 1 happens?
After applying Gamma correction & Color Correction using 3x3 Matrix (Assume I'm doing that correctly) still the colors are not reproduced correctly.
If i encode the gamma by using same Gamma Values from Previous stage using this formula, again the Colors are not reproduced correctly.

Corrected_R = 255 * (Input_R/255)^(Gamma_R)
So do i need to compute the Gamma again from the Color corrected output?

Any suggestion or References is greatly appreciated!


texts: Decimating a discrete time cosine function

texts: Let $x[n] = \cos(wn)\ \ $ be our signal. Let the other signal be $x_2[n] = \cos(2wn)\ \ $ . The spectrum of $x[n]$ and $x_2[n]$ are well defined and consist of impulses at the relevant frequencies with weight $\pi$. However, if we interpret $x_2[n]$ as the decimated version of $x[n]$, the spectrum of $x_2[n]$ must have impulses with weight $\frac{\pi}{2}$. I have the explanation but I would like to hear the alternatives, what is the best way to resolve  these seemingly contradicting facts?


texts:
I apologize in advance if this question is not suited to this forum.

I have a DSP background from school(graduate EE courses). I am good on theory, but have never worked on a real project. After graduating and spending 2 years working as firmware engineer, I come to realize that I want to pursue DSP and algorithms for my career rather than pure firmware coding. I want to develop a portfolio of projects that can show the future employer that I can do this stuff. But I am a little short on ideas. Do you have any project suggestions?

Thank you very much.


texts: Calculating pre-history of recursive filter from state space representation when optimising for initial z

texts:
For a recursive filter, suppose a set of $b$ and $a$ coefficients have been calculated. Assume a state-space representation for which an initial set of $z$-values have also been calculated as in Python's lfilter_zi function.

The author of this post demonstrates code that applies the filter using the state-space model.

We know that an equivalent filter representation is the recursive one

$y[n] = \sum_{i=1}^{N} a[i]y[n-i] + \sum_{i=0}^{N} b[i]x[n-i]$

My question is as follows: what are the pre-history values of $x$ and $y$ that, when used to cold-start the filter using the recursive formula, yield the same results as the application of the state-space model?

For example, for a 2nd order filter, I want to plug in the values of $x[-2]$, $x[-1]$, $y[-2]$ and $y[-1]$ in the recursive formula above for n=0 and n=1 and obtain the same results as when using the optimal $z$-values from the state-space model.


texts:
A radio signal is sampled at $f_{s}=5512$ Hz, and $n=2^{14}=16384$ points are used in the FFT to display it as a spectrogram.
The frequency resolution is $\frac{f_{s}}{n}=0.33642578125$ Hz/point, displayed in the spectrogram containing $\frac{n}{2}=8192$ pixels (Nyquist).
Assuming there is no DC offset, the first pixel displayed is $f_{min}=0$ Hz, the last corresponds to $f_{max}=\frac{f_{s}}{2}-\frac{f_{s}}{n}=2755.66357421875$ Hz.
So far, I understand.
But, what's the precise bandwidth?
Is this $\frac{f_{s}}{2}=2756$ Hz, or rather $f_{max}-f_{min}=2755.66357421875$ Hz?


texts: From Fourier (k space) to wavelet domain in MRI sensing

texts:
In compressed sensing MRI (cSENSE MRI) technology the idea seems to entail sampling from the Fourier domain (k space) in a way that, when transformed to the wavelet domain ("sparsification"), the sparsity is maximized.

The inverse recovery problem becomes exact provided that $\mu$ is small:

$$\mu \left( \mathcal F W^\top \right)=\max_{i,j}\vert \langle W_i, \mathcal F_j\rangle \vert$$

i.e. the dot products of the column of the Fourier transform and wavelet transform are minimal ("mutual incoherence").

The idea seems to stem from this paper by Candes, Romberg and Tao.

In the wavelet domain the coefficients include scale and translation, while in Fourier space the coefficients belong to different frequencies without temporal support.

I would like to confirm that there is indeed a triple step: First (partially) filling in k-space with Fourier coefficients; second, randomly sampling these coefficients; and third, transforming them to wavelet space in each MRI image acquisition (as opposed to being an assumption based on the nature of MRI images).

And if this is the case, how does the step from Fourier to wavelet takes place (a reference would be OK).

Or, contrarily, whether the signal is primarily analyzed as wavelets?


texts: How to do De-Houghing of a Hough transform'ed Image?

texts:
I'm working with code found at Rosetta Code for creating a Hough transform. I now want to find all the lines in an image. To do so I need the ρ and θ values of each of the peaks in the Hough space. A sample output for a pentagon looks like this:

How can I find a single [θ,ρ] coordinate for each of the 'hot spots' visible in the Hough space?


texts: How the select the representative frame when detecting vehicles

texts:
I am detecting vehicles from the video (taken with a camera), and the vehicle detection part works fine.

My problem is that I should always extract only one representative frame per vehicle, and not the whole sequence.

As an example, if I have a few seconds of video (and possibly multiple vehicles), my approach currently extracts multiple frames (for every vehicle). Instead, as an output I would like one representative frame in which the vehicle is whole, i.e. well visible.

I have already implemented an entropy technique on it which makes it better but am still getting too much frames of the same vehicle. Could you suggest a technique which I can use to extract only one frame which contains the whole vehicle and ignore all other frames which contain that (same) vehicle.

I use the following entropy technique in my problem

double entropy = 0.0;
for (int i=0; i<histNorm.rows; i++)
{
    float binEntry = histNorm.at<float>(i,0);
    if (binEntry != 0.0)
        entropy -= binEntry * log2(binEntry);
}


texts: Principal Component Analysis as a step between STFT and MFCC

texts:
For computing Mel Frequency Cepstral Coefficients you can use already calculated STFTs as a basis and perform the Mel frequency mapping on it.

My question: Does it make a significant difference if I calculate the STFTs, perform a PCA-transform on them and then calculate the MFCC compared to computing the PCA-transform at the very end on the MFCCs?

In both cases we have two dimensionality reductions (PCA and MFCC) and I'm not sure weather the order makes a difference.


texts: Discrete Fourier style transform with non-sinusoidal/arbitrary waveform components

texts:
A discrete Fourier transform will produce frequency samples sufficient to recreate the original time sampled signal from sinusoidal waveforms.
That works great as long as the original signal is actually composed of sinusoidal waveforms. Let's assume that we know a signal is composed of non-sinusoidal waveforms, under which assumption its decomposition may be simpler, such as triangular waveforms. Does a discrete Fourier style transform exist which will decompose the sampled signal into triangular, rectangular or other waveforms? Does such a transform exist for any arbitrary waveform?
For example, here is the FFT response of a 150 Hz pure sin wave sampled at 1024 Hz (courtesy numpy and matplotlib):

Here is the FFT of a 150 Hz pure sawtooth, again sampled at 1024 Hz:

It's desirable here to have a similarly "clean" transform using the assumption that the signal is composed of sawtooth waveforms. So my question, does a transform exist that would decompose the signal into some non-sinusoidal waveforms?


texts: Detect the Beginning and Ending of a Signal

texts:
I have a signal that looks more or less  like square wave signal with different phases (Like in the picture below).

However, the signal fluctuates a lot, we can see in the image below:

At first, I want to  detect the beginning and ending of the signal.
Second, I want to detect the start and the end of the ups and downs.

Is there a simple way to do it? I am new to signal processing domain

I am using pyspark since I have a lot of data, but any help on how to start is welcomed.

Thank you


texts:
Apologies if this makes no sense or seems rather simple, I'm a Matlab newbie.

Basically, I've got a .wav audio waveform which I am resampling to various different frequencies and then trimming the length to a second (44100 samples). However, when resampling to certain tones, such as the note C3 or A5, I'm not left with enough samples to trim the note to an acceptable length.

Obviously this means that the sample needs extending, which I've worked out can be done taking a sample period s1 to s2 and adding that onto the end as many times as is necessary.

extendedSignal = x(1:s1-1);

for i = 1:4
    extendedSignal = vertcat(extendedSignal, x(s1:s2));
end

extendedSignal = vertcat(extendedSignal, x(s2:44100));

However, when it comes to finding a continuous period in which to sample, I cannot seem to find one for any note which does not sound like it is being repeated. Am I missing something really obvious?