texts: How to simulate measured data with antialiasing filter

texts:
I wish to simulate measured data for developing signal processing methods. All properly measured data will have been through an antialiasing filter. How do I generate such simulated data? I have looked on this site and not found a simple answer (perhaps I have missed something).

For example the naive and wrong way to do it is a follows. Suppose I wish to simulate an exponential decay (my actual equations are nonlinear and much more complicated) using this equation

$y(t)  = Exp(-4.t) (0.3  sin(775 t)+cos(775 t))$

If I sample at 1500 samples per second for 4 seconds I get a time history (only part shown) of

I can take the DFT of this and get a spectrum. The theoretical spectrum, based on working out the Fourier transform of the above equation is

$H(f)=(5.99061 + i 0.159155 f)/(15214.4 + i 1.27324  f -
 f^2)$

where f is frequency in Hz.

Comparing the modulus and phase gives

Clearly aliasing has created the difference. (Although the theoretical spectrum will be a little different due to the presence of the antialiasing filter.)

What options do I have for generating simulated time histories that look like they have been through an antialiasing filter?

My thought so far is to interpolate the data. Resample at a very high frequency, and pass through a digital antialiasing filter. However, the resampling at a high frequency will result in aliasing at that high frequency so I am not technically correct. However, if my time history has very little high frequency then this may be good enough. Am I on the correct lines?


texts: Blind equalization with constant modulus algorithm in GNU Radio

texts:
I'm currently simulating the effect of multitap channel on BPSK as well as its remedy i.e. equalization.

The modem is expected to work with satellite TT&C transponder.
The communication channel is fixed line-of-sight and can be safely regarded as non-fading. Typically, satellite TT&C doesn't employ the use of a pilot signal.

My understanding of equalization is pretty rudimentary and so I have decided to stick to equalization blocks that are already available in GNU Radio (i.e. the constant modulus algorithm (CMA) equalizer block) before I write one by myself. From my understanding, this is a blind equalization algorithm, which I think is appropriate since there is no pilot signal involved.

Now, this block expects the number of taps, the modulus of the signal, an equalizer gain coefficient and a number of samples per symbol. In the simulation below, I'm using a channel with 8 taps $$\begin{pmatrix}1.0\\ 0.25-0.25j \\  0.50 + 0.10j\\ -0.3 + 0.2j \\ 1.0\\  0.25-0.25j\\  0.50 + 0.10j\\ -0.3 + 0.2j\end{pmatrix}\text.$$

