texts: Getting sound pressure level (SPL) for a frequency band from Welch's derived PSD

texts:
I'm using Welch's method (pwelch, in octave) to get the PSD of a recorded sound signal (i.e. pressure level vs. time). I then took the PSD data and took the sum of the amplitude in each octave band (1/3 octave). My question is: how can I convert the PSD data to a sound pressure level (per band, or per frequency if I do it before 1/3 octave sort/sum). For an overall SPL, I can take the pre-PSD estimation data and calculate via the equation Lp = 20 log10(Prms/Pref).

Any suggestions?


texts: Proving that FFT (and by extension any correlation based estimator) is unbiased

texts:
I have a pretty simple question which i should have been able to answer. Just wanted to check if people here have a better solution to it.
I am trying to establish that FFT as a frequency estimator (imagine a single sinusoid + AWGN model) is unbiased. I guess we can agree that the answer to that depends on the number of grid points on the frequency axis or put in a different way it depends on how many zeros you pad to your signal. And no matter what, you will always be limited by your grid resolution. There are many papers out there which sort of interpolate using samples around the detected peak to refine the estimate and that pretty much is unbiased. This question can then be extended to other sort of correlation based estimators such as range estimation for a radar wavefrom etc.
I have a similar setting that i am working on, skipping the details for brevity but the signal model and algorithm are essentially what i described above. I plotted the resulting errors and they do come out to be unbiased. I was just wondering if i can somehow prove mathematically that the estimator is indeed unbiased. The reason i want to do that is because i eventually want to compare with the CRB.
If any papers/lectures or ideas come to mind, do share them. Thanks.


texts: interpolate 1D signal depending on 3D parameter space

texts:
I have a 1D array of data d(x,p) in which a number of "bumps" or "dips" appear and/or move in the spacial dimension x, depending on 3 parameters p = (p1, p2, p3). I know d(x,p) on a regularly sampled grid in the 3-D parameter space, but now want to interpolate between those samples for arbitrary values of the parameter vector p.

I know how to do a 3D interpolation of 1D data, however this does not really do what I want: suppose a bump moves to a different position upon changing one parameter, but does not change it's amplitude, interpolating the signal would just cause a lower amplitude and wider width of this bump.

I tried fitting the bumps and then interpolating their parameters, however this also does not work since both the number of bumps, their position and amplitude changes.

So I thought that there might be ways of interpolating in Fourier space, however I have very little knowledge about this kind of signal processing. Are there standard methods of achieving this or could someone point me to relevant works or at least keywords to look into?

Here's some made-up data that should illustrate the problem:


texts: Is there a relation between abs(fft(u.^2)) and abs((fft(u)).^2

texts:
Is there a relation between
(1) abs(fft(u.^2))
(2) abs((fft(u)).^2
where u is a general variable.
I have accidentaly found that if (1) has a power law +2, (2) has a power law +1. Can someone help me?
Thanks


texts:
My scope has 200 MHz bandwidth and can sample at 4GS/s

I'm sampling a 138 MHz RF signal at 500 M Sa/s storing 14 M points.

When I plot the spectogram in matlab I get the following
spectrogram(Data,12000,8000,[],fs,'yaxis')

Finally here is the frequency spectrum

While the signal does show up at 138 MHz, there appear to be aliases at 50 MHz, 190 Mhz and complex conjugate aliases at 86 MHz, 172MHz and 224 MHz. Does anybody know why this could happen?

edit:


texts: Meaning of tight frame in the context of texture synthesis

texts:
I am currently reading "A parametric texture model based on joint statistics of complex wavelet coefficients", by Portilla and Simoncelli. In this article, they use the "steerable pyramid" to represent the image. This uses a set of filters to decompose the image in a set of something around 16 subbands in the frequency domain, varying in spatial frequency and orientation.

Problem: They mention (page 55) that their decomposition forms a "tight frame". I know the definition of a tight frame in general linear algebra setting, or say orthogonal wavelet decomposition. However, in those cases the set of vectors that form the frame has cardinality of at least the dimension of the space considered. It seems highly unlikely that the space of images they want to consider is at most 16, and thus highly unlikely that this is what they mean. Can anybody tell me what is meant by the notion of a "tight frame" in this case?

Background: I am a mathematician that recently started working in neuroscience and now reading the mentioned article, published in the international journal of computer vision. Therefore the problem is probably just a matter of different languages used in different fields.

Thanks a lot!


texts:
I have a sinusiodal signal x in time domain.I have done the convlotion as followings using matlab:
y=conv(x,x)  and y=conv(x,fliplr(x)).
Observing the result by Matlab, both of the above signal giving the same result in frequency domain. Why do you think is that?

P.S: FLIPLR(X) returns X with row preserved and columns flipped
    in the left/right direction.

X = 1 2 3     becomes  3 2 1
    4 5 6              6 5 4


texts:
I need to detect a short (10ms) VHF pulse that occurs every 1.5s at a specific frequency (e.g., 150.20 MHz).

I have an SDR that support IQ sample rates up to 10MSPS but due to computing constraints 2.5MSPS is probably the max I can handle (currently working with 1 MSPS with 150 MHz as the central frequency).

I am new to DSP and started simple: Collect 500ms of samples, do a standard DFT using python's mlab.magnitude_spectrum(), throw away everything under an empirically defined threshold, followed by a simple peak finding routine. Then repeat with the next 500ms.

That works quite well but as the pulse strength drops (move away from the pulse source) this fails to find anything.

I have searched around a bit on how to tackle this and come across suggestions around measuring kurtosis, averaging spectra, autocorrelation, etc. If anything averaging seemed to make things worse but given I'm new to the area some guidance would be appreciated as to how to best tackle this as well as a good way to set the threshold. Is there something in the time domain I can do? Python (pseudo) code would be a plus.


texts: i am trying to understand how does the image rotation occur? i studied and came to know that we just multiply the image coordinates with rotation matrix which gives us the value of new coordinates. ok that is fine i get this. But main thing which is confusing me is the pixel values. i read that the pixel coordinates which we get after rotation can be non decimal hence interpolation is needed to get value at that point for this we rotate image back to original coordinate and find value by interpolation. But if we will rotate it back to original coordinate the the pixel will regain its original coordinate and at that coordinates value w are already having the value then why to do interpolation.please somebody clear my confusion.


texts: Resources for the Basics on Signals and Electromagnetic Waves in the radio frequencies

texts:
What are some good books and videos, to learn about radio waves?
For someone (me) who has never studied Math, Physics, or anything people generally consider ‘technical’.


texts:
I have long been faced with the confusion regarding entropy and would be obliged if the following are answered in less technical jargon. Following the link Different kinds of entropy raises the following questions

Entropy- It is desired that the entropy of the system be maximized. Maximizing entropy means no symbol is better than the others or we do not know what the next symbol / outcome would be. However, the formula states a negative sign before the summation of the probability logarithms. Thus, it means we are maximizing a negative value!! Then if an original raw signal is quantized and the quantized information's entropy is calculated and found to be lesser than the original entropy would imply loss of information. So,why do we want to maximise entropy since it would mean that we are maximizing the uncertainty of the next symbol whereas we want to be certain about waht the next occurence of the symbol would be.
What are the differences between Shannon's entropy, topological entropy and source entropy?
What exactly is the significane of Kolgomorov complexity or Kolgomorov entropy. How is it related to Shannon's entropy?
What information does mutual information between two vectors convey?


texts: Difference between cross correlation definitions (signal processing / random processes)

texts:
I'm trying to understand the similarities/differences between the following definitions for cross correlation)
Signal Processing: (or it's discrete equivalent)
$$R_{XY}(\tau) = \int_{-\infty}^{+\infty} x^*(t) y(t+ \tau) dt$$
Random Processeses$$R_{XY}(t_1, t_2) = E\begin{Bmatrix}
X_{t_1} Y^*_{t_2}
\end{Bmatrix}$$
Are these related?


texts: What to do for the right FFT in Matlab (two peaks and incorrect amplitude)?

texts:
This is the FFT I'm using. But two peaks (did not expect the one in the beginning) are occurring and the amplitude is not as specified (expecting a value of 2). Any help is much appreciated!

Amp = 2;
freqHz = 10000;
fsHz = freqHz*2+1;
dt = 1/fsHz;

sine = Amp*sin(2*pi*freqHz*(0:dt:1-dt));

transform = fft(sine,fsHz)/fsHz;
magTransform = abs(transform);

faxis = linspace(0,fsHz/2,fsHz);
plot(faxis,fftshift(magTransform));
xlabel('Frequency')


texts: How do you derive the new mean and variance when two gaussian functions are multiplied together?

texts:
I'm trying to understand the contents of this web-page which suggest that multiplying two Gaussian functions together produce a new Gaussian function with a new mean and variance.
Say I have the two following Gaussian functions:

The web-page state that the new mean and variance to be:

I start by multiplying the two functions, but I am unable to get the result into Gaussian form such that I can extract the new mean and variance. Any tips on how these values were calculated?


texts:
Since my algorithm needs to detect long, coherent edges, the estimated local magnitude P for each pixel is normalized with respect to the magnitudes in it’s local neighborhood W (i.e., a 5x5 window).

$P_{normalized} = \dfrac{P - \mu _w}{\sigma_w}$

where $\mu _w$ and $\sigma_w$ denote the average and the standard deviation of edge magnitudes in the pixel neighborhood $w$.

This is given in the paper "GradientShop: A Gradient-Domain Optimization Framework for Image and Video Filtering" by Bhat, and I need to implement the algorithm described there. I am not getting, what does normalizing the gradient mean exactly ?


texts: Why Isn't the ML Estimator (MLE) in MIMO Spatial Multiplexing Obtained by the Least Squares Solution?

texts:
In the simplest scenario of MIMO spatial multiplexing: $$\mathbf{y} = H\mathbf{s} + \mathbf{n}$$
where:     $\mathbf{s}=[s_0,s_1,...s_{M-1}] \\\mathbf{y}=[y_0,y_1,...y_{N-1}]$

$\mathbf{n}=[n_0,n_1,...n_{N-1}]$ is zero mean complex normal independent Random vector

$M$ transmit antennas, $N$ receive antennas where $N\geq M$

The maximum likelihood solution is:
$$\tilde{\mathbf{s}} = \arg\min_{\mathbf{s} \in QPSK^M} ||\mathbf{y}-H\mathbf{s}||^2$$

Now, my textbook says that due to the fact that in this scenario $H^{*}H$ is not a diagonal matrix, the least squares solution followed by proccesing is not optimal, but I can't see why this is true.
To my understanding, denoting the LS solution as  $\hat{s}=(H^{*}H)^{-1}H^{*}\mathbf{y}$
I can write:
$$  ||\mathbf{y}-H\mathbf{s}||^2 = ||\mathbf{y}-H\hat{s}+H\hat{s}-H\mathbf{s}||^2 = \\||y-H\hat{s}||^2 + ||H(\hat{s}-s)||^2 + 2Re[(y-H\hat{s})^{*}H(\hat{s}-s)]$$

Now, by substituing the LS solution the $Re$ part vanishes.
And the problem reduces to
$$\arg\min_{s\in QPSK^{M}}|\hat{s}-s|^2$$
So I don't get why the ML is not optimal is this case.
Can anyone please clarify?
Thanks

edit: screen shot from the book


texts: Periodically-varying phase between transmitted sound and received sound

texts:
I have a function that can find the phase difference between two signals, this function uses Hilbert transform to do so.
Here is the python code for it

def phase_shift(carr_wave, rec_wave):
     assert len(carr_wave)==len(rec_wave)
     carr_comp = hilbert(carr_wave)
     rec_comp = hilbert(rec_wave)
     c = np.inner( carr_comp, np.conjugate(rec_comp) ) / math.sqrt( np.inner(carr_comp,np.conjugate(carr_comp)) * np.inner(rec_comp,np.conjugate(rec_comp)) )
     phase_diff = np.angle(c)
     phase_diff = abs(phase_diff)
     return phase_diff

I am inserting a sample code of how this works and how many samples it need to be near accurate

f = 20000
offset = 0.25
samples_in_bit = 55
t = np.linspace(0,1,44000)
signal_1 = np.sin(2*math.pi*f*t)
signal_2 = np.sin(2*math.pi*f*t + offset)
phase_difference_1 = phase_shift(signal_1,signal_2)
phase_difference_2 = phase_shift(signal_1[:55],signal_2[:55])
print phase_difference_1,phase_difference_2

output :
0.249981325271 0.249951020276

I am doing BPSK at a frequency $= 20000\textrm{ Hz}$, duration of one bit $= 0.00125\textrm{ sec}$, with a frame rate of $44000$ per second. I am modulating 1800 bits which are all zeros.

Now I plot the phase difference of this modulated wave and actual carrier wave (Both will be same as I am sending all zeros) with perfect synchronization of both waves, the plot looks like this, which is obvious.

NOTE: I am finding phase difference at each bit (from start to end sample in one bit duration) and plotting against bit number on $x$-axis and phase on $y$-axis.

Now let's say if there is some $20$-$30$ samples of miss synchronization between carrier and modulated wave there will be a an offset phase where the plot looks like this. Which is okay that its phase difference is constant through out all bits.

After modulation, when I transmit the sound and record it (Using a good mic). I found that the phase difference plot between the transmitted wave and the recorded wave looks like this,almost a periodic function.

And this one is the same plot, but I have plotted absolute of phase differences.

To exactly know the correct sync between recorded and transmitted waves I am using cross correlation techniques which give a correct synchronization correct to $30$ frames(samples). So even if synchronization instant is wrong by $30$ frames (I am having $55$ frames in on bit duration)  I should at least get a constant phase plot just like it did on my modulated wave (before transmitting).

My question is why is the phase difference between transmitted and recorded waves constantly changing ? And that too it is not a random change. It is kind of periodic.
Does this normally happen in all signal communications? or is there any problem in my transmission?
Or can I do any filtering after I receive the signal?

EDIT :

As the problem for this is found to be frequency offset between sent and recorded waves.

Is there a method or an idea to estimate the frequency from the recorded wave, which can be used to find the phase using Hilbert function.?

I am including the frequency spectrum of my sent and recorded wave as shown below.


texts: Filter coefficients for colored noise Voss Algorithm

texts: If I use a generic filter for generating colored noise like pink,brown,white then how do I modify this statement and how do I know what are the coefficients to be used in AR model for different noise.Signal is white gaussian noise and y gives different noise based on the value of a.


texts: How Could I extract the signal in the below Image

texts:
This is my first post so sorry if I make any mistake.

The entire signal length in the screen shot below is 10 second.
I ran a first order low pass filter (digital) at 10Hz cut-off frequency, and I obtained the orange/yellow line in the picture. Please ignore the red and the blue line for now.

I wonder if there is an algorithm that I can extract all of the lowest values of the green plot in real-time and reconstruct/interpolate data along the way like the black line that is shown below?

The green signal (50Hz) is the DFT result of a raw 1kHz sine data that is sampled at 10kHz. I really hope that we can find a light-weight algorithm for the filter, since the entire code has to run on a microcontroller.

I would love to hear your suggestions. And let me know if anything is needed for clarification.


texts: Region of convergence for discrete system

texts:
I am confused by two questions on the region of convergence (ROC) when applying the $Z$-transform.

Consider a discrete-time linear system $$x[n+1] = Ax[n]+Bu[n], \qquad y[n]=Cx[n]+D$$ whose transfer function is given by $$H(z) = C(zI - A)^{-1}B + D$$ Is there a region of convergence for $H(z)$? Is it supposed to be the whole complex plane?
Also from the Woodbury matrix identity, $(zI-A)^{-1}$ can be written as $\sum_{k=0}^\infty \frac{1}{z^{k+1}}A^k$. How should I derive the ROC of the latter series? Does it matter if the ROC of series is different from the ROC of original transfer function H(z)?

If only considering causal sequences, is it possible to have two sequence with same $Z$-transforms, but with different ROC? To be specific, consider two causal sequences $\left\{x_1(k)\right\}_{k=0}^{\infty}$, $\left\{x_2(k)\right\}_{k=0}^{\infty}$, is it possible  denote the two $z$-transform as $H_1(z), H_2(z)$ respectively. If $H_1(z) = H_2(z), \forall z\in ROC_1 \cap ROC_2$, can I claim that $x_1(k) = x_2(k), \forall k$?


texts: One of the best ways to numerically integrate a signal?

texts:
I need to get position $x$ (integration) from velocity $v$. One could use 1st order Euler integration as

$x_{t+1} = x_t + \Delta v_t$

However, doing so leads to errors proportional to sampling time delta.
Do you know any more accurate solution please?


texts: Time domain and frequency domain analysis equivalence

texts:
Let us imagine an LTI system with physically realizable input (ruling out fancy mathematical functions and the concomitant complexities and paradoxes) completely known from -$\infty$ to $\infty$.  We want to calculate the output. We can analyse it in time domain using the linear constant coefficient differential equation and in the frequency domain using the Fourier representation of the signal and frequency response of the system.

The outputs calculated from both the methods should be identical. In time domain, we get (transients + steady state)  response of the system,  In frequency domain, the sinusoids are eternal and hence transients are generally neglected.  But since we are talking about the same system with the same input, the outputs calculated from both the methods should be the same.

Will they be identical in spite of neglecting the transients in the frequency domain?


texts: Is it possible to generate depth map from graycode sequence without using printed calibration patterns?

texts: In a camera-projector system, if only the camera's intrinsics are known, is it possible to generate the depthmap fully automatically through a graycode sequence without any manual calibration process (like using a printed calibration pattern then move it away etc.) ?


texts: Presenting an all-pass filter of a high order as a cascade of all-pass filters of the 1-st order

texts:
Are all-pass filter of the n-th order just a cascade of N all-pass filters of the 1-st order?

There are some nice statements regarding zeros, poles, and gain for the 1-st order. Can those statements be applied to the n-th order all-pass filter too? For instance, the Pole is a reverse complex conjugate of the Zero.

Obviously, cascade of N all-pass filters of the 1-st order will produce the all-pass filter of n-th order, so the 'necessary' condition is true. But is the 'sufficient' condition true as well?

There are some edge cases of course, for example, if we consider two kind of all-pass filters with abs(H(exp(iw)))=1/2 and abs(H(exp(iw)))=2 then their multiplication (i.e. cascade) will give a true all-pass filter. But to me it's still an edge case. I'm more interested in non-constant Amplitude Frequency Responses.


texts: How to choose a windowing function when performing a Chirp $\mathcal Z$-transform?

texts:
I've got some data sampled at 32 Hz. The only interesting stuff is happening between about 0-8 Hz and I've got to do real time batch processing and detect events with a delay of no more than a second. So that means I can only perform FFTs with 32 bins or less, to remedy the lack of resolution I'm using a Chirp $\mathcal Z$-transform (czt in Matlab).

I was wondering if there were special considerations to be made when you choose a windowing function on the original time data before when performing a CZT instead of an FFT?


texts: Sketching a tri signal with horizontal shift

texts:
I'm unsure how to sketch a tri signal with horizontal shift.

Eg , let's say we want to sketch tri((t-4)/4). Can I reduce it to tri( (t/4) - 1 )? This would mean sketching a tri(t/4) signal and shifting it to the right by 1. Or should it be a tri(t/4) signal shifted to the right by 4?

I've found out on wikipedia that the rect signal is defined as rect((X-Y)/T) where you shift the center of the signal by varying Y. But I'm not sure if sketching tri applies the same concept.


texts: 1D filter for speed input and noisy position sensor

texts:
I was tasked with designing a filter to smooth out a translating table. The system, we command a speed and a desired position, and measure it with two noisy position sensors. The commanded speed is also not always perfect, but it's usually within 10%. The slide can only move left and right.

My initial thought was to simply use a moving average (via circular buffer), and grow / shrink the array based on the speed. For high speeds, I would have a smaller array (10 or so measurements), and for lower speeds, I would have a large array (maybe 40 measurements). When stopped, I would have it grow to 100 or so. The reason for this is to decrease the phase offset at high speeds where resolution is not as important, and decrease the noise at lower speeds, when accuracy is more important. However, this doesn't really take into account the speed, so I feel like I'm throwing away useful information.

I'm figuring there has to be a better way. Any thoughts on how I would approach this problem more elegantly?

Just so we're on the same page, "high speed" is only about 3 inches per second, and I am sampling at 20Hz (but each time, I get two measurements). There is a ramp up/ ramp down portion of the translation such that whenever the table approaches the desired position, it slows to a crawl of 0.5 inches per second.

Thanks!


texts: Chosing polyphase interpolation coefficients

texts:
I'm writing code for a 3x polyphase interpolator using a total of 9 coefficients. This is organized as 3 parallel FIR branches with 3-taps each that are summed at the 3x sampling rate.

The problem I have is the impulse response looks correct (I get my 9 filter coefficients as output) but when I use a real signal that is a smooth Gaussian area near the peak has 'digital noise'. I have confirmed by debug printouts the code is calculating what I want it to. It seems the combination of the input signal & the coefficients really does generate a non-monotonic output.

as reference I used the info at this site: [polyphase filter faq][1]

What am I missing? How do I avoid this? My coefficients are:
[4066   6261    8133    9391    9834    9391    8133    6261    4066]

The inputs sample rate is 100Mh and the output is 300Mh. I choose the filter coefficients to have a bandwidth < 50Mh.

Adding followup information. I looked at the step response and it oscillated every 3 output samples. The output to a step should be the sum of the coefficients in a branch. I noticed the 3 coefficients of the 3 branches did not all add to the same value. By deduction, I changed the coefficient 6261 to 5878. Now the coefficients are:
coeff = [4066   5878    8133    9391    9834    9391    8133    5878    4066] and the response is much better.

See new plot:

Thanks,
Brian


texts: What is convolution of two sine waves (tones)?

texts: Convolution of two sine waves (or tones as called in audio) is theoretically not defined as the integral is infinite. Taking finite duration windowed sine waves and doing there convolution computationally always contains a fundamental frequency equal to that of the lower frequency sine wave. I am not getting an intuitive understanding of this. Anybody who can share more understanding in this?


texts: how can the channel fading model (Rayleigh, or Rice, or Log normal distribution) infer the way channel estimation is done?

texts:
For a dense urban environment with rich scattering and no Line of sight (LoS)
path, the Rayleigh distribution can be used to model the fading.
Rician fading is a stochastic model for radio propagation anomaly caused by partial cancellation of a radio signal by itself.

How can the knowledge of the statistical distribution of the noise impact the way the channel estimation is done ?

I need to model a type of channel with the knowledge of its statiscal distribution: for that i use a Discrete FIR Filter block, and then impaired with the AWGN block.


texts: Why the time step series affect the noise in FFT results?

texts:
I used Matlab to calculate a sinusoidal function:

Cy2 = 0.0226809015035037 - 0.235834694520537*sin(5.06307730580011 - 2*pi*2*t);

If I choose a set of Matlab generated data:

dt=5.305164999999997e-04;
t=[dt1:dt1:245];

Then conduct FFT on y(t), I have a nice and clean result (black line):

However, if I use my experimental data's time steps 't_exp' (download here)

t=t_exp %t_exp is the time steps corresponding my experimental data

I get much more noisy data (black line):

I just wonder how come the extra noise while the function is exactly the same, and only a different time step series are chosen.

Conversely, I got noisy experimental data with such time step series (t_exp). Is it possible to make it cleaner by modifying the time step series? I almost mis-interpreted my exp data due to the noise.

Appreciate any help.

==========================

Appendix:
my code to do FFT:

function [f, P1]=fftzp2(t,y)
zpf=1;
%zpf=2;
L=length(t);
dt=t(20)-t(19);
Fs = 1/dt;
Y = fft(y,zpf*L);
%Y = fft(y,L);
P2 = abs(Y/L);
P1 = P2(1:ceil(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:ceil(L/2))/(zpf*L);
%% energy conservation check
Ent=sum(y.^2)*dt;
df = Fs/L;
Enf=sum(abs(Y*dt).^2)*df;
if Ent-Enf>1e-9; disp('energy not conserved'); pause;end
end


texts: Calculating the peak power of a signal on python

texts:
As far as I've researched, the energy and power of a given (discrete) signal are given by
$$E = \sum_n \left|x_n \right|^2$$
$$P = \lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_n \left|x_n \right|^2$$
Where N is the lenght of the given signal.
Im working with a signal embedded in some non-gaussian noise, and I want to calculate the ratio of the peak power of the signal and the power of the noise (see label of Fig 2 on https://arxiv.org/pdf/1701.00008.pdf).
Now, the power of the signal I already know how to compute, but I have not been able to find a precise definition of what peak power is, and I have not found a library on python to help me on that.
Is there a library on python that does this, or is there a concrete definition of peak power I can use to calculate it?
Thanks in advance!


texts: How to determine if a pixel is inside a contour?

texts:
I have following picture

which I extract the contour from. The contour looks like following

I am wondering, for any pixel in the image, how do I determine if it is inside the contour or not? Say "inside" here means any pixel in the body.

If possible, I'd prefer simple, direct, and elegant algorithm. Thanks.


texts: Image cleaning - How to remove line art from a colored image?

texts:
I have a set of line art drawings, which I have cleaned and colored (flood fill). I want to remove the original line art from the images and to leave only the colors.

This image illustrates pretty much what I want:

All of my line art is of a solid black, without anti-aliasing, and at high resolution. The resulting image has all black pixels replaced by nearby colored pixels. Colored pixels from the original image are left untouched.

I made a rough Python program to generate the "after" image. Notice that the colors are offset to the right and bottom of the original black areas (which causes distortion of my shapes when the line art has varying width). Also, my program is very slow; hence, I'm looking for a smarter solution to this problem.

Edit: An algorithm similar to those used to compute the morphological skeleton of the image's black pixels will greatly improve this operation.


texts: I'm working on a multiuser wireless system simulation and have run into something regarding sampling rates. From everything I've known for a matched filter receiver, the correlation step is performed followed by symbol rate sampling. Say the symbol rate is $T_s$. My question is: do we ever sample at higher than the symbol rate? What if I was able to sample at $\frac{T_s}{5}$? It seems to me then I'd have five samples of the same symbol, which is in the same spirit of repetition coding? Is there a difference between repetition coding and just sampling faster?


texts: Could someone advise me on some module for making the TDoA position system for my research work? I need to have measurements with an accuracy of a few meters. It equals a bandwidth of about 20 MHz. I guess there are some modules of 2.4 GHz with such bandwidth and given propagation time and also raw output without any high-level stacks like OSI. The only output of the threshold device after the matched filter maybe. I also need it to have a possibility to set my own partial signal (for example I want to use M-code).


texts: Implementing short-time Fourier transform

texts:
I have a $512$ sample frame, lets call it X. The frame is further divided into four sub-frames of $128$ samples each. Each of these sub-frames have a different frequency.

For STFT, I use MATLAB's gausswin(N) function where N is equal to $128$ in my case to generate a Gaussian window, let's call it W. I zero-pad my windows W to make it equal to $512$ samples.

I now multiply my original frame X point by point with W to get y. I then take FFT of y to get Y. It works till here.

My question is, naturally the next logical step should be to shift my zero-padded window W to right by 128 samples and repeat the multiply and FFT procedure again.

But how do I accumulate these results?
Concretely, how do I convert it into a time-frequency plot?


texts:
I'm using gist.github.com/sixtenbe/1178136#file-peakdetect-py to find maximum peaks in my spectrogram.

I used

_max, _min = peakdetect.peakdetect(array, lookahead=300, delta=0.3)

I get the following peaks. x_axis = time, y_axis = frequency

As you can see in the image, the peaks are not evenly spread out. All the points are skewed towards frequency axis. My array is in this format: [1, 2, 3 , 0 , 1, 3, 4]. It's 2D time+frequency array flattened into 1D array.

What can I do to show the image like this where the peaks are spread out evenly.


texts: Isolate recurrent pattern vs. evolving background sound

texts:
Here is a sound example for what follows.

Let's assume we have a signal

$$s(t) = r(t) + e(t)$$

where:

$r(t)$ is a signal which is recurrent with a given period, i.e. in my example $r(t) = $ the exact same snare sample, which is played every 1 second, starting at 00:00:00,500
$e(t)$ is an evolving sound (in my example it's the synth pad soft sound)

I assume the signal is mono, and that no compression or other master effect has been applied.

Question:
How to recover $r(t)$ from the mix $s(t)$, by analyzing which "pattern" is recurrent in $s(t)$? (the period 1 second is considered as given)

In one word, how to isolate the snare sound from the mix?

I would like to do it by slicing the signal each 1 second:

slice1.wav
slice2.wav
slice3.wav
slice4.wav
slice5.wav
slice6.wav
slice7.wav

and analyze what part of signal is recurrent (= the snare) and which is not (= the synth)  among all slices ?

Is there a separation algorithm based on analyzing recurrent parts / non-recurrent parts? Using FFT or something else?


texts: What is the DTW result (applied to the MFCCs of an audio file) proportional to?

texts:
I have computed the MFCCs for two audio files and applied the DTW to both sets of coeficients for both files. Now each such pair the two files is described by one single value (computed by the DTW). Let's call this value Diff.

If the audio files differ more, then Diff is greater, which is great. (so Diff is proportional to the difference between the audio signals - great!).

My problem is that Diff seems to be also proportional to the length of the files (?). Therefore, for two similar long audio signals, Diff might be greater than for two different short signals.

First I thought that, as DTW computes the shortest path for two arrays, the final value should be proportional to the lenght of the two arrays (i.e. the number of frames for each of the signals). So I tried dividing it by the product and by the average of the two lenghts but the results was discouraging (too big for some files that did not differ so much, too small for others).

Then I thought that if it is proportional to the lenght of the files, then dividing it to the average number of frames (in which I divided the two signals for computing the MFCCs) would do the trick... but it didn't...

Finally, I thought that, as the MFCCs are computed from the result of the FFT, the final DTW result might be proportional to the sum of the FFT magnitudes, so I tried dividing it by this sum (of the spectrum elements, returned by the FFT algorithm). Again, the results were not what I hoped them to be...

Right now I'm out of ideas... What is the DTW result depending on, what is it proportional to?

EDIT:
Here are some values, after pichenettes' suggestion to try the geometric average of the lengths of the files:

same sound . avgLengths = 106000. Diff = 409. Diff / avgLengths = 38
same sound . avgLengths = 191000. Diff = 609. Diff / avgLengths = 31
diff sounds. avgLengths =  80000. Diff = 437. Diff / avgLengths = 55
diff sounds. avgLengths = 193000. Diff = 692. Diff / avgLengths = 35

38 > 35 :(

(the sounds from 2. and 4. are almost similar, and so are the ones from 1. and 3. In 1. the trailing silence for one file is a little longer)


texts:
I would like to know how $\frac T2$ came for time period in the image below.Because I think the frequncy is 1 so time period should be 1.Is that correct? or can anyone help me..

The article from where I got this image is here


texts: Can this be considered wide sense stationary?

texts:
I was discussing this problem with one of my classmates.
The picture shows a recording of the heart rate during before and after sleep.

Can the whole process be considered wide sense stationary? (I say yes because the mean is approximately constant, but how can I estimate the autocorrelation?)
Can the whole process be considered ergodic in autocorrelation and mean? (I'm not sure at all about this)


texts:
I have a filter design problem in which lowpass filters have to be designed with varying roll-off. E.g: 2.7 db/ocatve or 1.4 db/ocave. See picture below.

This is easy to achieve with a FIR filter, yet its computational cost is too large. Therefore, if at all possible, I would like to do this with an IIR filter. Sadly, most IIR filters have multiples of 6dB/octave rolloff.


texts: What is the difference between entropy and entropy rate and which one is important in information theory?

texts:
By entropy I understand the uncertainty or randomness. But if the uncertainty is high then what is its implication and advantage in transmitting information or as an information source or is high entropy preferred for processing information?

Then there is entropy rate, which I cannot understand due to heavy mathematics. All that I can follow is that entropy rate = entropy/ number of unique symbols.

Please forgive if these questions seem trivial, but it is quite hard to filter the answers from textbooks which are mainly focused on theoretical definition and not on how one can use entropy in real applications. Below are my confusions.

If the entropy of an information generating source is high, that is considered to be good. But what is the role of entropy rate?
Does entropy have a lower and upper bound or a range? Based on equi-probable symbols, we say then 0.5 is the max entropy as equi-probable symbols convey maximum uncertainty about the source. Is the range of entropy is 0--1?
What information is gained or loss when entropy is discussed from information theory viewpoint?


texts: Difference between delta modulation and differential PCM

texts:
Delta modulation encodes the difference between a sample and the previous sample in only one bit.
Differential PCM encodes the difference between a sample and its predicted value. Is this difference encoded in one bit or more than one bit? If encoding is done in more than one bit, what is the advantage of differential PCM over delta modulation?


texts: Gabor filters for digital image processing

texts: For the application of Gabor filters in feature-based (e.g. lines and edges in images) digital image processing, it is desired that the DC component is 0. Why is this so? What effect does it have mathematically if the DC component of the filter is not 0?  Thanks in advance!


texts: Scipy Butterworth highpass filter results in distortion and bloated output

texts:
Trying to use Scipy Butterworth filter to filter frequencies below 54 Hz from audio files. Audio files are stored as WAVs. Sample rate is 44100.
When I try using this filter the result is distorted (very). Also, I assume with fixed time and sample rate that input and output files should be exactly the same size. Input file is 441 Kib, but output is 1.8 Mib. But the duration is unchanged.
Problem is with the filter or with application of the filter. I'm opening and saving WAV files with scipy.io.wavfile and have tested without applying the filter -- what I get out is what I put in, as expected.
Here's my code:
import scipy.io.wavfile
from scipy import signal

# Scipy reads audio data from WAV files
rate, data = scipy.io.wavfile.read('./test.wav')
# rate is the sample rate, data is the data
# NOTE: Sample rate of my input is 44100
assert rate == 44100

cutoff = 54  # want this to be 54 Hz
nyquist = 0.5 * rate
normal_cutoff = cutoff / nyquist
order = 5

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html
# Create a filter and return coefficients
b, a = signal.butter(order, normal_cutoff, btype='highpass', analog=True)

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html
# Apply filter to audio signal
filtered = signal.filtfilt(b, a, data)
scipy.io.wavfile.write('filtered.wav', rate, filtered)

Have also tried filtered = signal.lfilter(b, a, data) for application of filter, and have tried analog=False. Behavior is similar in both cases.
What am I doing wrong? Does this have to do with the sample rate? Bit depth? I'm stuck.


texts:
I'm currently using MATLAB's fdatool for filter design. Using that tool, I can easily design different kind of filters. For example, let's take a band-pass FIR filter with 10-40 Hz passband, and 5-10 Hz and 40-45 Hz transition bands. Usually, I design the filter with the selection "least-squares", which, if I understand correctly, uses the aforementioned method to find the best impulse response according to filter spesifications. To actually filter the signal, I use the command filtfilt, which does zero-phase FIR filtering.

Now, an alternative way to implement the filter would be to take the FFT of my signal, set frequencies outside the range 10-40 Hz as zeros, and then take the IDFT.

Is there any practical/theoretical difference between these two approaches? Will the frequency responses (magnitude and phase) be the same?


texts: How to construct Neyman Pearson statistic to detect discrete signal in AWGN with unknown location?

texts:
I have a discrete datastream of length $10^5$ samples. Somewhere embedded in it is a rectangular pulse of unknown amplitude (prior probability of data containing signal is 0.5) and width 100 samples. I need to be able to detect this, so I think my  statistical model should be:
\begin{align}
\mathcal H_0 &: Y_i = W_i \sim \mathcal{N}(0,1) \\
\mathcal H_1 &: Y_i = W_i + A\big(H[p] - H[p-100]\big)\sim \mathcal{N}\bigg(A\big(H[p] - H[p-100]\big), 1\bigg)
\end{align}
where $H[p]$ is the Heaviside function, and $p$ is unknown, but can be distributed uniformly over the sample length. I need to construct a Neyman-Pearson statistic from this model, but I'm unsure how to proceed. Any suggestions? Thanks in advance.
Edit 1: I am looking for an approach that requires $\mathcal{O}(N)$ or $\mathcal{O}(N\ln{N})$ operations, so my first guess is I should be able to reduce the test statistic, or at least be able to evaluate it, in a convolution sum/integral over the data. Also, the prior of $1/2$ implies that roughly 50% of the time all such operations are in vain (especially the moving average approach), so I need an computationally effective way to detect this.


texts: Calculate Q factor of a Low Shelf and High Shelf filter

texts: I'm developing an app that has an EQ using eq10q plugin. My eq has 10 bands. I'm using peaking filters for the middle bands and a low shelf filter for the first band and a high shelf filter for the last band. I calculated the Q factor for the middle bands by dividing the middle frequency by the bandwidth desired for that band but I don't know how to calculate the Q factor for the High Shelf filter and the Low Shelf Filter. How should I do it?


texts: Is there a way in MATLAB to approximate an FIR filter using an IIR Filter?

texts: I am a student trying to understand how this is possible. I've done some research on the web with not many sources showing up with solutions to this problem that are usable. What techniques are used to give an approximation of an FIR filter using IIR filters with say the same order?


texts: Does FFT magnitude spectrum graph consist of spectral lines or points

texts:
Hi I am using zedgraph to plot an FFT of a signal using c#.The FFT spectrum gives coefficients vs the frequency bins.I am confused whether the plot should consist of spectral lines(ie verticle lines) corresponding to each frequency bins or data point values?
Thanks in advance.


texts: NF of attenuator is same as the attenuation value as the noise floor remains same but signal level gets attenuated? But NF of amplifier depends on the device as it increases both signal and noise?


texts: How to restore an image from its filter output?

texts:
I'm wondering if there is a way that I can reconstruct an image, given the output vector obtained from the convolution with a certain filter.

For example I convolve my image with a filter from MR8 filter bank, so I got a feature vector v. Which is :

v = img ** filter

so can I define a function that takes v and filter as an input, and restores the img with a certain error ?

img_restored = reconstruct(filter,v);

What have I tried

I searched for deconvolution, I found two things :

Wiener Deconvolution : But it only works for de-blurring (I guess)
deconv() function of MATLAB, but it takes two vectors as inputs, and gives one -complex- number as output.

Edit

Another approach could be storing all feature vectors for pixels (centers of sliding windows) and try to restore only a single pixel, then obtain an image by concatenating these pixels.

Thanks for any help !


texts: What is the significance of the DFT in SC-FDMA systems?

texts: Long term evolution (LTE) systems use single-carrier frequency division multiplex (SC-FDMA) in the uplink. In such a system, there is a DFT block and a sub carrier mapping block before the IFFT block. What is the significance of DFT and sub carrier mapping? And how do these blocks make the system single carrier?


texts: DFT zero-padding of signals starting before n=0

texts: If a signal starts before n=0, what part of the signal should be used to compute DFT after zero-padding? For example, x(n) = {1, 2, 3, 4, 5}, where x(-2) = 1 and x(0) = 3. If this signal is zero-padded to N=8, the new signal would be x'(n) = {1, 2, 3, 4, 5, 0, 0, 0}, where, again, x'(-2) = 1 and x'(0)=3. However, to compute the DFT of this zero-padded signal x'(0), should I use {1, 2, 3, 4, 5, 0, 0 , 0}, which starts with x'(-2), or should I use {3, 4, 5, 0, 0, 0, 1, 2}, which starts with x'(0) with the first few terms (n<0) folded over? In terms of the DFT X(k), I think its magnitude would be the same for the two approaches, but its phase would be different, right?


texts: Cross correlation coefficient using MATLAB

texts: Actually I want to calculate the cross correlation coefficient of any signal on comparing with all mother wavelets like Morlet, Daubechies using MATLAB.


texts: What is a proper chunk size for audio analysis?

texts:
Here's my code:

import pyaudio
import wave
import sys
import spl_lib as spl
from scipy.signal import lfilter
import numpy

CHUNK = 4096   # This is the sample size
               # math.pow(2, 12) => RATE / CHUNK = 100ms = 0.1 sec
FORMAT = pyaudio.paInt16    # 16 bit
CHANNEL = 1    # 1 means mono. If stereo, put 2
RATE = 48000

error_count = 0
numerator, denominator = spl.A_weighting(RATE)

pa = pyaudio.PyAudio()

stream = pa.open(format = FORMAT,
                channels = CHANNEL,
                rate = RATE,
                input = True,
                frames_per_buffer = CHUNK)

print "Listening"

while True:
    try:
        # read() returns string. You need to decode it into an array later.
        block = stream.read(CHUNK)
    except IOError, e:
        error_count += 1
        print(" (%d) Error recording: %s" % (errorcount, e))

    decoded_block = numpy.fromstring(block, 'Int16')
    y = lfilter(numerator, denominator, decoded_block)
    decibel = 20*numpy.log10(spl.rms_flat(y))
    print('A-weighted: {:+.2f} dB'.format(decibel))

stream.stop_stream()
stream.close()
pa.terminate()

The code is calculating dBA.

I set CHUNK as 4096 because it is approximately 100ms and I thought that's a reasonable length. I thought chunk sizes should be power of 2 because all the web articles I read had 1024. This threw an error about every 1 second:

[Errno Input overflowed] -9981

I thought it was a memory issue but when I ran "top" or "free -m" I had enough memory. Then I changed the chunk size as 8192(2^13) and 1024. I expected chunk size of 8192 would throw more error as it will have calculations with larger buffer. However the result was interesting. It didn't matter whether it was 256 or 8192. They were all small enough for Raspberry Pi 2. However, what was more interesting or even exciting was that small difference of chunk size would bring dramatically different results.

With a chunk size of 258, error happened for every loop when 256 was error-free.
I carefully made an assumption: Why does a chunk size should be a power of 2? So I tried with 4800 which is a tenth of 48000 and then it didn't throw any error.

So I'm really confused. What is it to do with chunk size? Does it necessarily need to be power of 2? Why did 256 didn't throw an error when 258 threw error on every loop?


texts: Converting Received RF packets to binary file that RF packets sampled by dongle rtl-sdr

texts:
Im using #SDR software that I transmit packets from cc1350 launchpad to the air and I want to validate that if those packets are valid or not by checking its CRC.
my problem is actually the step to convert the demulated data that ranges between -1,1 to binary 0/1.
I transmit packets on frequency 868MHZ, 100Kboud(same as bit rate) , deviation 100KHZ. it's GFSK modulation. so in order to do my own receiver, Im using matlab to build my own receiver and check if the received packet are valid by its CRC, so Im sampling the transmitted packets by a dongle on sampling frequency 20480000HZ , so I see the sampled (IQ RAW DATA) in matlab shows my the packets like a peaks!
attaching down photo of my sampled data and it contains packets (there're peaks when there's packet )

So now Im doing my demodulator that its outputs a data ranges between [-1,1] , the outputs for instance could be a vector like this [1 -1 -0.5 0.5 1 0.7 -0.07 0.8].
I build my demodlator as follows this rule : [1 is actually 1 in binary, -1 is actually 0 in binary] .
my problem is now mapping the output of my demodulator to binary vector of 1/0 accordingly if there's value 1 then it will be mapped to 1 in binary, if there's value of -1 it will be mapped to 0 in binary.
but what about if there's a value ranges between [-1,1] where should I map it in order to get the minimum error of wrong mapped data?
I used according to likelihood theory threshold zero and said if my data that's output of my demulator:
if its values bigger than zero => 1 in binary  , if its values less then zero => 0 in binary ..
but I'm not succeeding to see my data packets correctly, this means Im not succeeding after converting the demodulated data (output of my demodulator that I built )to binary data not getting the same correct data of my packets that I transmitted .. so not enable to see the start/end of my packets !
Any help what should I do in order to convert optimally my data (output of demodulator) to binary data?
Why I need to convert to binary data? because I sent my packet by #SDR software as binary data and my transmitter using GFSK modulation .... so then at the receiver I want to re-map my recieved data to binary data in order to be able to detect my received packet data (when it start, end , detect CRC)
Any help what should I do ? thanks alot.
Attaching for your reference my sampled data by a dangle(it contains the packets as shown in the photo above) it's .dat type:
enter link description here


texts: Detect presence of a specific frequency (or lack thereof)

texts:
I know very little DSP beyond an intro undergrad course. I was looking to write a C algorithm that takes ADC samples and determines if a specific frequency is present (in the range of 100KHz); not caring about any other components of the signal.

I'm not sure where to start. I looked at Goertzel algorithm; would this be a recommendable approach?

I was wondering if it is feasible to analyse a signal at 100KHz with a ARM micro controller, or will I need a special DSP breed?

Should I even be tackling this with a PLL or something similar? If so, does anyone have some design references?

Thanks!


texts: Question about Difference of Gaussian (DoG) algorithm

texts:
I am recently learning about Computer Vision and I am having a trouble understanding Difference of Gaussian (DoG) algorithm. I get how the algorithm works in high level idea, but I am trying to implement my own and I am confused about some steps.
For instance, I am trying to create 5 blur level for each octave, and I am confused about which filter and sigma value applying to which image. Using Matlab, for the first octave, I created a filter and applied:

sigma = 0.5;
gauss = fspecial('gaussian', [5 5], sigma);
blur1 = imfilter(img, gauss, 'replicate');
dog1 = img - blur1;
%Next level
blur2 = imfilter(blur1, gauss, 'replicate');
dog2 = blur1 - blur2;

I am not so sure if this is how I need to apply? Do I apply gaussian filter to previously applied image? I also saw code using k*sigma. I am not sure what k means and how to apply? Oh and what value should I used for sigma? Is it in [0, 1] range or can be bigger than that? Could someone help me on this? Thank you very much.
Thank you so much!


texts: Looking for low-noise audio data to benchmark some noise-processing techniques

texts:
I'm looking for a source of very low-noise audio data for use in benchmarking various methods of random noise correction methods.

My intention is to have a very clean signal to which I can add known random noise such that I can test the efficacy of various methods of dealing with that noise.

It would need to be in some lossless format (preferably .wav).

Does anyone have suggestions?

Thanks!


texts:
I have been reading about "all-digital" or "software-defined radars". Lately, I was trying to simulate an FMCW short range radar in Matlab. One thing that struck me during the simulation is that to achieve a range resolution of around a 100 cm or less, the signal bandwidth is quite large -- in multiple GHz. This means, if one wishes to realise a digital implementation, one needs a really high speed ADC - multiple GigaSamples/second. Such ADCs are either not available and if you find something towards the top end of the high speed ADCs out there, the processing requirements are going to be ridiculously high. And of course, if you are clocking your baseband processor in GigaHertz, the power requirements of the device (typically an FPGA) are again going to be impractical.

My questions:

How, if at all, high-resolution, high-bandwidth digital radars are realised in practice?
what would make sense in case of FMCW radar is to mix the transmitted signal and the target return signal and digitize only the beat frequency for Inverse FFT and further analysis. In that case, could you still call it a digital or software defined radar?


texts: For a signal sampled at $f_s$, the frequency resolution (or "bin width") for an $N$ point FFT is $f_s/N$. Does this mean that the $k$th bin will contain energy from sinusoids within $0.5f_s/N$ on either side of the center frequency of that bin (ignoring spectral leakage for now)? Is this energy averaged? Is it the same as taking a $2N$ point FFT and then averaging every two bins?


texts: When I am in the DT Fourier Domain, and I want to come back to the time domain, which pair do I use? Asking because both pairs have the exact same "form" in the Fourier domain, and that is what is tripping me up. Do I just use pair 2 when I want to go into the Fourier domain and then I have to use pair 1 when I want to come back?? Or is there more to it than that?


texts: Reversing the order of up and downsampling

texts:
When we are dealing with the problem of changing the sampling rate by a noninteger factor like L/M, we first need to the do the upsampling, and then the downsampling.

My question is that under which condition can we reverse this order?


texts: Measure and simulate non-linear distortion

texts:
We can capture frequency response, phase response and reverb in recording of impulse response. Then we can use it to simulate audio system and/or space by convolving it with dry audio signal.

Is there similar technique for measuring and simulating non-linear distortion?

I thought of measuring how the amplitude of output signal corresponds to amplitude of input signal. Then I could simulate this system by waveshaping. Will it work? Is there any better method?


texts: How to get info (mean,std) out of normplot in matlab?

texts:
Lets say we have the following data and then we have a normplot:

load census;
h=normplot(cdate);
[mean std]=normfit(cdate);

I know i can get the mean(1890) and std(62.0484) of the data using normfit fxn, but thats not what i need.

As you can see in the plot, the straight red line is the ideal for the data.

i want to know the mean and std of the line.

I think we can do it from the slope of the line.. using handles of the figure probably with h(3).Xdata and Ydata.., but somehow i cant figure it out.

Can you please help me out.

Thanks


texts: I am aware of matched filter and its application. Now, wondering if there is any application of inverse matched filter? What I mean my inverse matched filter is that convolution of matched filter and the inverse matched filter would lead to close to delta function.


texts:
I am reading about audio processing for my university project and I found out that DCT basicly filters out some frequencies. I wanted to use for example high-pass filter and then DCT for compression, but isn't dct a filter?

What's the difference?


texts:
Suppose $\mathcal{L}$ be invertible system with memory. Does $\mathcal{L}^{-1}$ have memory necessarily?

Intuitively I think the answer is "yes". There are many examples showing that. For instance $\mathcal{L}(x(t)) = x(t-2)$ and $\mathcal{L}(x(t)) = x(\frac t 3)$. Another example which seems problematic to me is
$$\mathcal{L}(x(t)) = \int_{-\infty}^{t}x(\lambda)d\lambda$$The inverse is $$\mathcal{L}^{-1}(x(t)) = \frac{dx(t)}{dt}$$Does differentiator have memory? Of course the main question here is about memory of an invertible system which has memory. Note that here $\mathcal{L}$ can be nonlinear as well.
For clarity, I add some related definitions from Oppenheim's book:
Invertible system: A system is said to be invertible if distinct inputs lead to distinct outputs.
Causal system: A system is causal if the output at any time depends only on values of the input at the present time and in the past.
Memoryless system: A system is said to be memoryless if its output for each value of the independent variable at a given time is dependent only on the input at that same time.


texts: Apply wavelet transform to analyse EEG signal

texts:
I would like to apply The Morlet wavelet transform to analyse my EEG signals. I have many short signals each is only 1 min long. and they all recorded in 30Hz.
I have Two questions:

In the Morlet wavelet,  What is the best scale (alpha) to use in my case ?
About the edge effect: How can I know / compute what portion of my data will be corrupted due to the wavelet "edge effect" ?


texts: Deconvolution Using Response to an Heavy Side

texts:
I'm measuring a "charge" signal in function of time from an amplifier.
Here is a measured signal (x-axis is the time in some arbitrary units, y-axis is the charge in ADU):

I would like to get the "true", original signal, by deconvolving the instrument's response from the measured signal.

Here the response of the instrument to an Heaviside step function:

I guess this should help me characterize the instrument's response, and do the deconvolution.
I read about SVD and Bayesian unfolding, but these don't seem to do the trick (I need to built a transfer matrix first).
Could you suggest me any tool/algorithm to do this?

Thanks in advance for the help.


texts: Convert complex valued sinusoid to real valued sinusoid

texts:
This is the homework problem: convert $x[n]=je^{j\pi n/8}-je^{-j\pi n/8}$ to a real valued sinusoid.
I understand that $\sin\theta=\dfrac{e^{j\theta}-e^{-j\theta}}{2j}$
In the solution, the answers claim that $x[n]=\dfrac{-e^{j\pi n/8}+e^{-j\pi n/8}}{j}$, and I don't understand how to get from the original
$$x[n]=je^{j\pi n/8}-je^{-j\pi n/8}$$
and arrive at
$$x[n]=\dfrac{-e^{j\pi n/8}+e^{-j\pi n/8}}{j}$$
after which it is easy to see  $x[n]=-2\sin(\pi n/8)$


texts: Wrong values calculating FFT with EEG Bands using Numpy

texts:
First of all I have to say I am very new to these matters. I am trying to apply FFT algorithm to some values obtained by EEG bands.
I found the following code but I am unable to make it work correctly.
import numpy as np
fs = 200
#data = np.random.uniform(0, 100, 1024)  # 2 sec of data b/w 0.0-100.0

data = np.array([(-0.27,-0.09,0.16,0.01,0.67,0.65,-1.59,0.11,-0.36,0.69,0.94,0.43,-0.53,-0.57,0.98,-0.97),(-0.74,-0.56,0.76,0.58,1.79,2.22,-4.96,-0.25,-1.19,1.57,2.95,1.47,-1.00,-1.06,2.97,-3.43),
(-0.89,-1.50,1.76,1.72,2.09,4.04,-7.97,-1.30,-2.09,1.11,4.82,2.57,-0.01,0.05,4.64,-6.44),
(-0.84,-2.78,3.00,2.88,2.00,6.36,-11.46,-2.36,-3.22,-0.01,6.96,3.84,1.85,2.11,6.64,-10.31),
(-0.39,-4.27,4.26,4.57,1.09,8.77,-14.78,-4.09,-4.40,-1.87,9.03,5.12,4.86,5.31,8.76,-14.77)])

print("After data")
# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))
print("After fft_vales")
print(fft_vals)
# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)
print("After fft_freq")
print(fft_freq)
# Define EEG bands
eeg_bands = {'Delta': (0, 4),
         'Theta': (4, 8),
         'Alpha': (8, 12),
         'Beta': (12, 30),
         'Gamma': (30, 45)}

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:
   freq_ix = np.where((fft_freq >= eeg_bands[band][0]) &
                   (fft_freq <= eeg_bands[band][1]))[0]
   eeg_band_fft[band] = np.mean(fft_vals[freq_ix])

print(eeg_band_fft)

I am getting the following errors when it goes into the for.
/home/xxxx/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3372: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
/home/xxxx/.local/lib/python3.8/site-packages/numpy/core/_methods.py:170: RuntimeWarning: invalid value encountered in double_scalars
ret = ret.dtype.type(ret / rcount)

What am I missing or what's wrong?
As I said above, I am very new to this matter so please, try to explain everything as simple as possible.
This is the result I get:
band        val
0  Delta  64.642945
1  Theta        NaN
2  Alpha        NaN
3   Beta        NaN
4  Gamma  25.851368

This is how the file I am reading looks like. It is a file  from a 16 electrodes reading with nearly 200k entries.
I have only taken the first entries of the file in order to test this code snippet.


texts:
I am trying generate guided waves on simple metallic plate through piezo-ceramic material, I place two piezo on the surface of the plate at the different location and through Arbitrary pulse generator, I generate the Gaussian pulse of 150KHz and excite the one piezo. This excitation generate guided waves which is sensed by the other piezo. but the problem is, the output signal has very high noises and very low SNR ratio, I hardly even see the signal on the oscilloscope, but whenever I switch to the normal sine wave of the same frequency, and shut off the arbitrary and burst option, then the output signal will have very high resolution and SNR value.

actually I have mechanical background, so I coundnt figure out the source of these noises. As u can see the blue signals, the output one has very low amplitude, where each box has value of 50mV, Moreover, when I remove the probes, these noises still persist there. Can u please guide me what data I should provide to avoid these noises.
Please someone guide me how to reduce signal noises,


texts: How to test I/Q modulation with only one branch (I or Q)?

texts:
Given is a direct-conversion I/Q up- and downconverter system. Receiver and transmitter share the same (10MHz) reference and hence the LO frequency is identical and there is an (unknown) phase difference between them. The transmitted RF signal is given as:

$$
x_{\rm rf} = x_i \cos(\omega_c t) + x_q \sin(\omega_c t)
$$

Now I look at only one channel (e.g. the I channel) at the receiver (that means that any I/Q imbalance of the receiver is irrelevant). At the I branch, this signal is downconverted with an arbitrary phase shift phi and filtered with H:

$$
y_i = H( x_{\rm rf} \cos(\omega_c t + \phi) ) = \cdots \\
\approx H( \cos(\phi) x_i/2 - \sin(\phi) x_q/2 ) \\
= \frac{\cos(\phi)}{2}H(x_i + x_q \tan \phi) \\
= G(x_i + \alpha x_q) \\
= G(x_i) + \alpha G( x_q) \\
= G_1(x_i) + G_2(x_q)
$$

The filter H is contains all filtering effects (mainly anti-aliasing filter).

I want to find a linear relationship between the time-domain samples $y_i[n]$, $x_i[n]$ and $x_q[n]$. A straight forward way is to bring this equation in a vector-/matrix form, model $G_1$ and $G_2$ as FIR filters with $Q$ taps and solve it via Least Squares:

$$
\mathbf{y}_i = \mathbf{X}_i \mathbf{g}_1 + \mathbf{X}_q \mathbf{g}_2
$$

This works nicely and I get the result as expected (in my case, the error between $\mathbf{y}_i$ and the received samples is -45dB). The caveat is I get a set of $2Q$ (real-valued) coefficients but according my derivation, $G_1$ should identical to $G_2$ up to a constant factor. If I would do the same for the $Q$ channel I get again $2Q$ coefficients and a total of $4Q$ real-valued coefficients if I were to model the I/Q relationship that way.

Clearly I can model everything directly in complex domain (with $Q$ complex-valued or $2Q$ real-valued coefficients) but then I cannot look at only the I channel without effect from the Q channel.

Hence a more proper way is to look at this relationship:

$$
\mathbf{y}_i = \mathbf{X}_i \mathbf{g}_1 + \mathbf{X}_q \alpha \mathbf{g}_1
$$

However, this is not a least squares problem any more. For that reason, I sweep $\alpha$ within a reasonable range, compute the LS solution and take the $\alpha$ with the minimum error (this seems to be convex). But if I do this, the error -35 dB (as opposed to -45dB above).

I really do not understand why this approach would not work.
Why? And what is the proper way to model this?

Background info: I am testing my own direct conversion receiver. If I transmit data in only one channel (I or Q) and look at only one output channel (e.g. I) I get the performance that I expect based on linearity-, noise- and phase noise measurements of the receiver. If I look at I/Q at the same time (and model the the effect of the receiver as complex-valued FIR filter) then the error is more than 5 dB worse. For debugging I would like to avoid any interaction (such as I/Q imbalance) between I and Q and only look at the I channel output. But as soon as I transmit a full complex-valued input and look only at the I output my result is again at least 5 dB worse than predicted by noise, linearity, phase noise.

As another step, I use an Arbitrary Waveform Generator (Rohde & Schwarz SMW200A) used as transmitter and a VSA (Rohde & Schwarz FSW) as a receiver where I see the same issues. Both instruments are calibrated and the transmitter has superior image rejection. For these reasons I assume I made a conceptional or algebraic mistake in my derivations/assumptions.


texts: Which IIR form has lowest round-off noise?

texts:
We are making some audio effects using IIR, and seeking ways to optimize its round-off noise. After searched some materials, most of them tells transposed form II has lowest round-off noise, and many IIR libraries are implemented in it. However, some other materials tells direct form I is the best. So which one should I choose?

Here are specifications of our working signal and currently applied strategy:

direct form I;
sample rate in 48000 or 44100 Hz;
input/output data in single-precision float, inner state storage and calculation in double-precision;
running in CPU.

At current, this strategy hears good. But we are still seeking on reduce double-precision calculation to single precision, if possible.


texts: Image registration / alignment using on iOS and Android (using C#/Xamarin)

texts:
I'm trying to align two images of a the same scene that are taken with different cameras - identifying what the scaling factor is that should be applied to one of the images to align them as closely as possible.
I'm assuming I need access to one of the various feature-detection algorithms such as SIFT or BRISK, however I need to do this on both iOS and Android platforms (using Xamarin Forms).
Looking at other similar posts on this forum, https://simpleelastix.github.io seems to do what I want, although is only currently compiled for Windows/Linux and MacOS.
I've also tried to use an OpenCV port, EMGU using a tutorial that covers image matching here, however including the entire OpenCV library in my app will increase the size by ~200MB which is rather a lot for a simple image match.
Does anyone know of an easy to use library, ideally in C#, that can identify the basic image transformations to match two images?
Grateful for any help,
Oliver.


texts:
I have been studying about the sampling theorem and it seems that even though we sample at a frequency with the nyquist criterion, the harmonics (due to sampling process) remain within the nyquist frequency range, as illustrated below:

If we sample a 5MHz sine wave with 20 MHz sampling frequency (nyquist frequency being 10 MHz), its fundamental and harmonic frequencies (up to 4th) lie at:

fin=5 MHz (fundamental)

1st Harmonics: (around fs)
fs+fin=25 MHz
fs-fin=15 MHz

2nd Harmonics: (outside nyquist)
fs+(2*fin)=30 MHz
fs-(2*fin)=10 MHz

3rd Harmonics:
fs+(3*fin)=35 MHz
fs-(3*fin)=5 MHz (Inside nyquist)

4th Harmonics:
fs+(4*fin)=40 MHz
fs-(4*fin)=0 (DC)

I understand that when doing a high level simulation in matlab, we do windowing or coherent sampling to overcome this harmonics falling into the nyquist range or signal band by choosing a prime number of cycles.

But practically, if we give an arbitrary frequency to an ADC, which usually has a sample/hold circuitry in the front, aren't these harmonics generated by it?

Or do we have the control of the harmonics' power so that we can adjust it to a minimum (depending on the application).

(This question also relates to my previous question - Higher order harmonics during sampling)


texts: Duration of downsampled signal in analog frequency axis

texts:
Suppose, I sampled a signal, $x(t)$ and let $x[n]$ be the sampled signal of $x(t)$. Let then I downsampled $x[n]$ by a factor of $M=3$. Let $x[n]$ after being downsampled by the factor of $M=3$ yields $y[n]$.
Now, my question is: if I plot the frequency spectrum of $y[n]$ along y-axis with frequency (in Hertz) along x-axis, then, will I get only $M=3$ copies of frequency scaled version of the original spectrum(spectrum of the original signal, $x(t)$)  or I will get infinite copies of the original spectrum ?


texts: Transient response of system with single pole $0 \le p < 1$

texts:
$G(z) = \frac{1-p}{z-p}$

  If the value of p satisfies $ 0 \leq p < 1$ there are no oscillations
  in the transient response.

Question: Why is that $\uparrow$  true? I know roughly what a transient response is but how is the relation between the position of the pole and the nature of the transient response?

Edit: In another paper by the same authors:

  Enforcing the stability of the controlled system means ensuring that
  the pole p is non-negative and less than 1.

Now I'm completely out, stability does not depend upon the position of the pole as long as they are inside the unit circle of the z-plane I have thought?

Paper 2:

http://www.iste.uni-stuttgart.de/fileadmin/user_upload/iste/zss/publications/supplementaryMaterial/2014-ICSE-Control-PaperReviewCopy01.pdf

Look at the sentence short above eq. (10) on page 5

Paper 1 is not available on the internet as fas as I know. It's Antonio Filieri et al. "Software Engineering Meets Control Theory"


texts: Dimensional reduction from DWT with threshold

texts:
I have been trying to find out how can the discrete wavelet transform (DWT) be possible to reduce dimension of data.

Then I saw the question which is seemingly related to my work:

Feature extraction/reduction using DWT

But after seeing the post, I have a question in my mind.

Is it possible to say that setting DWT coefficient values to zero which is lower than threshold is the dimensional reduction?

I mean, suppose we have 5 dimensional vector
<1,4,3,5,2> which can be thought as result of DWT.
And after setting some  vector's  values to zero in case the threshold is 3, we get <0,4,3,5,0>.

I think it seems inappropriate to say that this is the dimensional reduction. Because I understand the dimensional reduction is the actually reduction such <1,4,3,5,2> -> <4,3,5>

So I searched some paper to give some light on my curiousness. Some researches say that they select DWT coefficients for dimensional reduction.

But I wander how it can be possible. Because,,, what if we want to reconstruct original signal or data from DWT coefficients?

I guess ,by selecting some values from DWT coefficients (<1,4,3,5,2> -> <4,3,5>), they already lose the meaning of original data. It should be <0,4,3,5,0> not <4,3,5>. So I think they should post some encoding-decoding technique to replace selected DWT coefficients (<4,3,5> -> <0,4,3,5,0>), but they don't.

Am I misunderstanding?


texts: How to make sense of audio file's $\textrm{dBFS}$?

texts:
I'm trying to build a song visualizer using the FFT method. When I analyze the audio, I get the $\textrm{dBFS}$ value for each frequency range.

From my understanding,

$\textrm{dBFS}$ is in the range $[-\infty, 0]$, where $0$ is the most loud.
The formula for percentage loudness is $10^{\rm dBFS/20}$.

However, the $\textrm{dBFS}$ values I'm getting are consistently below $-20$, i.e., $10\%$ loudness, even at loud parts of the song. Why is this so? Am I calculating the loudness wrongly?

Also, how can I calculate the overall loudness at a point in time? Is it the sum of all the values, or the maximum?


texts: Analog Filter Prototype to Direct Form Second Order Cascade

texts:
I've been using the Faust programming language a lot lately for experimenting with DSP, and I've been digging into their implementation of an Elliptic (Cauer) Lowpass filter. Here's an example of such a filter with order three.

If you're unfamiliar with Faust, hopefully the example there is still fairly clear: they've designed a 3rd order Cauer filter with the properties listed in the comments and have implemented it as a second order direct form filter feeding into a first order direct form filter with the coefficients as listed.

I understand the analog prototype they've designed with [z,p,g] = ncauer(Rp,Rs,3); in Matlab/Octave, and I think I have a good understanding of how to factor the transfer function into cascaded first- and second-order filters. What I don't understand is their use of poly in Octave to find the coefficients there, and why the frequency of this elliptic filter is only governed by the last coefficient in the final first order filter? Intuitively it seems like so many more of those coefficients should depend on the frequency of the filter, especially if we have to consider multiple sample rates. Did they skip some steps here or make some assumptions that I'm missing? Any explanation to my confusion here would be greatly appreciated, thank you!

Update: I find it interesting that in Matlab/Octave I can compute [z,p,g] = ncauer(0.2,60,3); sos=zp2sos(z, p, g), deriving coefficients for second order sections without ever specifying the cutoff frequency. Maybe that's where I'm most confused?


texts:
I have to fix the amplitude of measured IQ receiver data. According to my information about the data, I can estimate that one signal needs to be corrected to the other in amplitude. The phase offset/imbalance can be ignored in this case. I would like to implement this in MATLAB. At the moment my I/Q signal looks like this:

I would like to correct it, that is shows a circle (would be the ideal result). But until now, I found no direct way to do that. I have no possibility to train or measure the amplitude imbalance from the I/Q receiver direclty. I have only the data. What have I to do, to normalize the Q component that the chart shows a circle?

I have already tried to calculate the amplitude offset: o=sqrt(Q^2/(1-I^2)), but this doesn't work. Both signals have a small DC offset (really small, about $10^{-4}$ in relation to values between $0$ and $1$), but this is not relevant.


texts:
As known, Bandwidth is the difference between the upper and lower frequencies in a continuous band of frequencies.
OK, what means when saying that bandwidth range is 15 ~ 20 KHz ?  Does that mean that the bandwidth should be 5KHz ?


texts: Question regarding using cardinality in pattern recognition

texts:
So I've read this article for a project im working on:

Ortiz-Cataln: Cardinality as a highly descriptive
feature in myoelectric pattern
recognition for decoding motor
volition, Chalmers Uni, published in "Frontiers in Neuroscience", 2015

They use cardinality in pattern recognition and claim that it is "highly descriptive".

I can't understand how they have used it or why it would be of any use, if my understanding of cardinality as simply a measure of the size of dataset is correct. Can anyone explain to me what they mean and how they have used cardinality?


texts: Impulse response convolution and normalization2

texts:
when I take inverse Laplace transform of a system transfer function \

Lets say LPF whose TF is

$$\frac{Y(s)}{X(s)} \triangleq H(s) = \frac{W}{s+W} $$

the inverse Laplace/impulse response is

$$h(t) = We^{-Wt}u(t) $$

where $u(t)$ is the Heaviside unit step function:

$$ u(t) \triangleq \begin{cases}
1 \qquad & t \ge 0 \\
0 \qquad & t < 0 \\
\end{cases} $$

Now to see the system response to a square wave $x(t)$ with $|x(t)| = 1$, I need to convolve

$$ y(t) = h(t) \star x(t) $$

Now if you look at $h(t)$, the maximum amplitude of $h(t)$ is

$$\max{|h(t)|} = W $$

Then $y(t)$ is amplified by $W$ times.

So what is happening? how do I normalize this?

Should I normalize this by $\max{|h(t)|}$ or
Should I normalize this with

$$W_{z_1} W_{z_2} \cdots W_{z_n}/(W_{p_1} W_{p_2} \cdots W_{p_n})$$

(product of zeroes)/(product of poles) of transfer function?

Why is this even happening?


texts: How to learn image processing / computer vision in a practical way?

texts:
I am new to the field of image processing / computer vision and I find it very hard to get my hands on the subject.
I read a lot about the basics/theory but I feel in no way able to use this knowledge in a practical way.

Do you have any tips on how to translate my theoretical knowledge into practical knowledge?
Are there any "cookbooks" for image processing that are not bound to a software/library? (I plan to use Halcon but only find resources for Photoshop or OpenCV.)

Thanks in advance!


texts:
I have an image of bricks of different sizes. What I want to do is when I click on the corner of the brick it should mark the lines around the brick. Can anyone please help me how can I do that...
here is the image


texts:
I am trying to understand Haar transform of simple vector of numbers.

Some references I found say, that first level transform of $ X = (4,6,8,10,13,9,3,3)$ will be $\sqrt{2}(5,9,11,3, -1,-1,2,0)$. Others say, it will be $(5,9,11,3,-1,-1,2,0)$.

Which one is true? Or are they both correct in some way and am I missing something?

And important question, which one do I use when I am performing Haar transform on image and this vector represents one row/column of pixels in image?


texts: Which kind of machine learning algorithm is Hidden Markov Model?

texts:
I have used Hidden Markov Model algorithm for automated speech recognition in a signal processing class. Now going through Machine learning literature i see that  algorithms are classified as "Classification" , "Clustering" or "Regression".
Which bucket does HMM fall into? I did not come across hidden markov models listed in the literature.


texts: Finding Signals when the baseline varies greatly by signal set

texts:
I have milliwatt consumption data from home appliances.  I'm collecting this data in hundreds of households.  Based on the machine model, the baseline value for a given household can vary from 0 to 500 milliwatts.  The baseline non-zero values vary by 5 to 15 milliwatts and the data is reported every 13 seconds.  The vast majority of data is baseline, and there is only signal when people are running the appliance.  Perhaps 10 hours a week is real data.  The signals are very pronounced, going from baseline to values exceeding 100,000.  I have 1/3 of a billion records (growing by several million records a day) so signal detection must be automated.

I really just need, start and end time of the event.  We wrote python code that looked for the jump in signal value and return to baseline.  The trouble is that each machine has a different baseline value.

Are there simple tools/techniques that can determine from recent context the end of a pronounced signal?  My preference is python but I'm happy to use anything that works.  Thanks in advance for any feedback.


texts: In general, the output length of an FIR filter is given by (N+M-1) samples, where N stands for signal length and M stands for the number of filter coefficients. Can the same formulae be applied for IIR Filters? Let's say I have 4 input samples, 3 feed-forward coefficients, and 2 feedback coefficients. Then can I tell that my output length is 4+(3+2)-1? It is also known that the impulse response of an IIR Filter is not of finite duration. Does this thing affect the output length of an IIR filter, in any way?


texts: Can a Kalman Filter Be Applied Using Measurement-Space Dependent Sensors?

texts:
I'm currently attempting to apply a Kalman filter to track the angular position, velocity, and acceleration of a bike wheel, and I'm having a lot of trouble, so I want to check if I'm even applying the right tool for the job.  As you will see, this is not a "classical" Kalman implementation, and I want to check if I can still use it.

I have a wheel with n reed switches evenly spaced around the circumference of the wheel.  The way that I have the filter set up is that every time one of the reed switches is triggered by a magnet sitting on the bike's fork, a new measurement is fed to the filter that has the position and velocity (z = [x_prev + wheel_circumference/n ; wheel_circumference/(n*dt)] where x_prev is the previous position, wheel_circumference is the circumference of the bike wheel, and dt is the time since the last measurement).  (I have tried this system using only position measurements, and I am still having trouble.)

Unlike the classic application of the Kalman filter, the measurements in this system occur periodically in position, not in time, i.e. the sensor(s) sends a measurement to the filter iff the wheel's position reaches a certain position.  Compared to a traditional sensor, this implies a heavy correlation between position, velocity, and inter-measurement period.  Does this correlation negatively affect the Kalman filter?


texts:
I want to implement a spatio-temporal filter for my images which are a set of silhouettes of a cycle of walking. So I gather all the images in an array and convolve the 3D array with my spatio-temporal kernel, I mentioned below.

but my result is'nt correct. Here is my code please tell me what is my mistake.

import cv2,numpy as np
import matplotlib.pyplot as plt
import glob
from PIL import Image
from scipy import ndimage

def gkern(l,sig):
    """\
    creates gaussian kernel with side length l and a sigma of sig
    """

    ax = np.linspace(-(l - 1) / 2., (l - 1) / 2., l)
    xx, yy = np.meshgrid(ax, ax)

    kernel = np.exp(-0.5 * (np.square(xx) +  np.square(yy)) / np.square(sig))

    return kernel / np.sum(kernel)

t = np.linspace(0,20,23)

h1 = np.array(gkern(5,5))[None,:,:]
h2= np.array(np.sinc(t-2))[:,None,None]

img_array=[]

for filename in glob.glob('t/*.bmp'):
    img = np.array(Image.open(filename).convert('I;16'))
    img_array.append(img)

v = np.array(img_array)

def fft(array):
  fft = np.fft.ifftshift(np.fft.fftn(np.fft.fftshift(array)))
  return fft

def ifft(array):
  ifft = np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(array)))
  return ifft

def conv_3D(array, kernel):
  conv = np.abs(ifft(fft(array)*fft(kernel)))
  return conv

st_kernel = conv_3D(h1, h2)
f_kernel = ndimage.convolve(v,st_kernel)

plt.figure(0)
plt.imshow(f[0,:,:],cmap='gray')


texts: need general resource about image matching

texts:
i'm working on image matching mainly SIFT and SURF for my thesis. i read Moravec's and harris corner detectors and understand them quiet well in my own opinion.A COMBINED CORNER AND EDGE DETECTOR[1]

I even understoond harris-lapace detector but had no luck with understanding it's descriptor.Indexing based on scale invariant interest points[2]

i understood SIFT and SURF too.

what i didn't really understand is some of these topics : harris-laplace descriptor , how does hessian matrix help to detect good interest points , RANSAC(have general idea but need details for image matching) , Haar wavlet( also have a general idea, don't know how it's output look like and how i can work with that)

that's pretty much it.
i would really appreciate If anyone can help me with finding general reference for these topics.

(i have sent this question on stackoverflow. they suggest here would be better place to ask it and get proper answer)


texts: Are pulse shaping and windowing functions the same in concept?

texts:
I read a digital communication problem and they interchange terms as windowing function and pulse shaping and refer to it as same entity. For example the following windowing function with Transition time $T_R$ used to smooth transitions between symbols is also refered to as pulse shaping

$w_T(t)= \left\{ \begin{array}{ll}
         \sin^2(\frac{\pi}{2}(0.5+\frac{t}{T_{TR}})) & \mbox{if $-T_{TR}/2 < t<T_{RT}/2$};\\
        1 & \mbox{if $T_{TR}/2 < t<T-T_{TR}/2$}\\
 \sin^2(\frac{\pi}{2}(0.5-\frac{t-T}{T_{TR}})) & \mbox{if $T-T_{TR}/2 < t<T+T_{TR}/2$};\\
\end{array} \right. $

My understanding is that a pulse shaping function limits signal bandwidth to transmission bandwidth.

How come the two are refered to as same?

Thanks


texts:
I am doing fmcw radar simulation. I have an output of mixer and from that I want to reconstruct the amplitude and phase of the signal. Therefore, I am taking hilbert transform and by the abs() and angle(), I can get the phase and amplitude information. Since the hilbert transform is in time domain, I want to map this time domain axis into frequency domain so I can see the value of amplitude and phase at different frequencies. FFT is not an option since by taking the fft, I will get the beat frequency.Therefore, how to map this time domain axis to frequency domain.

here are some of the relations between time and frequency. t = f/K, where

K = Bw / Tm
f = frequency,
Bw = bandwidth,
Tm = transmit time

I guess here, deltaf = 1/ t(2)-t(1) and fmax = 1/deltat, I am trying to map the time axis into frequency, but it gives dimension mismatch. Still there is something missing( some multiplyting factor or diving).Acoordig to my SV, its a demultiplexing process. That woud be really good if anyone can help me.

mixed1=(y.*(Y2));
hil = hilbert(y_filt);
figure, plot(t,abs(hil));
figure, plot(t,angle(hil));
deltaf = 1/t(2) - t(1);
fmax = 1/delta t;
f =?            %%% how to define frequency axis so i can plot amplitude & phase
figure, plot(f,abs(hil));
figure, plot(f,angle(hil));


texts: Why is the signum function $=2u(t) - 1$?

texts: Basically I couldnt find proof for this anywhere though it's a very simple and basic equation.


texts: I have a sequence which is made of complex exponentials with arbitrary amplitudes and initial phases.  No noise.  That is it is a sum of the complex exponentials of the form $A_i e^{jnw_i+p_i}$.  For example, let us say i ranges from 1 to 10.  What is the minimum necessary data (number of samples) to calculate their amplitudes, initial phases and frequencies exactly i.e. $A$, $w$ and $p$ s?  And how can we find them? $j$ here is square root of -1 and $n$ is the sample number.


texts: Filter Coefficients in MATLAB are always float data type

texts:
I will be using the filter coefficients generated by MATLAB in my code that will run on a microcontroller. The poor microcontroller is very bad in floating point Arithmetic. So What do I do? Shall I directly truncate all float values to nearest int? shall I do ceil? or floor? Is there any function that can automatically do this conversion for me?

Update
All my filter coefficeints are coming like: 0.465, 0.76, 0.23 etc, so If I use ciel then all will become zero. What should I do now?

Actual Filter Coefficients are:

> -0.00385638
>  0.004944457
>  0.01505063
>  0.018768283
>  0.009635631
> -0.01192891
> -0.035809426
> -0.045043857
> -0.023993426
>  0.032106934
>  0.111777547
>  0.18997355
>  0.238374966
>  0.238374966
>  0.18997355
>  0.111777547
>  0.032106934
> -0.023993426
> -0.045043857
> -0.035809426
> -0.01192891
> 0.009635631
> 0.018768283
> 0.01505063
> 0.004944457
> -0.00385638


texts:
Let $$D(s) + KN(s) = 0 \tag{1}$$where $D(s)$ and $N(s)$ are polynomials of $s \in \mathbb{C}$ such that $\text{Deg}(D) = n, \ \text{Deg}(N) = m$ and $n\ge m$. The root locus method tells us how the solutions of $(1)$ changes as we change parameter $K$ from $K=0$ to $K = \infty$.
I'm trying to understand these extreme cases. Let $K \to 0$ and we have $$D(s) = 0$$ so in this case, the set of the solutions is $A = \{s \in \mathbb{R} | D(s) = 0\}$.
Now let $K \to \infty$, if we choose $s$ such that $N(s) \not = 0$ then the answer will be infinity. So we should choose $s$ such that $N(s) = 0$. In that case, if $(1)$ holds, we also have $D(s) =0$ which means $N(s)$ and $D(s)$ have the same factor but this isn't the result that should be obtained. Curiously, if we rewrite $(1)$ $$\frac{N(s)}{D(s)} = -\frac{1}{K} \tag{2}$$ and let $K\to \infty $, one possible case that $(2)$ holds is $N(s) = 0$ which gives us $m$ solutions and this doesn't require $D(s)$ have the same factor as $N(s)$! Why this happens? And why the first solution is wrong?
Example: Let $D(s) = s^2 - 4$ and $N(s) = s + 3$. So $(1)$ becomes $$s^2 - 4 + K(s+3) = 0$$If $K \to \infty $ and $s = -3$ then $9 - 4 = 0$ which is clearly wrong. On the other hand, rewriting the equation $$\frac{s+3}{s^2 - 4} = -\frac{1}{K}$$ If $K \to \infty $ and $s = -3$ then $\frac{0}{9-4} = 0$ which is true, of course.


texts: Confusion regarding pdf of circularly symmetric complex gaussian rv

texts:
Considering a random variable $x$ that takes in values from a complex domain. Its real and imaginary components are totally uncorrelated. I am following this link and
also studying this document.

In the wikipedia link, the pdf for a single observation is given and $k =2$ for bivariate gaussian as my assumption is that the real and imaginary are totally uncorrelated and are gaussian respectively. There is a 2 in the denominator in the pdf but the log likelihood for the complex case does not have any 2 in the denominator. The log-likelihood for $k=2$ would be $$-N \ln| \Sigma| - (x_n) {|\Sigma|}^{-1}(x_n)^H - k N \ln \pi$$
where $\mu=0$ as I am assuming zero mean r.v $x$

Confusion 1:

I was thinking that the variance of $x$ was
$$\Sigma = \begin{bmatrix}\sigma_1^2 & 0\\ 0& \sigma_2^2\end{bmatrix}\,\text.$$
So for $N$ samples (observation), the joint pdf would turn out to be
$$P_x(x_1,x_2,...,x_N) = \prod_{n=1}^N\frac{1}{\pi \sigma^2_x} \exp \bigg(\frac{-{({x_n})}^H ({x_n})}{\sigma^2_x} \bigg)$$
where $\sigma^2_x = [\sigma_1^2,\sigma_2^2]$. Would there be a 2 in the denominator of the pdf and what is the correct pdf expression?
Confusion 2:

In the document, the expression for the pdf for complex case looks different from the wikipedia link. Are they the same but maybe I am missing some link between them? Which pdf should I use?
Confusion 3:

When simulating the r.v in Matlab with zero mean and variance 1, I am halving the variance because the real and imaginary part's variance should add up to 1. Then, would the pdf also contain half of the variance as $\sigma^2_x/2$?


texts: Different representations of images and their purpose

texts:
I'm doing a small overview of different image representations used in image processing and computer vision. For every representation, I would like to have (at least) the following information:

the name of the family of representations
some basic characteristics of the representation family (what information is retained in the representation)
what tasks are they typically used for
what are different concrete types of representation from this family
references for different representations and references for their usages

In addition to the representations I already know / have information about it (as an answer listed below), I am mostly interested in compressed domain (e.g. Fourier domain, wavelets), but knowing about other representations, e.g. edge-based, or whatever else I don't know about yet would be great.

If at all possible, I would appreciate if you could provide references (e.g. the article that proposed the representation, or a book explaining the representation), as indicated with the reference-request tag.

Disclamer

I know list-questions are usually discouraged, but we had something similar about segmentation and thresholding, and it worked well.

If the info I already provided is more suited as an answer, I could edit it out and post it as an answer. Also, I'm thinking maybe this would be a good question for community wiki.


texts: Why the delta at the end of the approximation?

texts: The equation to approximate an input signal with a unit impulse in Continuous Time, is shown below, before we take the limit $\hat{x}(t)=\frac{lim}{\Delta\rightarrow0}\sum^{\infty}_{-\infty}x(k\Delta)\delta_\Delta(t-k\Delta)\Delta$ <-- why is there a final $\Delta$ multiplying the $\delta_\Delta(t-k\Delta)$?


texts:
When deriving DFT from DTFT,we sample the frequency domain with uniformly spaced samples,hence adding periodicity to time domain.

But DFT requires a limited support: we take only 1 period.

Does that make circular convolution the same as a linear convolution with the signal's periodic extension?


texts: Detecting a specific pop in a real time audio signal

texts:
I am listening for a very definite pop in a real time audio signal. So far I have managed to get the audio signal in at a sampling rate of 44100.0 Hz and 2048 frames per second. I have visualized the waveform and computed and visualized the FFT. I am currently looking for distinct features to recognize this pop sound (it is the sound of a ping pong ball on the table).

I found a few papers which suggested features such as the zero crossing rate of a percussive sound like this http://www.csl.sony.fr/downloads/papers/2000/gouyon-dafx2000.pdf. However I am struggling a bit to identify the percussive envelop within the signal. I would like to identify a peak in a certain frequency range within my FFT as well, then combine the existence of an envelope with a peak in the correct range to determine whether a sound was heard.

I have two questions, how would I identify an envelope in the time domain, or how would I approach this problem. Further, how could I identify a peak in the frequency domain when I have a very fine mesh FFT (lots of bins and noise). I also have a lot of pink noise as this is a real time audio recording, so I have louder frequencies toward the bottom end than the top end, but I haven't managed to get those frequencies out of my FFT (I don't need to worry about converting back to the time domain).

Identify envelopes beside natural noise
Identify FFT peaks beside natural noise
Any experience with noticing percussive transients?

Thanks


texts: How to find $h[n]$ system response of this equation?

texts:
$2y[n-2]-2y[n-3]-4y[n-4]=x[n]-10x[n-1]-4x[n-2] + 4x[n-3]$  is the system that I'm looking for the response.

I transformed this system via using Z transform:

$$\frac{Y(z)}{X(z)}=H(z)=\frac{z^4 - 10z^3 - 4z^2 + 4z}{2( z^2 -z -2 )}$$

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

How to find ROC to inverse the Z transfom to $h[n]$ ?


texts: Equation Derivation: Help with intermediate steps

texts:
As per the attached image above, I don't know how the equation in the first line transformed to the form in second line.

I tried solving it by multiplication but the result was different.


texts: EbNo, SNR, and BER simulation for BPSK with convolutional code

texts:
Im currently simulating BER performance of CCSDS convolutional code. Im having some troubles in specifying EbNo and SNR in gnuradio blocks.
From matlab webpages on BER simulation for RRC pulse and Convolutional codes, the SNR was defined by the equation

$SNR|_{dB} = EbNo|_{dB} + 10log(k\times rate) - 10log(sps)$

For uncoded BPSK, rate = 1 and k = 1 bit/symbol and sps is the number of samples per symbol.
Noise variance was calculated as

$\sigma ^{2} = antilog(\frac{-SNR|_{dB}}{10})$

I'm trying to import these equations into gnuradio, starting with uncoded BPSK shown in the first flowgraph. The bpskPulseShapeRRC block is
a hierarchical

map[-1,1] -> char to float -> interp_fir_filter(int(sps),taps).

$\sigma$ was used as the noise voltage of the channel model
block. Running the flowgraph as is resulted in a BER curve that was quite off. Surprisingly, using a noise voltage of

$\frac{\sigma}{sps}$

(or setting the bpskPulseShapeRRC filter gain to sps) produced correct BER curve (I don't understand why. Any explanation is welcomed).

Coming back to CC FEC, the second flowgraph was executed after setting the coding rate to 0.5. Unfortunately, the BER values obtained were
far off as compared to BER curves from the CCSDS Green Books. My guess is, this has something to do with incorrect noise voltage calculation. Any suggestions?

Flowgraph for uncoded BPSK

Flowgraph for BPSK with CC FEC


texts:
Our current project requires us to do some analysis using Wavelet Transform. Can anybody suggest me a practical book, preferably with MATLAB or C examples.
I am currently reading some tutorials, but it is not giving me a feeling as I have for Fourier transform. I need a book having many practical examples with source code.

Really appreciate for your suggestions.


texts: Compute hilbert transform from DCT- FFTW3 REDFT00

texts:
My problem is similar to this but not the same.

Edit: Here is a description of Neumann boundary conditions for differential equations:
Neumann-type boundary conditions means derivatives are specified at the end of the domain. In this case I refer to the specific case where the derivative vanishes at both the ends of the domain. Henceforth, such boundary condition is referred to as Neumann type.

I have 1D wavenumber domain equivalent to a fourier space with integer multiples of a fundamental mode. At the end of the domain neumann boundary conditions are satisfied.

As such, my spatial real signal has basis functions of cosine types.  Due to neumann boundary conditions, to project the signal onto the fourier space, dct type 1 transform is used. Specifically, to accomplish this I use the REDFT00 transform of FFTW3 library.

Now, since it is DCT type 1, only half the length (plus 1) of spatial signal is stored in a 1D array with n bins. Subsequently, I compute the DCT type 1 store it in another array with another n bins. After that I'm confused how to compute the complex part of the hilbert analytic signal from DCT-type1.

It would be better (but not mandatory) to share a piece of code used with FFTW3 REDFT00.