I have configured the CMA to have 9 taps. I have used different gain taps in the range 0.001 to 0.05. Smaller gain values result into a very non-flat response at the output of the Costas loop(Ideally, the input signal to the decision device i.e output of Costa's loop should be flat). My questions are:

why is the BER performance so poor?
how to configure this block properly?

Additional comments

The number of taps: A satellite transponder aboard the satellite is connected to a number of RF components (eg. amplifiers, omux, imux, mixers just to name a few) which are specific to each satellite. In addition to that, a satellite ground station has a number of RF components as well. My point is there is no way of knowing for sure the number of taps in advance. How can blind equalization algorithms cope with this?


texts:
if we take a sinus with Fd=100KHz with fs=1000KHz sampling rate We do an FFT with 256 members we get energy leakage.
Why choosing Fd=113.28Khz solves this problem?
Thanks


texts: Voice Biometrics (aka Speaker Recognition) - Standard accuracy rate? State of the art tech?

texts:
I'm working on a voice biometrics project and am looking for some guidance on two key issues:

First, does anyone know if a standard exists for voice biometric system accuracy? I have been unable to find any specific numbers, but it looks like most publications aim to achieve an accuracy of at least 80%.

Second, what is considered the most "state-of-the-art" system? From what I've gathered, current work on feature extraction has largely ranked on MFCC's over LPC's. Acoustic modeling seems split between GMM's, HMM's, Vector Quantization, etc. How does Dynamic Time Warping strengthen accuracy?

If anyone is knowledgable, I would greatly appreciate some help.


texts:
I am having trouble understanding the relationship between a frequency function and it's Inverse Fourier transform.

The Frequency function is

$$\frac{1+0.8(e^{-j 2\pi f}+e^{j 2\pi f})+0.64}{1+1.4\cos(2\pi f)+0.49}$$

And the inverse Forier transform is
$$\frac{(-0.7)^{|n|}}{1-0.49}+\frac{0.64(-0.7)^{|n|}}{1-0.49}+\frac{0.8(-0.7)^{|n-1|}}{1-0.49}+\frac{0.8(-0.7)^{|n+1|}}{1-0.49}$$


texts: What is a speaker behaviour for out of range frequencies?

texts: I was watching this video. The experiment is to create a square signal and record what is outputted by a speaker. We see two diracs signal (+ and -) that looks like the derivative of the signal. Can anyone explain this phenomena? My intuition is that it is caused by some really high frequencies that are out of range for the speakers designed for human. But why the derivative?


texts: Error in measuring wavelength/frequency using autocorrelation

texts:
I am currently using these two fonctions to compute to autocorrelation of a signal in order to measure the wavelength:

import numpy as np
def Cross_corr(x,y, mode = 'normalized'):
    if mode == 'normalized':
        x_norm = (x - x.mean())/(x.std()*np.sqrt(x.size)) ##### on normalise et centre le signal x
        y_norm = (y - y.mean())/(y.std()*np.sqrt(y.size)) ##### idem pour y
    else:
        x_norm = x - x.mean()
        y_norm = y - y.mean()
    return np.correlate(x_norm,y_norm, 'full')

def autocorr(x, mode = 'normalized'):
    return Cross_corr(x,x, mode = mode)[x.size-1:]

After the global maxima, this gives me a nice first peak form which I can find the dominant wavelength of my signal.

The question is: How do you compute an estimation of the error of the wavelength obtained by autocorrelation ? (i.e. my original signal is not a perfect sinusoid)


texts: How to reconstruct 2D flat image from a series of pictures of a tube

texts: I have pictures of an oesophagus (ie a tube) taken under standard white light looking down the tube. I have about 10 images taken over 20cm in series. The image looks forward down into the tube and examples can be found here. Each image is taken serially down the tube at the same session and the centre is positioned largely the same in all the pictures. I would like to be able to combine these images into one (I guess 3D cylinder) in order to then create a flat 2D surface image of the surface lining the tube. Is this possible?


texts: Confusion in the spectrum of FSK modulated signal by using matlab fskmod module

texts:
I want to do a comparison of power spectral density of BFSK, 4FSK and 8FSK by using matlab baseband module "fskmod".Following is the code:

M = 2;
freqsep = 32;
nsamp = 4;
Fs = 64;
x = randint(1000,1,M); % Random signal
y = fskmod(x,M,freqsep,nsamp,Fs); % Modulate.
ly = length(y);
% Create an FFT plot.
freq = [-Fs/2 : Fs/ly : Fs/2 - Fs/ly];
Syy = 10*log10(fftshift(abs(fft(y))));
plot(freq,Syy)

With this configuration, I get a very good spectrum but when I change nsamp to 8 or any other value,multiple lobes start to appear between two frequencies. The confusion is in the understanding of the relation of nsamp and frequency separation. Could any one solve this confusion???


texts:
I am learning signal processing for my subject and have some problems with the unit-power information signal. I refer to this paper: https://arxiv.org/pdf/1906.03949.pdf
In this paper, they formulate the received signal as: $y=h_{sd} \sqrt{p}s+n$ with $h_{sd}$ is the channel, $p$ is the transmit power, $s$ is the unit-power information signal, and $n$ is the noise.
What is the unit-power information signal? As I understand, the unit-power signal is the signal that has the power of 1. However, how can I formulate it (e.g., using Matlab)? For example, with BPSK signal, the signal is formulated as:
$s_1(t) = A_c \cos(2\pi f_ct)$ and $s_0(t)=A_c \cos(2\pi f_ct + \pi)$ for bit 1 and 0, respectively.
With a given transmit power $p$, how can I turn $s_1(t)$ and $s_0(t)$ into unit-power signal?
Is there any suggestion or reference?


texts: MATLAB appears to be applying my filter to the spectrum twice, once on each side

texts:
So it looks like MATLAB is applying the filter to the "negative" and "positive" frequencies separately and I end up with two band pass regions with a null around zero. I want the whole spectrum filtered once.

[b,a] = butter(2,[.18 .3739],'bandpass');
hm = filter(b,a,Y);
figure;
plot(20*log10(abs((fftshift(fft(hm))))))

Any ideas?


texts: How to plot BER Vs SNR of a received signal

texts:
I am using USRP to transmit and received MIMO frames similar to Wi-Fi standard that I created on my own. After receiving the signal, it is processed as described in the standard.

While I plot the QPSK constellation points after frequency and phase correction, I have a plot (shown above) which is noisy. Now I want to plot (BER Vs SNR).  I know that we can calculate BER as: (Error Bits/Total No. of bits). But I am not sure how to calculate SNR from these complex symbols.

What I understand from online materials is that SNR is the noise variance that is estimated from around the constellation points (maybe some distance/deviation measurement form the true point. Or its is given by Eb/N0, where Eb is the energy of bits received and N0 is the noise. The question might be very silly, what I do not understand is how I should find Eb and No from these complex received values that is plotted on a constellation diagram. I am doing all the coding in python.


texts: Difference equation system identification

texts:
Suppose I have a difference equation modeling a vehicle like this:
$$d[k+1]=d[k]+a\cdot u[k]+b,\tag1\label{eq}$$
where $d[k]$ is total distance traveled at time $k$, $u[k]$ is engine input at time $k$ (e.g. some measure of engine exertion at that time, not really important exactly what it is), and $a, b$ are parameters that I want to estimate from data points that I measure, for example by least squares.
From the difference equation I see that the units of $a\cdot u[k]$ and $b$ need to be distance. But suppose that the data I gather are actually samples of the vehicle's velocity (not position) for some range of input values, so I can plot velocity vs. engine input.
My question: I'm confused because it seems that it wouldn't make sense to do least squares regression to fit $a\cdot u[k]+b$ to these points since that would violate the units, but that is what is done in an example in my course. After the parameters are determined, the example then uses those same parameters to model the distance traveled over time via $\eqref{eq}$, going off of some initial starting distance $d[0]$. Am I missing something in this example?


texts: Frequency response of each component of a system given its global response

texts:
Given the following block diagram, find the frequency responses $H1(f)$ and $H2(f)$. The frequency response of the whole system has to be $H(f)=(\alpha_0+\alpha_1e^{-j2\pi T_1f}+\alpha_2e^{-j2\pi T_2f})^{-1}$

The fact that there's a loop confuses me. I would express $Y(f)$ as $Y(f)=X(f)H1(f)+(X(f)-X(f)H1(f)H2(f))\ H1(f)+...$, but it doesn't seem to be correct. Could you give me some hints? Thanks in advance!

  Solution: $H1(f)=\alpha_0^{-1},\ H2(f)=\alpha_0^{-1}\alpha_1e^{-j2\pi T_1f}+\alpha_0^{-1}\alpha_2e^{-j2\pi T_2f}$.


texts: Find original Hz of converted audio-file?

texts: I have received an audio file which is 44.1kHz wav. My ears suspect this audio clip was originally lower than 44.1, and encoded to 44.1 later, much like re-encoding a 720p video to a 1080p video. Is there a way to find out if this is true? Like, analyse the audio clip to see.. how many Hz it uses..? I have literally no idea how this works.


texts:
I have been trying to wrap my head around CIC filters.

Say that I want to do filtered decimation by a factor D.

Option#1: boxcar filter. Each output sample is formed by taking the sum of D input samples, and multiplying/shifting to get a proper average. D-1 adds per output sample.

Option#2: 1-stage CIC filter with an accumulator running at the input rate, and a comb filter running at the output rate. D adds per output sample.

So what is the fuzz about? Does CIC only make sense for higher orders? If order N CIC is like convolving the rect() function N-1 times, that would be a B-spline, right? I have seen authors praising B-splines for their functional benefits. Being able to do them efficiently might be handy.

Thanks for any input
Knut

Edit: pseudo matlab code showing the alternatives as I understand them.

N=100;
x = randn(N,1);
D = 5;

%% decimated moving average, ~100 adds
for k = 1:N/D
  y(k) = (1/D)*(x(k*D) + x(k*D-1) + ... + x(k*D-D+1));
end

 %% recursive moving average, ~200 adds
y(1) = x(1);

for k = 2:N
  y(k) = x(k) + y(k-1);
end

for k = D:N
    z(k) = y(k) - y(k-D);
end

z = z(1:D:end);

%% CIC, ~120 adds
y(1) = x(1);

for k = 2:N
  y(k) = x(k) + y(k-1);
end

z = y(1:D:end);

for k = 2:N/D
  w(k) = z(k) - z(k-1);
end


texts: How does projecting a few IR dots capture an image three dimensionally?

texts:
The original Microsoft Kinect device IR projected a spacially encoded dot field and looked for distortion in the dot field with a single IR camera to create a depth map. But there were thousands of dots.

I've seen a few devices however, usually in motion capture, that only use a few dots. The one pictued below for example is used to determine the position, movement and rotation of a moving golf ball.

Do they use a similar process of looking for distortion of the projected pattern, or some other process?


texts: Why is the ROC of Laplace transform independent of imaginary part of s?

texts:
An integral is defined as converging if it yields a finite value upon application of limits of integration. It is divergent otherwise.

Now sticking to the mathematical notation of Laplace transform, we have for a causal function $x(t) = u(t)$:
$$
X(s) = \int_0^\infty x(t)e^{-st}dt = \int_0^\infty e^{-st}dt
$$

$$
X(s) = -\frac{1}{s}e^{-st}\Biggr|_{0}^{\infty} = \frac{1}{s} = \frac{1}{(\sigma + jw)}
$$

By multiplying numerator and denominator with $(\sigma - jw)$, we get:
$$
X(s) = \frac{(\sigma - jw)}{(\sigma^2 +\omega^2)}
$$

For laplace transform not to exist, the denominator must become 0. Hence in this contrived example, both $\sigma$ and $w$ must be 0.

Conversely, if $\sigma = 0$ and $w \neq 0$, $X(s)$ exists as $\frac{-j}{\omega}$ with a magnitude of $\frac{1}{\omega}$.

Unless my basics are messed up, why in the literature is $jw$ disregarded as influencing $X(s)$? In other words, why is the ROC dependent only on $\sigma$?


texts: Alternative for pcshow() function in matlab

texts:
I'm using Matlab R2015a and I want estimate depth
but an error occurs and I wonder that it doesn't know pcshow() function. Why it doesn't work and is there another alternative for that?


texts: Mapping Saw Wave That Is Between -1 And 1 To A Sine Wave

texts:
I am trying to map values of a saw-wave that goes from -1 minimum value to 1 maximum value as a periodic signal onto a sine-wave.
Repeating the matrix of this waveform that is 257 element in length oscillates as expected as a saw-wave, but if I change the function into a sine-wave -1 maps to 0, then 1 also maps to 0, and when two zeros are together when we loop the waveform with "repmat" matlab function, there appears discontinuities in the sine-wave function...How can I map this saw-wave, so the sine-wave perfectly oscillates with all the values used from the saw-wave that is between -1 and 1 values?

f=(-1:1/128:1)';

c1=repmat(f,1000,1);

c2=sin(pi*c1);   % -1 of saw wave maps to 0, 1 maps also to 0

sound(c2,44100); % two adjacent 0's cause discontinuities

I also tried the fft method which gives me the needed coefficients but if my waveform is complex this would become more complicated to adapt:

f=(-1:1/128:1)';

m=fft(f);

m(3:256)=0;
m(1)=0;

m2=ifft(m);

m3=m2/max(abs(m2));

m4=repmat(m3,1000,1);

sound(m4,44100);


texts: How should we conduct a perturbation test for an LTI unknown system by perturbing the input to see the output response i.e. How many different perturbation step size should be given to the input?


texts: How to find minimum length of a FIR symmetric filter if I am given DTFT

texts:
I am practising for upcoming exams and came across this question.
Let $h[n]$ be an FIR filter such that $h[n] = 0$ when $|n| > M$ and $h[n] = h[−n]$. A plot of $H(e^{j\omega})$ (DTFT of $h[n]$) is shown below.
The eventual goal is to implement this filter digitally. Then, with the given information, is it possible to determine the minimum
value of $M$? Yes or no, with a suitable technical explanation.
What I though was that if I knew the order of this filter I could have found the length of the FIR filter as well. However order cannot be extracted from this graph.


texts:
Say I want to use Kalman filters for predicting the price of items at a supermarket. I have a Kalman filter for each item (apple/beef/brooms/etc). I notice that some items are sort of related, like apple with orange, beef with pork.

My question is how can I make these Kalman filters interact with each other? For example let's say I get observations about the price of beef and they show that they are going up. I want to tell the Kalman filter for pork somehow that it should also update. How can I do this?

Any literature or examples of this usage will be greatly appreciated, thanks

EDIT:

I am actually trying to ask this question in the general sense, but I wanted to frame the problem in a way that would be simpler for others to understand. So I understand that the Kalman filter uses observations to update its state. Suppose I want to incorporate observations for some Kalman filter into another, is there an accepted way of 'transforming' these observations? Or is it more of an art?


texts:
I am wondering about the correlation between input size and number of coefficients given by a discrete wavelet transform.

I am using Daubechies wavelets to describe a 1D function and I'm using PyWavelets to implement it (which is analogous to the MATLAB toolbox).

I started by implementing it using Haar wavelets, which gave correct results and I understand exactly how it works. Let's say my input function has 16 datapoints, if I use Haar, what I get from a multilevel decomposition (wavedec) is something like this (the number of shifts in brackets):

V1[1], W1[1], W2[2], W3[4], W4[8]

This is all well and good. The V1 gives me the scaling function and the W1-W5 wavelets of different scale and dilation. My problem is when I use the next Daubechies (referred to as 'db2' in the toolbox, which is called the D4), and I get

V1[6], W1[6], W2[9]

I lose all my intuition. I have no idea where 6, 6 and 9 come from, and they change depending on the level I specify (not even sure what it means to specify a level) and of course the input size. How can I predict the number of coefficients, and what are some good resources for gaining better understanding of why?

Thanks!

EDIT: Clarification on V and W:

$V_n$ usually denotes the span of a certain scaling function, $\phi$, i.e. $\{\phi_{n,k}\}$, where $k$ is the shift and $n$ the scaling. $W_n$ is the same except for the wavelet function, $\psi$. I might have abused the notation a bit by referring to the vectors of coefficients by V and W though.

EDIT2: Code

Here is the MATLAB code to produce the above:

>> [C, L] = wavedec(1:16, 4, 'db1'); L
L =
     1     1     2     4     8    16
>> [C, L] = wavedec(1:16, 2, 'db2'); L
L =
     6     6     9    16

I actually used PyWavelets, where it looked like this:

>>> import pywt
>>> map(len, pywt.wavedec(range(16), 'db1')) # defaults to level = 4
[1, 1, 2, 4, 8]
>>> map(len, pywt.wavedec(range(16), 'db2')) # defaults to level = 2
[6, 6, 9]


texts: Why do we use the HSV colour space so often in vision and image processing?

texts: I see the HSV colour space used all over the place: for tracking, human detection, etc...  I'm wondering, why?  What is it about this colour space that makes it better than using RGB?


texts: Detecting presence of music in ambient sound

texts:
I am hoping to create an application that would listen to ambient sound and detect if music is being played. It is not important to identify the music being played; just detecting that some music is being played is enough.

I looked around for existing solutions but couldn't find any. Does anyone know algorithms that I can use to solve this problem? If source code is available, all the better. I found a couple of academic papers and implemented solutions suggested in them. But the results I obtained were not satisfactory.

PS: i) It would be a bonus if the algorithm is not computationally intensive; if algorithm is completely in time-domain that would be wonderful. ii) It is okay if the solution is not very accurate; occasional false-positives are okay.


texts:
I have gone to wavelets course and up to my knowledge, changing the scale and translation of wavelets produces multiresolution images.

is this correct or can someone tell me what is multiresolution analysis in images?

Does multiresolution directly describe wavelet or multiresolution can be achieved other than wavelets?


texts: Calculation of derivative projection of differential Radon matrix

texts:
I am trying to write a matlab code  for finding the differential Radon Matrix.
.

Above are formula  for Radon Transform and differential Radon Matrix  .
What I get is the formula for Radon Matrix

R = radon(I)
R = radon(I,theta)
[R,xp] = radon(___)

But I unable to figure out how to use this code for calculating differential Radon Matrix


texts: How to perform dechirping of SAR image in MATLAB?

texts:
I have the transmitted pulse and the raw SAR image before range compression (dechirping). I understand the pulse compression its a convolution along the columns but I am not able to make it work. I've tried this:

   h_compressed = conv2(fliplr(conj(g)), 1, Dr, 'same');

Where g is the transmitted chirp and Dr is the raw SAR image.


texts:
How can I compare the way that word is pronounced with a sample? For example, most of dictionaries have sample word pronunciation recorded. The goal is to provide feedback to the user so that he or she can improve pronunciation.

I wonder if there are API available for such task or existing programs that are able to do it.

I found similar question here How do I compare 2 audio files, returning a percentage of difference , but I am providing a definition what actual difference is.

Regards,
Eugene


texts: Accelerometer drift: What hardware-level signal conditioning operations are being performed under the hood of a MEMS accelerometer chip?

texts:
I'm in a debate with a peer who says that filtering accelerometer signals at the chip-level has nothing to do with mitigating the problem of drift
Often, additional software-level filtering is employed to smooth noise, depending on the application. But My question pertains to any filtering or processing that is happening on the chip, specifically as it pertains to drift (and not so much to smoothness of the output).
My assertion is that noise due to thermal, electrical, vibrational, and sampling jitter will invariably contribute to drift, and that this type of noise is filtered on the chip to help reduce the problem. By reading an Invensense product PDF, we know there's "signal conditioning" happening on the chip. But the documentation makes no real mention of what's going on in there. It seems to be proprietary, and reasonably so.
I don't need any proprietary information, but if anyone could speak as to whether noise filtering reduces drift, that would be tremendously helpful. It seems obvious to me, but perhaps I am mistaken.


texts: Basic difference between Multi-scale and Multi-resolution images

texts:
While working on Convolutional Neural Network (CNN) I saw many research articles working on the Multi-scale images. In the literature of image processing I saw there is also Multi-resolution analysis of image which is at different level of images. Whereas in Multi-scale it seems that size of images are not changed instead it gets more blur and edges are getting more clear of the objects.
I am not using any wavelet transform. So in terms of general images how can we explain to beginners about the basic difference between Multi-scale and Multi-resolution images?


texts: Which one of these signals will be least distorted by this communication channel

texts:
So I was solving a problem that resulted in a communication channel that heavily distorts high frequency but does not significantly affect low frequency and I got stuck at figuring which one of these signals would be the least distorted. I know I have to find the one with the lowest frequency but I'm not completely sure how. These are the signals:

$$x_1(t) = \delta(t);\quad x_2(t) = 5; \quad x_3(t) = 10e^{j1000t};\quad x_4(t) = \displaystyle\frac 1t$$

I just need to know how to find the frequency of one of the signals


texts: Relationship between fractional delay coefficients and frequency?

texts:
Some background: I am putting a sig-gen signal through an ADC, into an FPGA where I can modify it, out a DAC and up to a SpecAn.  All that works fine.

I am experimenting with fractional delay filters, coefficients, and IF values (which are just my sine wave coming in for testing).  I sample at 200MHz, and I have played with IFs of 25MHz, 50MHz, 66.667MHz and other in between.  What I do on the FPGA is take the inputted signal, and sum it with a delayed copy of a signal (think multipath).  The delay is being varied by using a fractional delay filter to create the copy (either 3 coefficients or 7 right now; all the below is for 7 coefficients).

My question is, what is the relationship between the number of coefficients I use, and the frequency of the IF?  I think I am missing something fundamental (and probably easy), but I can't see what.

Here is a plot of the summations of the 25MHz signal with 7 coefficients on the fractional delay.  The y-axis is the signal level (dBm) and the x-axis is the amount of delay.  I stepped from 0 clock delay to 12 clocks (60ns) in 0.4 clock increments (2ns).  You can see the in-phase additions at clocks delays 0 and 8, and the out of phase cancellations at 4 and 12 clock delays.

Next is the same sort of plot, but the IF is 66.667MHz and I go from 0 clock delays to 3 clocks delays (15ns) in 0.1 clock (0.5ns) increments.  Here we should be in-phase at 0 and 3 clocks delays, and out of phase at 1.5.

What you can see in the second plot is that there seems to be a lag in the prediction.  I actually have peaks just after where I would expect to see them (delays 0 and presumably 3), and my out-of-phase signal seems to appear after where I would expect it as well (1.5 clocks).  There also appears to be more unevenness and the signal isn't as steady as the first plot (which was nearly as perfect and smooth as I could expect).

I understand that the faster IF is going to be harder to predict, but I feel like there is a relationship that I am missing between the IF, fractional delay coefficients, and maybe even sample rate.

Can anyone enlighten me?


texts: What are the various definitions for $\rm SNR$, and their associated methods for measuring it?

texts:
The definition of $\rm SNR$ seems to be somewhat of a tower of babel in industry. What definitions of $\rm SNR$ are there (feel free to site application), and how exactly can it be measured for that applications?

My specific questions on $\rm SNR$ are:

How may we measure $\rm SNR$ for a communications system, if we have not yet been able to attain optimal bit sampling timing, and all we have to work with are all the signals in the receiver up to and including the envelope of the I and Q channels? See this post for context.
Once we attain optimal bit sampling, and softbits are attained, how best to measure $\rm SNR$ (or $E_bN_0$)? One way I use is:
$$
10\log_{10}\left[ \frac{\textrm{mean}\left\{\lvert s_n \rvert^2\right\}}{\textrm{var}\textrm\{\lvert s_n \rvert\textrm\}}\right],
$$

however I understand this is not suitable for low $\rm SNR$ cases. What other ways exist?


texts: Decoding Matrix Encoded Surround Sound Through Convolution

texts:
I wish to create a finite impulse response filter to decode a dolby digital matrix encoded 2 channel signal into 5 channels. These filters would then be used in a realtime pipeline on a Linux machine using jack and jconvolver to do the convolution.

The equation is as follows:

The reference for the above equation is: http://en.wikipedia.org/wiki/Matrix_decoder

It would seem this should be able to be done but I haven't the faintest how. Any help, and in particular working code that could generate the filters (Octave or other open source preferred) would be much appreciated!


texts:
I have encountered convolution of two different impulse signals.

x[n] = (1/2)^n . u[n-2] * u[n]
x[n] = u[n] * [n]

u[n] = discrete impulse signal
. = product operation
* = convolution operation

For the first one, I found this solution:

x[n] = 1/4 if n = 2
x[n] = 0 if n != 2

For the second one, I found impulse signal itself

Edit: Are my answers are true ? My professor told me that the answer for the first one is wrong, but he did not say the correct answer.


texts: Sinogram Fourier Analysis Doesn't Show Clear Bow Tie Shape

texts:
I'm working on finding the center of rotation of a set of test tomographic projections in order to perform 2D reconstruction. I'm trying to implement the  algorithm to find center of rotation as described in Vo et al. 2014, which involves Fourier transforming the sinogram and searching frequency space outside the "bow tie" double-wedge region for minimum artifacts that result from wrong center of rotation. However, when I Fourier transformed a sample sinogram generated from Shepp-Logan model, the resulting frequency space has lots of artifacts outside the double-wedge region even when center of rotation is correct. Through visual inspection I've also noticed that the amount of artifacts are sensitive to number of projections in the sinogram.

Below: FT of sinogram with 200 projections, from 0 degree to 180 degree

Below: FT of sinogram with 180 projections, from 0 degree to 180 degree

The FFT'ed sinogram in the paper shows distinct zero regions outside the double wedge area, which is mathematically true as it is proven here. However, I'm not sure why my FFT'ed sinograms do not show similar zero regions as in the paper.

This is my first time posting here, so please let me know what other information is needed for others to help me. Thank you.

Below is my code

P=phantom(256);

angs=linspace(0,179,200);

pjs=radon(P,angs);

pjs=rot90(pjs,3);

pjs_inv=fliplr(pjs);

pjshift_inv=circshift(pjs_inv,[0 0]);

sinogram=[pjs; pjshift_inv];

figure(1),imagesc(sinogram);axis image,colormap gray

ft=real(fftshift(fft2(sinogram)));

figure(2), imshow(ft),colormap gray,axis image,axis on


texts: Motion detection one threshold over dataset

texts:
We're working on a project where we have to detect the moments where motions happen in a medical scene. I have a data set of 20 videos that are closely resemble (Same camera & configuration used, same shooting/environment conditions etc...). Here is what I did for each video:

Compute the optical flow (Farneback) between 2 consecutive images (t-1, t)
For each couple of images (t-1, t), compute the Norm1 of the displacement vector of all the pixels and then get the mean of these values. Then, I represent each couple of images by the mean of the Norm1 of optical flow results.
By doing so, I could draw a signal that represents the variation of these values as you can see below (Example video #1)

We zoom on the first 100 seconds of this signal to see the following:

To be able to detect the motions in video #1, we can just use as threshold the value 0.2 so we can differentiate between motion and not. But to be able to handle similar cases, we have to set a threshold at application/project level not on a video level: the threshold should be applicable over all the videos and for any new video of the same type. Let's take a look over the same kinda signal for video #2:

Obviously, the threshold in this case is 0.25 which is different from 0.2, the one we used for video #1 . To have a robust solution, we shall have one threshold per application so my question is about how can we use these thresholds (choosed manually) to have a global threshold for this application?

P.S: We already tried to use the mean of all these thresholds but it didn't work.. It failed to detect some important motions. We're asking such a question here because we just want to know if there is a fundamental principle to follow in such cases.


texts: Electromyography: is this interference from the heart/how to remove it?

texts:
I'm using electromyography to evaluate the activity of lumbar erector spinae at various levels of effort.
For this I also intend to sample the activity at complete rest (in prone lying). In my first few attempts I noticed regulary occuring spikes in the raw signal.

could this be interference from the heart beat or is it something else?

In another question I found some advice on how to filter these artifacts out.
But since emg-sampling and signal processing are both new to me I wonder

if there are methods easier to understand and apply than independent component analysis?

In the question above it is mentioned that it's possible to record the heart signal simultaneously and remove it subsequently from the emg-signal.

Is it possible to record the signal of the heart in prone lying and with a normal EMG-system?

Here is the data from the graph. I sampled it at 2000Hz.


texts: What is the Reference in Control Theory?

texts:
I've started studying Control and I've come across some concepts that I'm finding  a bit difficult to understand.

For example, in the system showed here, there is a signal called "reference" as well as a feedback block. As far as I know, the feedback block has the function of getting the output as close as possible to the reference. But here is what I don't understand: what if the reference is, for example, an impulse, and I want the system to show a step when an impulse is at its input?

I'm getting confused with the difference between reference signal and input signal. Why would someone connect a feedback block in order to get the output to follow the reference? What's the advantage of doing this?


texts: Java - Convert Image to Binary (Two Colours/Black & White)

texts:
Is it possible to convert a .jpg image to binary using Java?

I have asked a similar question on StackOverflow but thought that maybe this may be a good place to ask, too.

I am currently able to convert an image to negative in Java, but am unable as of yet to figure out a way of converting it into binary.

Any suggestions - and even sample code - would be amazing.

Many thanks.


texts: Transmission capacity in "duobinary coding"

texts:
In "duo-binary signalling" how exactly is the transmission capacity increased?

Because , when we encode the binary bits, the bit duration remains the same, but the amplitude will change. Then how can the transmission capacity be better?


texts:
This is probably basic but as I am new to the field it confuses me a bit.
While looking at some solutions provided to a problem in the final step following happens:

$$\begin{aligned} \frac{1}{10}\sum_{l=-\infty}^{\infty}x[l]\sum_{k=0}^9 e^{-i2\pi(n-l)k\frac{1}{10}}
&=\frac{1}{10}\sum_{l=-\infty}^{\infty} x[l]10\delta_{10}[n-l]\\
&=\sum_{s=-\infty}^{\infty}x[n-10s]\end{aligned}$$

Now, according to the DFT table we have provided for the course, the transformation in the 1st step is defined like this

$$e^{2i\pi(k-k_0)/N}\longmapsto N\delta[k-k_0]$$

The delta with subscript was never defined before appearing in this task so searching online I found this:

$$\delta_N[k]=\sum_{s=-\infty}^{\infty}\delta[k-sN]$$
My questions are:
What exactly describes the $\delta_N[k]$?

In the 1st step regarding the delta, is one of the 2 sources wrong or could it be that in the transformation table it should be ''understood'' that the delta has to have the subscript depending on its period?
Is the last step just the result of applying the properties of the delta(with subscript) function on the x signal?


texts: Binarization, and then thinning/skeletonization

texts:
I have pictures of networks, like the one below, and my goal is to obtain the network's skeleton from processing those images.

My approach lies in two steps, first I convert grayscale image to binary image using local thresholding or Otsu method, and then a medianfilter (python function medfilt). Result is shown below.

Then, I use a thinning algorithm to extract the network's skeleton.
Here is the result of my implementation of the Zhang-Suen thinning algorithm.

Considering the quality of the first image, I'm pretty sure it's possible to do a lot better than that. I thus have two questions :

1) Taking the last image, how would you remove all the small lines perpendendicular to the ridges, and the small gaps, that are artefacts ?

2) What algorithm would you actually advise me to use for the above mentionned steps (conversion to binary image and thinning) ?

I'm mainly working in python.


texts: How to water level deconvolve a noisy signal if i have a zeros and poles file?

texts:
Thanks for your time and help!
I am working with Apollo project passive seismic experiment (PSE) data, and I have a large set of seismic records (on digital counts) and the corresponding file of poles and zeros.
I need to deconvolve these registers to obtain calibrated seismograms on SI units.
I have read that this can be accomplished with a water-level deconvolution of 10% but I am not sure of how can make this.
I ask for some orientation and if possible bibliographical references
I am working in MatLab


texts:
I need to draw waveforms for biometric data like ECG and EEG signals. When I have more samples than pixels at the X-axis, I need to draw a vertical line between the MIN and MAX sample-value for that period.

Now I was wondering, suppose the signal is highly volatile, and the MAX for the next point is lower than the MIN from the previous point, there will be a 'gap' in the waveform.

What is the correct way to connect these points? When drawn using standard polyline routines the pen will draw from the last position (MAX) to the start of the new line (MIN), which is a much larger distance, causing the signal look thicker than it would be if I just connected the previous MIN to the new MAX (if the trend is downwards), or the previous MAX to the new MIN (if the trend is upwards).

Is there a standardized way to do this, or is it just a matter of taste?


texts: How to implement simple noise filter based on the frequency domain information?

texts:
Brief introduction:
I am working on my university project. My task is to develop application that captures sound from mobile phone and removes/reduces noise (this simpilied version of task).

I never faced with DSP before, I spend a lot of time to figure out how to get audio signal from phone (no direct API is porvided), so I have only few weeks left to study DSP.
There are a lot of information on that field, my basic source of information was: http://www.dspguide.com/.

My work, to be accepted should include some work(algorithms) in the field of DSP (so I should avoid ready libraries).

What I need
simple algorithm(easy to understand) for noise reduction, I am very limited in time (I got only one day left, after that I should present my work)

refernces to resources for understanding FFT & inverse FFT (although I implemented that algorithms in code, I simply adopted code from the book)

what I have done
FFT, Inverse FFT
I tried to remove frequencies by applying FFT to signal in time domain and zeroing needed ranges in result frequency domain (ranges < 300 && ranges >3700 - as human voice is in range [300;3700] - so doing that I tried to remove sounds not related to human voice, but I got very bad results (https://stackoverflow.com/questions/24101814/why-ideal-band-pass-filter-not-working-as-expected)

Could you please suggest me how I can based on the results of FFT reduce noise (I don't need super efficient method, I need something that I can present in my work and explain how I did it)

Please help!


texts: Adding white noise to complex signal (complex envelope)

texts:
Good day everyone.
So, I run in to trouble - How to add white noise to signal samples?
For now I'm working with MSK signal - quadrature modulator and quadrature demodulator + Viterbi processor.  My signal samples at TX and RX side represents Q and I signals of modulator(real and complex part of the signal). Model works nice in ideal situation, without any noise and other interference.
So, I'm working with complex envelope of the MSK signal with amplitude equal to one. Currently, I would like to test my demodulator with some noise.
It looks like, I just need to add (+) samples from my noise generator to samples of transmitted signal. But I'm not sure if this is correct way,

Should I sum both Q and I channels to one sample of noise signal or not?
And I'm pretty confused how to determine signal to noise ratio after this. I would like to measure bit error ratio at some SNR points. As I read in some books changing $\sigma$ parameter of white noise generator will play a role here, but still how to determine SNR in numbers after this?


texts:
I'm working with cascaded second-order allpass filters. Each filter is characterized by conjugate pole/zero pairs. In a long (>50) cascade of identical filters, the pole locations of the overall transfer function begin to diverge from their initial, theoretical location. They seem to fan out in a ring around where the true pole should be. This occasionally leads to instability, as the poles can end up outside the unit circle.

I'm assuming this is some kind of round-off or accuracy error, but haven't had any luck locating resources on this problem. Can anyone point me in the right direction?

EDIT:

Here is a code snippet demonstrating the problem. Uses Numpy/Scipy and the zplane function found here:

http://www.dsprelated.com/showcode/244.php

b = [0.72986695907729304, -1.0592380171923862e-16, 1]
a = [1, -1.0592380171923862e-16, 0.72986695907729304]

bN = array(1);
aN = array(1);

Nf = 10

for nf in range(Nf):
  bN = array(convolve(b, bN), dtype=float32);
  aN = array(convolve(a, aN), dtype=float32);

z_casc, p_casc, g_casc = zplane(bN, aN)


texts: Origin of artifacts 12 Hz away from electrical line noise (50 Hz)?

texts:
I made electroencephalographic measurements from people while they were sleeping. When plotting the frequencies, I see that I do not only have an expected 50 Hz (and 150 Hz and 250 Hz...) artifact (from nearby electrical devices) but also smaller artifacts to either side that are approximately 12 Hz apart (38 Hz and 62 Hz, 138 Hz and 162, ...) (examples from two participants):

Can someone tell me what might cause those artifacts?


texts: SIFT About Difference-of-Gaussian function extrema?


texts: Why does a wave continuously decreasing in frequency start increasing its frequency past the half of its length?

texts:
I'm trying to programmatically generate a wave (sine or square) with a continuously decreasing frequency.
To do so, I simply multiply the starting frequency for a decreasing value, that linearly goes from 1 to 0, at intermediate steps throughout a given time span.
I was puzzled to see (and hear) that, exactly at the half of the time span, the frequency started to increase again.
I checked the frequency value but everything is ok.
I can use only the first 'half' of the wave, but it would be great to understand what I'm experiencing.
Any help appreciated !
Thanks


texts: Absolute value vs. power (to make a signal positive-only)

texts:
In order to make a signal positive-only, I have the choice between taking the absolute value and multiplying the signal by itself (i.e. taking the signal to the power of two).

Obviously, multiplying the signal by itself requires more circuitry (a multiplier) and the result has a larger dynamic range. Are there other aspects that I need to consider when chosing between the two methods?

Edit: the application is to detect negative and/or positive signal peaks just above the noise level.


texts: Linear phase FIR filter for impulse responses that don't appear symmetric

texts:
I would like to clarify some confusion I have about linear phase FIR filters of which do not seem to have symmetric impulse responses.

Starting with a simple case, a delay, $ h[n] =\delta(n-n_0) $ does not seem to have  have a symmetric impulse response. For example,

$$h[n] = [0, 0 ,0,1,0] $$ is not a symmetric impulse response. How is it that the property of symmetry for a linear phase FIR filter holds true, in this case?

For another example, say I have a filter defined as

$$h[n] = [0,0,0,0,1,0,0,0,0,0,1]$$

I believe this filter is linear phase too. But it's impulse response is not symmetric.

How does the symmetric property hold true? Does it have to do with zero padding? But even if we don't zeropad the filter, it is still linear phase? If we zeropad it by a large factor, is it also linear phase?

For one more example, say I have this filter:

$$ h[n] = [1, 2, 3, 4, 3, 2, 1]$$

This is obviously linear phase. When I zeropad it by a factor N, the impulse response is no longer symmetric, but it stays linear phase for all values of N. Why is this the case?


texts: Why is privacy scrambling a standalone block in a communications transmitter?

texts:
In some telecommunication systems, scrambling of transmitted data is applied to provide some level of security. For example, the communication system might look like:

Data -> Scrambling -> Compression -> Error correction encoding

Could one possibly match the scrambling with the error-correcting code to obtain some benefit? Why is scrambling typically separated from the other blocks of the transmitter?


texts: How to remove sum and frequency components in FFT spectrum during signal resampling.?

texts:
I am working on to remove non linear tuning of laser . I am using resampling technique based on triggering technique .
Firstly, I make collect a vector of maxima and minima from Tigger sine wave of frequency F1. Then I reconstruct signal sine wave (of frequency F2) with the maxima and minima vector from the trigger. FFT is performed after this step to identify noise free frequency components of the signal.
This kind of of resampling is working excellent. But after FFT,  besides the actual frequency component of the signal wave, it also gives rise other unwanted frequency components.
These are exactly sum and frequency components arising from F1 and F1,( i.e., F1+F2 & F1-F2 ).
For the clear understanding of my question please refer to attached image  and suggest me the possible solution to avoid these unwanted frequency generations.
thank you


texts: Mean Square Continuity of Random Process

texts:
Show that a stochastic process $X(t)$ is mean square continuous if and only if its autocorrelation function $R_X(t_1,t_2)$ is continous

$\Rightarrow$ Proof:

We have $E[(X(t)-X(t_0))^2]=R_X(t,t)-R_X(t_0,t)-R_X(t,t_0)+R_X(t_0,t_0)$, so if $t \rightarrow t_0$ the right side of the equation is equal to $0$, therefore $\Rightarrow$ implication is true.

$\Leftarrow$ Proof:

This is where I have trouble.

Is this implication true due to the fact that if $X(t)$ is continuous at $t_0$ in m.s then $\lim_{t \to t_0} m_X(t)=m_X(t_0)$ ??


texts: MATLAB flanger...zipper noise, crackling, help with interpolation values?

texts:
i'm writing a batch script in MATLAB to flange.

i'm getting some periodic noise, and i'm having trouble tracking down the source. it's definitely in the algorithm, because soundsc() will scale values larger than 1. and theres no memory or computing time issue, because it's running as a batch process.

any mp3 or wave file should do the trick

current script, with interpolation (sounds worse so far than non-interpolating, because i am probably using the wrong X for the interpolation.)

function result = dLinTerp(x1,x2,y1,y2,x)

% calculate decimal position of x
dx = (x - x1)/(x2 - x1);

% use weighted sum method of interpolating
result = dx*y2 + (1-dx)*y1;

and

    cur_frac=cur_sin*max_samp_delay;

    if (readIdx + cur_frac >= (max_samp_delay))
        fetchIdx = cur_frac - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_frac;
    end

    if (floor(fetchIdx) = max_samp_delay)
        fetchNext = 1;
    else
        fetchNext = ceil(fetchIdx);
    end
     integ = fix(cur_frac);
     frac = abs(cur_frac - integ);
    interp = dLinTerp(0.0, 1.0, delayBuff(fetchPrev), delayBuff(fetchNext),frac);
    out(i) = (amp*input(i)) + (amp*interp);

here is the original script, no interpolation

%flanger

[input,fs] = audioread('pianoriff_s.mp3');
input = input(:,1);
len = length(input);
out = zeros(len,1);

max_time_delay=0.007; % 7ms max delay in seconds
rate=0.2; %rate of flange in Hz
max_samp_delay= floor(max_time_delay*fs);
delayBuff = zeros(max_samp_delay,1);
amp=0.7;
readIdx = 1;

for i = 1:len,
    delayBuff(readIdx) = input(i);

    sin_ref = sin(2*pi*i*(rate/fs));
    cur_sin=abs(sin_ref); %abs of current sin val 0-1
    cur_delay=ceil(cur_sin*max_samp_delay);

    if (readIdx + cur_delay >= max_samp_delay)
        fetchIdx = cur_delay - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_delay;
    end

    out(i) = (amp*input(i)) + (amp*(delayBuff(fetchIdx)));

    if (readIdx+1 >= max_samp_delay)
        readIdx = 1;
    else
        readIdx = readIdx + 1;
    end
end

soundsc(out,fs);


texts: general n-dimensional name of Pixel and Voxel

texts:
In 2D images we call the smallest picture elements pixels, in 3D images we call them Voxels (for volume pixels).

I was wondering if there is a word that can denote them both (pixel and voxel), and possibly be used to mean any n-dimensional picture element.

For example, I want to write something like:

  Given an n-dimensional image with N pixels/voxels

Where I want to replace pixels/voxels with just one word that would also denote these elements in other dimensions.


texts: Instantaneous phase calculation problems

texts:
I am currently struggling with a problem when calculating the instantaneous phase of a wavefield using a 1D Hilbert transformation.

The scheme is as follows.

1D Hilbert transformation of wave field $U(x,z)$
$$
q_z(x,z) = H_z[U(x,z)] = \int_R U(x-\xi,z)\frac{d\xi}{\xi}
$$

Calculation of the instantaneous phase $\phi(x,z)$

$$
\phi(x,z) = arctan\left(\frac{q_z(x,z)}{U(x,z}\right),
$$
where $x$ and $z$ are the coordinates and $H_z$ is the 1D Hilbert transformation over the vertical direction.

My problem is, that $\phi(x,z)$ which is implemented using the intrinsic function $atan2()$ contains phase shifts [$\pi \rightarrow -\pi$].

Using the instantaneous phase for further calculations in which its gradient is derived hence leads to imaging artifacts as shown in the attached pictures:

TOP LEFT: input wavefield $U(x,z)$

TOP RIGHT: Hilbert transformed wavefield $q_z(x,z) = H_z[U(x,z)]$

BOTTOM LEFT: instantaneous phase $\phi(x,z)$

BOTTOM RIGHT: derived propagation angle with artifacts

.

My question: I am able to overcome this problem using Gauss filtering or Median filtering. Yet this does not wotrk properly for more complicated wavefields and higher frequencies. Is there another way to overcome those sharp transitions in the area of my wavefield? (not the noise around)

Thanks alot!


texts: Bandpass overshoots at the beginning and end of signal: why?

texts:
The picture below shows the original data in the upper figure and the filtered data in the figure below. The filter applied to the original data is a simple bandpass implementation:

$H(z)  = \frac{0.2 + 0.1z^{-1} - 0.1z^{-3}  - 0.2z^{-4}}{1-0.98z^{-1}}$  (Relative Spectral Processing Filter)

As you can see, at the beginning of the data and at the end the filtered signal overshoots / undershoots sharply.

So far I experimented with a) adding random noise to the original signal and b) removing leading / trailing zeros but those are just my own ideas without any reference for a thesis. So my question is:

How would this problem be handled in digital signal processing?

+++++++++++

If you're interested: Removing leading / trailing zeros reduced the problem (original overshoot is ~40, now it's around 12 which is still the highest value in the filtered signal) and adding random noise reduced the problem even further but the signal is noisy now (which  is a unfavorable tradeoff, isn't it?)


texts:
What is the Fourier series representation of a cycloid?

The parametric representation of the curve is as follows.
$$
t=\dfrac{\theta-\sin\theta}{\pi}\\
x=\dfrac{1-\cos\theta}{\pi}
$$
The period is $2$, so the coefficients of the complex exponential Fourier series should be
$$
c_n=\dfrac{1}{2}\int_0^2\!x(t)e^{-jn\pi t}\,\mathrm{d}t
$$
and that's where I got stuck. I performed substitution
$$
t=\dfrac{\theta-\sin\theta}{\pi}\implies\mathrm{d}t=\dfrac{1-\cos\theta}{\pi}\,\mathrm{d}\theta\\
c_n=\dfrac{1}{2}\int_0^{2\pi}\!\left(\dfrac{1-\cos\theta}{\pi}\right)^2e^{-jn(\theta-\sin\theta)}\,\mathrm{d}\theta
$$
but don't know how to integrate this monster. Is there any other way?


texts: Looking for cycles of periods longer than the input signal length

texts:
There was a similar question asked here, however I would like to focus on some specifics of this problem. Let me present the Python code sample to illustrate the situation when only a part of the period of the input signal is analyzed using FFT:

N = 1200
H = 4
fig = plt.figure()
sig = np.zeros(N)
for h in range(H):
        sig += [(h+1)*np.sin(2*np.pi*(h+1)*x/N) for x in (range(N))]
ax1, ax2 = (None, None)
for h in range(H):
        if (h):
                fig.add_subplot(H, 2, 2*h+1, sharex=ax1)
        else:
                ax1 = fig.add_subplot(H, 2, 2*h+1)
        plt.plot(sig[:int(N/(h+1))])
        if (h):
                fig.add_subplot(H, 2, 2*h+2, sharex=ax2)
        else:
                ax2 = fig.add_subplot(H, 2, 2*h+2)
        fft = pd.DataFrame()
        fft['coef'] = np.fft.fft(sig[:int(N/(h+1))])
        fft['freq'] = np.fft.fftfreq(int(N/(h+1)))
        fft['per']  = 1/fft['freq']
        fft['mag']  = fft['coef'].apply(lambda z: np.sqrt(z.real**2 + z.imag**2))
        rside = fft['freq'] > 0
        plt.bar(fft[rside]['per'], fft[rside]['mag'])
fig.show()

The parameter H specifies the number of composite sines, their magnitude and the number of divisions of the original signal. For each division we get 1 row of subplots: the input signal in the left column and the period spectrum in the right column. The X axes in each column are shared, so the plots are aligned. Only positive frequencies are considered. The example with H = 4 and the signal length at N = 1200 should produce the set of plots like on the picture below.

The FFT results in the top row are spot on, but on other plots from the right column we can see that the spikes representing the correct dominant cycles are missing, but some additional smaller spikes appear at periods much smaller than the smallest one used to generate the test signal.

In a general case when given a sample of discrete observations or a time series data I would like to have a way to test how would some bigger cycles fit in the spectrum, ie. if we consider the Mercury's distance from the barycenter in that other thread someone suggested that orbital cycles of other planets should be considered, all of which are well known and some are longer than 8 years (the size of the input data). I was searching online, but due to lacks in theoretical (and practical) knowledge on the matter I probably wasn't using the right keywords, so I couldn't find any helpful information. These are the questions I would like to ask:

Can the small cycles that FFT 'found' in the period slices provide clues about the bigger ones that were actually in effect, ie. could there be a correlation in frequency or phase between the small fake cycles and the real big ones?
Is there a way to expand the period domain or in other words assume there are bigger cycles at play and include them in the spectrum? Can Fourier Transform be somehow manipulated to do so?
Lastly is Fourier Transform the right tool or at least a starting point to deal with such issue?


texts: MATLAB - How do I get the Filter Function to get the output Energy Spectrum?

texts:
The image was from a question I posted on reddit, so ignore the title.

I'm trying to write a program that plots these graphs

I know that the first plot is just a plot of the source-signal so the "time" domain of the voice signal. I have that down.
the 2nd plot is the frequency domain signal, which I should be able to get with fft() fcn in matlab? but the plot it gave doesn't look right.

weird fft plot of source signal - me - a female just saying "aaaaahhhh" for about 15 seconds

And after I get the right filter function, I do conv??

Thanks in advance.


texts:
The rectangular function is defined as:
$$\mathrm{rect}(t) = \begin{cases}
0 & \mbox{if } |t| > \frac{1}{2} \\
\frac{1}{2} & \mbox{if } |t| = \frac{1}{2} \\
1 & \mbox{if } |t| < \frac{1}{2}. \\
\end{cases}$$

The triangular function is defined as:
$$\operatorname{tri}(t)
=
\begin{cases}
1 - |t|, & |t| < 1 \\
0, & \mbox{otherwise}
\end{cases}
$$
It is the convolution of two identical unit rectangular functions:
$$
\operatorname{tri}(t) = \operatorname{rect}(t) * \operatorname{rect}(t) \quad
=  \int_{-\infty}^\infty \mathrm{rect}(\tau) \cdot \mathrm{rect}(t-\tau)\ d\tau$$
Zero-order hold and First-order hold use these functions. In fact, it has:
$$
x_{\mathrm{ZOH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{rect} \left(t-n\right) \
$$
for Zero-order hold, and
$$x_{\mathrm{FOH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{tri} \left(t-n\right) \ $$
for First-order hold. Since $\operatorname{tri}(t) = \operatorname{rect}(t) * \operatorname{rect}(t)$, I would like to know if this is just a coincidence or if, for the Second-order hold the impulse response is
$$\operatorname{tri}(t)*\operatorname{tri}(t)=\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right)*\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right).$$
Is it true also for a general $k$-th order hold? Namely, put
$$x_{\mathrm{K-TH}}(t)\,= \sum_{n=-\infty}^{\infty} x(n)\, \mathrm{g}_k \left(t-n\right) \ $$
where $\mathrm{g}_k \left(t-n\right)$ is the impulse response of the $k$-th order hold, I would like to know if its impulse response is
$$\mathrm{g}_k \left(t-n\right)=\left(\operatorname{rect}(t) * \operatorname{rect}(t)\right)* \dots * \left(\operatorname{rect}(t) * \operatorname{rect}(t)\right),$$
k times.


texts: Two different definitions of cross-correlation

texts:
I have come across two different definitions of the cross-correlation function:

$$
R_{XY}(\tau) = \int_{-\infty}^{+\infty} x(t) y(t+ \tau) dt
$$

and

$$
R_{XY}(\tau) = \int_{-\infty}^{+\infty} x(t) y(t- \tau) dt
$$

Are these two definitions equal?

Also is $ R_{XY}(\tau) = R_{YX}(\tau)$?


texts: Algorithm for identifying beat grid in music with non-constant tempo

texts:
Are there any algorithms for automatically fitting a beat grid to a song (tempo matching) when the song has a slowly changing tempo or BPM?

The output would be a vector of onsets.

Here's a picture of fitting a grid to a waveform for a part of a song with a constant BPM:


texts: Corner detection using Chris Harris & Mike Stephens

texts:
I am not able to understand the formula,

What is $W$ (window) and intensity in the formula mean,

I found this formula in opencv doc http://docs.opencv.org/trunk/doc/py_tutorials/py_feature2d/py_features_harris/py_features_harris.html


texts: what is nyquist rate of $h(t)\cdot h(t)$ and $h(t)\circledast h(t)$

texts:
Let's say we have $h_c(t)$ as a continuous-time signal with bandwidth $B$ and we would like to sample it.
To be able to reconstruct it correctly, the sampling rate must be greater than $2B$.
Now assume we have $h_{c1}(t)=h_c(t)\cdot h_c(t)$ and $h_{c2}(t)=h_c(t)\circledast h_c(t)$ what is Nyquist rate of these two signals? Is it $4B$ for both of them because each $h_{c}(t)$ needs $2B$?


texts: How to calculate the Fourier transform of absolute-valued sinusoid?

texts:
The signal $x(t) = |1+a\sin(\omega t)|,(a>0)$ is a continuous waveform. In order to extract the frequency parameter $\omega$, I conduct the FFT of it and obtain its spectrum showed as follows. The frequency $\omega$ is set as 1 Hz. Spikes occur at the fundamental frequency and its harmonics. I'm curious that why this phenomenon exists. How to calculte the Fourier transform of $x(t)$?

Spectrum of $X(f)$

The waveform of $x(t)$


texts: Inverse DTFT of phase shifted complex exponential

texts:
I have been working on this problem for a few days now and I think this is the closest I have gotten.
I am getting an Answer of zero and I would like to know if that is correct and if someone could check my work please? Thanks for your time, here is what I have so far (also Im quite sure evaluating at n = -10 isnt correct, could someone please give me the correct process?)
Not looking for an answer, just how to start/continue to process. Thanks
I put all my steps in the attatched picture. (sorry, Im new to Stack Exchange and couldnt figure out how to format correctly)
Thanks for your time :)
-Dom


texts: What is the meaning of positive and negative shift in signal correlation?

texts:
If $x(t)$ and $y(t)$ are real valued signal then from this page cross correlation between them is defined as

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau) dt\quad\text{+ve shift}\tag{1a}\\
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\quad\text{-ve shift}\tag{1b}\\
\end{align}

I tried to solve a problem where $x(t) = e^{-at}u(t)$ and $y(t) = e^{-bt}u(t)$.

  From $(1a)$: Positive shift

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau)dt\\
&= \int_{-\infty}^{\infty}e^{-at}e^{-b(t-\tau)}u(t)u(t-\tau)dt\tag{2}
\end{align}
For $\tau\geq0$. Limits of integration will be from $\tau\leq t\leq\infty$. Hence from $(2)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{\tau}^{\infty}e^{-at}e^{-b(t-\tau)}dt\\
&=e^{b\tau}\int_{\tau}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{-a\tau}}{a+b}
\end{align}
For $\tau\leq0$. Limits of integration will be from $0\leq t\leq\infty$. Hence from $(2)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{0}^{\infty}e^{-at}e^{-b(t-\tau)}dt\\
&=-\frac{e^{b\tau}}{a+b}\int_{0}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{b\tau}}{a+b}
\end{align}
Therefore
\begin{align}
R_{XY}(\tau)=\begin{cases}
                    \frac{e^{-a\tau}}{a+b}; \tau\geq0\\
                    \frac{e^{b\tau}}{a+b}; \tau\leq0
             \end{cases}\tag{3}
\end{align}

  From $(1b)$: Negative shift

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t+\tau)y(t)dt\\
&= \int_{-\infty}^{\infty}e^{-a(t+\tau)}e^{-bt}u(t+\tau)u(t)dt\tag{4}
\end{align}
For $\tau\geq0$. Limits of integration will be from $0\leq t\leq\infty$. Hence from $(4)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{0}^{\infty}e^{-a(t+\tau)}e^{-bt)}dt\\
&=e^{-a\tau}\int_{0}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{-a\tau}}{a+b}
\end{align}

Corrected my analysis after @Fat32 answer

For $\tau\leq0$. Limits of integration will be from $-\tau\leq t\leq\infty$. Hence from $(4)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{-\tau}^{\infty}e^{-a(t+\tau)}e^{-bt)}dt\\
&=e^{-a\tau}\int_{-\tau}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{b\tau}}{a+b}
\end{align}
Therefore
\begin{align}
R_{XY}(\tau)=\begin{cases}
                    \frac{e^{-a\tau}}{a+b}; \tau\geq0\\
                    \frac{e^{b\tau}}{a+b}; \tau\leq0
             \end{cases}\tag{5}
\end{align}

Hence from $(3)$ and $(5)$ as can be seen the result match. If my analysis is wrong do correct me and if its correct then what does these negative and postive shift signifies?

On substituting $u = t - \tau \Rightarrow t = u + \tau \Rightarrow dt = du$ in $(1a)$, we get $(1b)$ as
\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(u+\tau)y(u) du = \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\\
\end{align}

Hence

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau) du = \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\\
\end{align}


texts:
Practical Infinite-Impulse-Response (IIR) filters are usually based upon analogue equivalents (Butterworth, Chebyshev, etc.), using a transformation known as the bilinear transformation which maps the S-plane poles and zeros of the analogue filter into the z-plane. However, it is quite possible to design
an IIR filter without any reference to analogue designs, for example by choosing appropriate locations for the poles and zeroes. Can somebody explain the latter design  of digital IIR filters (i.e. without any reference to analogue design) for the following simple example?
For a digital system with sampling frequency of 60 MHz, design a digital IIR filter with two complex conjugate poles at 23 MHz, and one zero at 18 MHz.
This is basically an equalizer for a lossy channel. The filter is flat at lower frequencies (DC attenuation), with a peaking at higher freq and then drops rapidly. For that, only knowing the poles and zero locations should be enough which defines the DC attenuation, bandwidth, and boost (peaking) of the filter. The amount of boost or DC attenuation does not matter as they can be tweaked by changing the poles and zeros locations. So I don't think any further information is required here. But if so, simply make an assumption.


texts: Why can Sum of Squared Differences (SSD) seen as cross-correlation?

texts:
In this slides, on page 4,  the SSD can be seen as cross-correlation:

http://www.lira.dist.unige.it/teaching/SINA/slides-current/pyramids.pdf
http://www.cse.yorku.ca/~kosta/CompVis_Notes/ssd_cross-correlation.pdf

The SSD is defined as:

$$SSD(m,n) = \sum_i \sum_j [g(i,j) - t(i-m,j-n)]^2,$$

and expand the SSD. It looks like this:

$$SSD(m,n) = \sum_i \sum_j [g(i,j)^2 + t(i-m,j-n)^2 - 2g(i,j)t(i-m,j-n)], $$

where the squared terms can be assumed as constant. As a result, the cross-correlation remains:

$$R(m,n)= \sum_i \sum_j [g(i,j) t(i-m, j-n)].$$

Question: Why can we assume that the squared terms are constant?


texts:
I have this recording from my DAW that showed a strange form of distortion: the signal clipped at many places, but instead of simply being clipped (loss of information), somehow the integer value of the signal overflowed and has been stored in the file. It resulted in worse-sounding distortion, but the information is still there, as a mirrored image in the opposite side of the signal.
This is plainly evident when looking at the waveform. You can clearly tell when the signal clips because the rest of the wave is visible between the peaks.
How should I approach fixing this file? I am very comfortable with SOX, ffmpeg, the C language, SuperCollider and various other tools.


texts: Reducing latency due to FFT windowing for a real-time signal

texts:
I need to compute the amplitude of a certain frequency band as the input to a PID algorithm.  My system looks like this:

[Sine Generator] -> [Actuator] -> [Physical complex coupling] -> [Sensor]
  ^                                                                 |
  |                                                                 |
  +--[PID Control Algorithm] <---- [FFT Amplitude computation] <----+

The PID controller adjusts the amplitude of the sine wave until the sensor and FFT register the correct amplitude in the output signal.  In order to have a performant system, I need to compute the amplitude as quickly as possible so that the PID controller can make adjustments based on the most recent events.

In normal analysis, I would use data from the future to calculate the amplitude at a given point in time, but that information is not available when I need it.  As a result, the input data is all from the past. (Also, during the initial start-up of the system, I have no data at all!) The number of missed cycles is dependent on my window size - for accuracy, I would like it to be as large as possible, but for latency, I would like it to be as small as possible (I have abundant processing power, that's not  a significant constraint yet).

Due to windowing, the dominant samples used as input to the FFT are several cycles before the present time and the most recent, most important samples are only very slightly weighted.  The obvious solution would be to adjust the windowing function such that the most recent cycles are heavily weighted - essentially, using just one side of a windowing function.  However, this does not seem to preserve the desired effects of the windowing function, as explained (somewhat circuitously) in this question: FFT with asymmetric windowing?.

I've also tried inverting and mirroring the data, which works when the input signal is a pure sine wave (of varying amplitude) at a zero-crossing, but introduces artifacts at other times.  I've begun looking into the Goertzel algorithm and sliding DFTs, but those seem to only improve the difficulty of the computation and not the latency of the result.

How can I reduce the lag associated with sampling a real-time signal?


texts: Fringe movement analysis to extract displacement

texts:
I am using a Michelson interferometer with the interferogram recorded on a CMOS sensor. I have taken a video of the fringes moving when a displacement of 50 $\mu m$ is imposed to one of the mirrors.
The frame rate is 17.9980 and it is already the maximum I can have.
For a 50 $\mu m$ displacement I should see around 79 fringes passing on the sensor, as each fringe that passes corresponds to one full wavelength (632.8 nm). However, I cannot obtain clear "peaks" that allow me to identify if a fringe is moving.
In the picture you can see the frame with the fringes (top left) and the signal in the central pixel column (top right) for each frame. Then, I transformed the signals using FFT in Matlab, and I obtained the spatial frequency (always the same, so it should be correct) and the phase. Theoretically, I should be able to extract the displacement in two consecutive frames checking the phase difference, so that $$\Delta x = \frac{\lambda}{2} \frac{1}{2\pi}\Delta\phi$$
and the total displacement being the sum of the $\Delta x$ for each couple of frames.
Alternatively, I thought that I could extract the displacements looking at the integral of the signal $$\Delta S_k = S_k - S_0 $$, which is the difference between the signal at the $k$-frame and the initial signal at frame $0$. When the two signals are in phase, the integral of the difference signal should be close to zero, while when the two signals are not in phase the integral increase to a maximum when the $\Delta \phi=\pi$ (bottom right)

In the following picture there is the wrapped phase (top), unwrapped phase (centre) at the main frequency, and the integral of the difference signal at each frame (bottom). Imposing a 50 $\mu m$ in the same direction, the sign of $Delta \phi$ should be always the same, but apparently it is not. As you can see, from the difference signal the peaks are not easy to determine, there are "double peaks" and small peaks. Theoretically, I should find 79 peaks but I got more.

Is there an easier way to count fringes?
I guess that the peaks would be easier to identify if I move the mirrors more slowly?


texts: How to overcome effect of intensity variation in image processing?

texts: I want to process an image for which I need the method to be robust to lighting conditions. For that the only thing I can think of is normalising intensity. I did this by dividing the image with average intensity and then scaled the values to be between 0-255. This method doesn't seem to work well. What other approach can I follow to make my processing robust to variation in lighting conditions?


texts: Deconvolution of Synthetic 1D Signals - How To?

texts:
I convolved a square wave with a Gaussian wave using linear convolution. Can I get the original square wave back by deconvolving my output with the Gaussian function?

I took the FFT of both signals, divided and then took the IFFT to get back the square wave. The output I got looked like random noise. I thought this might be due to the fact that I am doing a division and some denominator values might be very low, hence causing this error. I tried to set up a threshold and obtained the result depicted below.

How can I improve the result?

Edit 1 : Adding the code used to generate Gaussian signal. Was using matlab.

fs      =  200;           % Sampling frequency

t1 = 1/fs : 1/fs : n1/fs; % Where n1 is the length of the square wave signal

s2   =  4*gaussmf(t1, [ 0.4 4.5 ]);   % Generating Gaussian signal

Edit 2   :  Adding code used to generate square wave. I am afraid, it's a clumsy code!

fs      =  200;                           % Sampling Frequency
t       =  1/fs : 1/fs : 3;
n       =  length(t);

sq(1 : round(n/3) )  = eps;
sq( (round(n/3) + 1) : 2*round(n/3) ) = 3;
sq( (2*round(n/3) + 1) : n ) = eps;

sq( (3*round(n/3) + 1) : 4*round(n/3) ) = 3;
sq( (4*round(n/3) + 1) : 5*round(n/3) ) = eps;
sq( (5*round(n/3) + 1) : 6*round(n/3) ) = 3;

sq( (6*round(n/3) + 1) : 7*round(n/3) ) = eps;
sq( (7*round(n/3) + 1) : 8*round(n/3) ) = 3;
sq( (8*round(n/3) + 1) : 9*round(n/3) ) = eps;


texts:
I'm trying to pick a method to increase the bin-width on several discrete Fourier transforms I'm incorporating into some scientific analysis code. I have several questions about how this could work using Welch's method, and in particular how that would affect the values of the spectral densities and correlation functions one could compute from the Fourier transforms.
Because I'm writing this in C++, and because I care about the amplitude of the spectral density detected (not just correctly detecting a particular frequency) I'm thinking it would be better not to use a windowing function with the sub-transforms; in particular, if the amplitude at a bin matters, then I have to figure out how much averaging the different windowed periodograms together changed the frequency in question, which maybe is simple for someone using functions from a DSP library but may not be so simple for me. If this is a bad idea, or in particular if you're really not supposed to use Welch's method without a windowing function, I'd appreciate understanding that better.
Are there any issues with simply solving for Welch subsample length and number, given a fraction of overlap between Welch subsamples and a desired window width? Aside from the computational expense of this process being less and less rewarding as more resamples are added, are there other issues with shortening the transform in this way?


texts:
I would like to know when I can determine that a system is stable. I've seen Stability of a system but this is very theoretically.
My background is that I've set up a measurement system to detect the signal changes of a gas sensor. Means, when I add such and that of a gas the sensor will change accordingly.

But when do I state that the signal I observe is stable? Per eyeballing it isn't that hard but I would like to use a better, more reproducible characteristic. A common thumb of rule is fine.


texts: A system that perfoms Fourier Transform operation - is it an LTI system?

texts:
If a system takes input as the time domain signal and outputs the frequency domain signal, is such a system an LTI system?
 For if the input time domain signal can be represented as a linear combination of sinusoids, then since the operation is transformation to the frequency domain, the outputs will also be the same Linear Combinations of the individual outputs (frequency domain representation of the individual sinusoids).. right? so we can say it's Linear?
 Since the system performs the same operation independent of time, it would be Time-Invariant too right?


texts: To estimate the uncertainty of Slicing holes from an image?

texts:
I am thinking to estimate the uncertainty of holes photographed from the bottom and top side.
The difference of the same hole photographed from two sides is distributed as Gaussian.
The length of the plane is $10.00 \pm 0.01 cm$ where are holes distributed hexagonally with $50 \mu m$ difference in distance.
I have attributed these factors in affecting the uncertainty of the hole in photographing it in small magnitudes:

length of the plane
viewing angle of the camera (by Taylor expansion getting 3.6% systematic uncertainty)
slicing of holes from the plane

I am uncertain how much can the algorithmic approach of detecting the border of the hole in the image affect the uncertainty of the diameter of the hole.
Some steps of algorithmic approach:

Canny algorithm for edges so allowing finding circles
circle's diameter can be measured as maximal distance between two points in one circle

The hole plane of 10cm is divided into smaller regions for the algorithmic approach, like the following for one sample where the Canny is applied:

which is a simplification, since the original data is 26.1 Mb (tiff) in a single sample picture.

How can you approximate the uncertainty in slicing a hole from a plane?


texts:
I have found this interesting question here:

Which of these two waveforms has a richer spectral content, and why?

Intuitively speaking, I think that the first one is richer, but I cant really justify it! Cheers for any input.


texts: Understanding why (at least) 2 points of the FFT of a real signal are real

texts:
I am aware that the FFT of a real signal is Hermitian, i.e. $\text{FFT}(f)[i] = \text{FFT}(f)[-i]^*$, where $\text{FFT}(f)[i]$ is the i-th component of the FFT of $f$.
The logical consequence of this is that $\text{FFT}(f)[0]$ is real. Another way of seeing this is that the coefficient with index 0 is the "DC" component, or the "mean", and that this is of course real for a real signal.
Another consequence of course is that one needs to know only "half" the FFT to know it in full (i.e., it should be enough to know the coefficients number $0..N/2$, where N is the number of samples).
I was looking for a "real FFT" algorithm for my microcontroller and found the "official" ARM FFT. Of course they do it effectively, using an output of size N real numbers rather than N complex numbers, and outputing "only half the complex FFT" instead of the full (redundantly hermitian) FFT. This saves 50% space.
However I was surprised to see (in the discussion here https://github.com/ARM-software/CMSIS_5/issues/1091 ) that the FFT of the real signal actually has 2 purely real output values (at least, this is how they pack the output): the coefficient of index 0 (as expected) and also the coefficient of index N/2. I do not understand why the coefficient of index N/2 is real. Any hint why?
Details from the discussion there: "Except the first complex number that contains the two real numbers
X[0] and X[N/2] all the data is complex. In other words, the first complex sample
contains two real values packed."

As pointed out in comments, what I am interested in is actually some DFT properties. So that should be something like ($l$ and $k$ the DFT and function indexes, [] the "index slicing"):
$ DFT(f)[l] = \frac{1}{N} \sum_{k=0}^{N-1} f[k] e^{i 2 \pi k l / N}$


texts: Could there be any reason to prefer convolution-based calculation of autocorrelation?

texts:
Theoretically both of ways of calculating autocorrelation function are identical: strightforward convolution and Fourier-based method where we use FFT/iFFT in practice. And as it is well known, the computational complexity is $\mathcal O(N^2)$ and $\mathcal O(N\log N)$ for them respectively.

So the question is: are there any reasons why one could prefer convolution over the FFT-based method?

I could imagine only the memory-related argument (e.g. when we don't have any extra memory for calculations), but thinking of a regular PCs, at the time when memory constrains start to play role, the time complexity of $N^2$ will already be a killer for the first method.

Are there any particular numerical issue of the FFT that one needs to be aware of here (and choose convolution to avoid it)?


texts: What is the relation between input and output PSDs given system transfer function $H(s)$

texts:
If I have the system transfer function $H(s)$ in the complex frequency domain, how would I relate the input/output power spectral densities?
I have come across the relation $P_{out}(f) = |H(f)|^2P_{in}(f)$ in the frequency domain, where $P_{out/in}(f)$ refer to the input and output PSDs. Would I be able to use this same relation in the complex frequency domain as $P_{out} = |H(i\omega)|^2P_{in}$? Although I suppose that would mean the PSD would be in complex frequency domain as well?
This is all very new to me so any clarification or resources that I could look at would be greatly appreciated.


texts: Alternative to scipy maximum_filter to find peaks?

texts:
I'm using scipy maximum_filter together with binary_erosion to find peaks in my spectrum.

I just wondering if there are any other better alternatives?

I have looked at endolith's project. But the peaks from it are not spread out evenly in $x$ and $y$ axis.


texts:
Ok I need to understand how some harmonic calculations work to apply some transforms.

If I have a periodic functiion f(x). And I have the values of that function for several values of x (say x goes from 0 to 359) so I have

f(0)=3
f(1)=5
... etc

So if I'm asked to take the max of those values without the first harmonic. How can I do this? and what does it mean?

Thanks

Edit: Here is the pdf of the method I'm implementing, it's the only one I found publicaly available. There are tables of functionals on the paper which are applioed to texture analysis of images. The result of the second set of functionals is a periodic function and on the third set of functionales there are a couple of functionales where you need to get i of xi so that xi is the max or min value without the first harmonic. http://pdf.aminer.org/000/067/295/texture_classification_with_thousands_of_features.pdf


texts: Looking for an algorithm that computes "Similarity" between two signals

texts:
I'm new in this forum! ... and I'm not an expert in signal processing! ... So I hope to find help here:-)
I was wondering if someone could point me to an algorithm/technique that is used to compare time dependent signals. Ideally, this hypothetical algorithm would take in 2 signals as inputs and return a number that would be the percentage similarity between the signals.
In other terms I’m looking for a technique that allows to establish a degree of “similarity” between the two signals even if there is a delay or a shift between them. Obviously, if this technique could be implemented in a recursive way or with a small number of point ...it would be great! :-)

Thanks


texts: How to initialize an IIR filter in Matlab

texts: I designed an IIR filter and want to do zero, step and projection initialization.  There is mention in different reference papers but no explicit examples for Matlab.  How do you do IIR initialization, for the previously mentioned types, in Matlab?


texts: Image Denoising with Better Edge Preservation

texts:
I have the input image :

and the output of vein detection for the leaf using a Gabor filter, but the output is really noisy:

I tried using Total variation denoising however the results are not good:

However I don't want to loose the fine details in the leaf's veins, so a median filter won't suit my problem


texts: Which are the rules of choosing a target image for histogram matching?

texts: With histogram matching we can transform an image's intensity distribution. But this idea is very subjective. Sometimes for choosing a target image we have options no better than a guess. But is their any rule of thumb of choosing an appropriate target image for desirable improvement in input image? Is it possible to construct an artificial image and use it as a target image to obtain a specified improvement in input image? Is there any mathematical approach to obtain this? How things will change for a color image?


texts: Why Is the DFT of a Signal Symmetric About Its Central Bin?

texts: When i take an N point DFT of a signal, it comes out to be conjugate symmetric about the point N/2 . Could someone please tell how to understand this intuitively or mathematically ?


texts: Why are MFCCs of two equivalent signals completely different?

texts:
I have an audio file, of which I calculate 16 (actually 15 because I omit the first one) MFCCs in R. When I stream this file via VLC-Player and Icecast2-Server, receive it in Java (with the Player of Javazoom-Lib) and then pass it on to R again for MFCC calculation I get completely different values? Has anyone an idea why this could be?

Additional info:

When I write the received data to a file again and view it next to the original file they look almost the same (Waveform + Spectrogram).
The file has a length of 3 seconds and contains the sound of a passing vehicle (car).
MFCCs are calculated for centered 44100 samples.
I need valid data in order to automatically classify cars and trucks with a SVM.
This classifier is trained with features calculated from 150 audio files (each
containing one vehicle).
Since R had problems with reading and processing mp3-files directly I first converted wav to mp3 and back again (in order to "simulate" the loss of information due to compression).

Furthermore, I automatically detect the vehicles in a continuous audio stream and therefore (for the particular example data above) the waveforms for calculating the MFCCs have an offset of about 4400 Samples. Does this matter with an analysis window of 44100 samples?


texts: what is the differnence between 8x8 dct coefficients and 25 x 25 dct coefficients?

texts:
As I knkw the DCT is discrete cosine transform that mostly used in time domain analysis to frequency domain. Also 8x8 is used in nomaly in img processing. I don't know why that is defined by 8x8 .anyway 8x8 is used in processing.

BTW, I am want to know what is the difference between 25x25 and 8x8 in dct?
Is this just for quality?


texts: Can we detect multiple point target within a narrow-band transmit signal?

texts:
As I understand matched filter is optimal linear filter for maximizing the Signal to noise ratio (SNR) in presence of additive stochastic noise. Consider a radar/sonar system with narrow-band transmit signal and I guess generally people use matched filter for detection, specifically for range resolution. As per my understanding, under this condition we cannot resolve if there are multiple target within the signal, my first question: is it correct? We can detect multiple target within the signal if we use linear chirp with reasonable bandwidth. Multiple target within the signal refers to point targets within transmission time, (the return comes from more than one target).

My second question (sorry not a technical question) : Is anyone in the community is interested to resolve targets within transmission time  with narrow-band transmit signal in the expense of high (O(N)=$N^3$, $N$ is number of samples) computational cost?


texts: Integration with infinite/finite limits as a form of summation in DSP notation for discrete signals

texts:
While reading some DSP literature, I've noticed integrals of the following form.  An example integral (Equation 7.17) is given on pg. 124 of a book that I've found very challenging to read (Seismic Inverse Q filtering):

$u(t) = h(t)\int\limits_{ - \infty }^\infty  {u(\tau ,t)d} \tau$

This is exactly how the author has written the equation.  Gauging the author's intent, I think that $u(\tau ,t)$ is really another function, so in the above, let's re-write the equation as the following:

$u(t) = h(t)\int\limits_{ - \infty }^\infty  {p(\tau ,t)d} \tau$

In the above, $u(t) \neq p(\tau ,t)$, so I think that this clears up some of the notation.

This appears to be notation associated with a continuous signal $p(\tau, t)$ that is of infinite length.  For a discrete signal $p[\tau, t]$ with a finite length $N$, how do I evaluate the integral?

Case #1

Is the integration treated as a summation?

$u[t] = h[t]\sum\limits_{i = 1}^N {p[i,t]}$

Case #2

Is the integration treated as a "classical" integration using the trapezoidal rule?

How do I know which of these cases apply for the discrete signal case?  Is this really an issue of notation?

Another Example

Another example can be found on pg. 128 of the same book (Equation 7.21):

$\tilde U(\tau ,\omega ) = U(\tau ,\omega )\frac{1}{{\Lambda (\tau ,\omega )}}\exp \left[ {i\int\limits_0^\tau  {\left( {{{\left( {\frac{\omega }{{{\omega _h}}}} \right)}^{\gamma (\tau ')}} - 1} \right)} \omega d\tau '} \right]$

When the integration is being done over $\tau'$, is this integration a summation or a "classical" integration?

Yet another example

On pg. 9 of  "R. Wang, Introduction to Orthogonal Transforms: With Applications in Data Processing and Analysis. Cambridge University Press, 2012", the author states that the energy in a continuous signal $x(t)$ is given by:

$E = \int\limits_{ - \infty }^\infty  {{{\left| {x(t)} \right|}^2}} dt$

The energy in a discrete signal $x[n]$ is:

$E = \sum\limits_{n =  - \infty }^\infty  {{{\left| {x[n]} \right|}^2}}$

But shouldn't a discrete signal have finite lower and upper bounds of summation?  What is really meant here using infinite lower and upper bounds?  Shouldn't the trapezoidal rule be used here?

So here is the gist of my question: Using this notation, is integration done using the trapezoidal rule for discrete signals, or is it simply a form of summation?  What is the general use of this type of notation in DSP?

Given an equation for a continuous function integration, how does this relate to integration in the context of a discrete function?


texts:
If you don't have a lot of money to spend but you want to estimate the frequency response of different loudspeaker setups what is the best way to do it?

I bought a very nice microphone and built a box with foam on the inside. My idea is to place the loudspeaker and the mic in the box and measure the power level at different frequencies by playing a sine tone at those frequencies.

Or is it better to have the microphone and loudspeaker in an open space where the mic is 10 inches away from the loudspeaker?

I'm sure neither of the 2 methods above are valid for making such measurements (if so please explain why), but how should I do then?

Also, should I use sine tones or is it better to use another test signal like an "impulse", "white" noise or a chirp signal?

Hoping to get some good answers / advice

Thank you!


texts: Matlab: Could anyone suggest why Iam getting a different geometry for QPSK constellation diagram

texts:
I'm developing a Multi-user MIMO communication model in MATLAB where a base station communicates with 2 users. The model uses OFDM (QPSK) modulation. Channel model, precoding and covolutional encoder is also used. Iam sending single stream of data to each user and I'm trying to plot the receive constellation diagram for my analysis. At the receiver side, I have single antenna users. Ideally, I'm supposed to get a rectangle lattice configuration for QPSK, however am getting weird geometrical designs which am not able to understand why. I tried to check at the transmitter side how the constellation would be and it seems to be okay. I would be glad if anyone can help me out to debug this issue.
[]


texts: Why is the bandwidth for this QPSK system larger than expected?

texts:
Below is my gnuradio-companion graph.

As you can see the sample rate is 100k, I am modulating with QPSK, and am filtering the waveform (RRC) with the "Excess BW" parameter, at 0.35.  The spectrum is shown below:

The bandwidth seems to be about 2*160kHz.  I would expect it to be about the baud rate, but it is far from that.  What am I doing wrong?


texts:
I am trying to understand the difference between FFT magnitude and PSD (Power spectral density) with pwelch function of Matlab. In my example, Val1 and Val2 will output the same thing?

% N = 94144

% fs = 44100

segmentLength = N; % window size

[pxx,f] = pwelch(signal,segmentLength,[],segmentLength,fs);

%pwelch(data,segmentLength,overlappedSamples,segmentLength,fs);

Y = fft(signal);

Y_mag = abs(Y);

Val1 = 20*log10(Y_mag);

Val2 = 10*log10(pxx);

PS: I have already seen the other posts and I understand very well the difference between PSD and FFT but here I am interested in the diffeeence in code. Thanks


texts: Phase differences between DFT and FFT (in MatLab)

texts:
I have implemented a (direct) DFT in MatLab (following script) and compared it to the built-in FFT routine. The magnitude response seems to be identical (excluding some possibly round-off errors), but in the phase I get different results.

I can't find the source of these differences and I would like some insight on that.

script (only the case of even vector length is presented here):

%% Calculate needed variables
DFT = zeros(length(data), 1); % Pre-allocate result vector
n = (1:length(data))/length(data); % Calculate the n/N factor
n = -1j * n * 2 * pi; % Calculate the (-j * n * 2 * pi)/N factor

dftLength = (length(data)/2) + 1; % Calculate the length of unique Fourier coefficients

% Calculate the coefficients
for i = 1:dftLength
    exponent = exp((i - 1) * n); % Calculate the exponent
    DFT(i) = exponent * data; % Multiply the signal with the exponent

    % Add the conjugate symmetric part of the spectrum
    if i ~= 1
        % Skip first bin (it's DC)
        DFT(2 * dftLength - i) = conj(DFT(i)); % Add conjugate symmetric bin
    end
end

You can see the resulting plot here


texts: Finding filaments in high dimensional space

texts:
I have unit vectors in a high-dimensional space (~300) that are non-uniformly distributed over this hypersphere. There is physical meaning when the vectors are clustered together, but I'm hypothesizing that there is additional meaning when the vectors form a "path" along this hypersphere.

For the purposes of this question, let's assume that the deviation between the geodesic across the hypersphere and a straight line are are approximately equal (that is, I'm assuming the distances are small relative to the radius).

The question becomes, "From a large 300-dimensional point cloud in Euclidian space, how do I detect groups of points that are aligned in any direction?"


texts: Difference Between a 1st Order SG Filter And a Least Squares Moving Average

texts:
I have been studying SG filters and i recently found another filter which seem to be commonly used in financial data smoothing which is the least-squares moving average, this filter is also called moving linear regression or time series forecast, here this filter is defined as :

  The endpoint moving average (EPMA) establishes an average price by fitting a least squares straight line (see Linear Regression) through the past N days closing prices and taking the endpoint of the line (ie. the line as at the last day) as the average.

Which is similar to the SG filter in the sense that both locally fit a polynomial function using the method of least squares, however the SG filter appear smoother and more precise than the EPMA, so why are they appearing different ?

I might be wrong but my guess is that a 1st order SG filter using the coefficients show'n below :

will start convolution where N = 14 down to N = 0 while the EPMA will reverse the order of the convolution starting with N = 0 to N = 14, i don't know if i'am clear.


texts: Calculating an incoherence property from sub-optimal sampling patterns

texts:
EDIT (after comments and subject matter review)

CS is based on a choice of a sensing basis $\Phi$ relative to a representation basis $\Psi$. Using an "Incoherence Property" $\mu$ that measures the correlation between vectors from each basis it has been shown that the CS paradigm can achieve good reconstructions if the incoherence between $\Phi$ and $\Psi$ is low. Examples of low $\mu$ come from: The spike basis ($\Phi$) and the Fourier basis ($\Psi$), the noiselet basis ($\Phi$) and the wavelet basis ($\Psi$), a random matrix ($\Phi$) and a fixed basis ($\Psi$).

updated question is:

Is there a way to calculate the $\mu$ when the user is not able to chose $\Phi$ ahead of time to sample the signal, i.e. when the samples are already provided through some arbitrary sampling scheme that is sub-optimal vis-vis CS. My aim is to show/quantify, via $\mu$, that the arbitrary sampling scheme I have been given is sub-optimal if one is expecting CS to work.

ORIGINAL

With respect to Matrix Completion and Compressive Sampling (CS) I'm trying to understand how to calculate an incoherence property μ between two bases Φ and Ψ. Getting this incoherence is important because if Φ and Ψ are highly correlated there is little chance of succesfully reconstructing a signal from sparse samples. It is stated here (page 3) that μ is given by:

where n represents the number of elements in a matrix M - say an image signal.

I understand Φ to be a sensing basis and Ψ to be a sparse representation basis. I am using uniform random sampling to get a set of a samples from M. But I would also like to try arbitrary sampling patterns. (I'm not using wavelets or noiselets or Fourier coefficients - although I would like to try this eventually)

My question is:

How do I actually obtain Φ and Ψ?