texts: Translating Bessel filter from python to C - odd behavior past the cutoff frequency

texts:
I recently translated the code from scipy.signal.bessel to C in order to implement a bessel filter in a piece of analysis software I am working on. Everything went well, I was able to reproduce exactly the filter coefficients, and the filter has the right cutoff and produces numerically exact (ie, my C implementation and the direct use of scipy.signal.filtfilt agree numerically) results up until the cutoff frequency. After the cutoff frequency, though, the C implementation levels off. The image below shows the comparison of the PSD after the two filters (200kHz cutoff, 8th order, forward then back filtered) are applied to the same signal (x - frequency in $Hz$ - y - amplitude in $pA^2/Hz$). Blue is my C code, green is scipy:

My first thought was that it might be numerical noise. I use 64-bit doubles for everything in the C code, and my first guess was that maybe the python implementation uses 128-bit doubles, since the point at which the python code levels off is suspiciously close to 1/2^106. But the C code levels off well above 1/2^53, and none of my digging through scipy code suggests that they are using long doubles.

Another possiblity is that in my C code implementation, I have not implemented Gustaffson's method to ensure that the forward and backward filtered results are the same, which is done in the python version. Could that be the source of the higher noise in high frequencies?

Any ideas what is happening? I am fairly confident that the filter is implemented correctly since it gives correct results up to the cutoff frequency.

Numerically it hardly matters - leveling off at 1e-10 vs 1e-30 is irrelevent to the actual signals. I would just like to understand where the difference comes from.


texts:
I am studying H264 video these days. I know video encoder will do frame reorder for GOP with B pictures.I am not sure if my understand is correct!

GOP input to encoder (it's raw data but I just map it as GOP for better understand).

I(1) B(2) B(3) P(4) B(5) B(6) P(7) B(8) B(9) P(10) B(11) B(12) I(13)

GOP output from encoder.

I(1) P(4) B(2) B(3) P(7) B(5) B(6) P(10) B(8) B(9) I(13) B(11) B(12)

GOP input to decoder.(the same as output from encoder)

I(1) P(4) B(2) B(3) P(7) B(5) B(6) P(10) B(8) B(9) I(13) B(11) B(12)

GOP output from decode. (the same as input to encoder)

I(1) B(2) B(3) P(4) B(5) B(6) P(7) B(8) B(9) P(10) B(11) B(12) I(13)


texts: Fourier Synthesis With Varying Frequencies Over Time

texts:
I'm doing some fourier synthesis by making bins of complex numbers, filling them out with the amplitude and phase information and using IFFT to turn that into time domain samples.

I'm able to make sine, saw, triangle and square waves, and other more exotic periodic wave forms.

My question is, how do you do fourier synth for sounds that have frequency amplitudes that change over time?

Do you do it in small windows and append them together? If so, how do you make sure the phase's line up for each frequency from window to window?

Thanks so much for any info!


texts: Why does this plot correspond to this function?

texts:
I have a function:

$$
u_1(\mathbf{x})=
\begin{cases}
1, & \text{if} \;\;|x_1|\le r_1, |x_2| \le r_2 \\
0, & \text{otherwise}
\end{cases}
,\mathbf{x}=[x_1,x_2]^T \in \mathbb{R}^2
$$

And I'm told that if $r_{1}=r_{2}=5$, then plotting it in Matlab with stem3 should look like this:

Maybe I'm just too exhausted at this point, but I don't see how this is the same function. It looks like it's equal to 1 between 10 and 20 rather than between -5 and 5? Can someone please explain what's going on here?


texts: Reading Waveform from microphone - signal shape problem

texts:
I read 2kHz sinusoidal signal that I generate using youtube and mobile phone by shotgun microphone and sound card in my laptop. Results (python):

But why are the values zero at the beginning and then the values are very large and around 8 [s] the signal has a "good" shape?


texts: Array Signal Processing: What are snapshots in Direction of Arrival estimation?

texts:
I do have an overview of basic radar signal processing chain [Data cube, 3 stages of FFT, Range, Range-Doppler map, coherent/non-coherent integration, coarse DoA etc.].

I am now trying to understand an advanced DoA estimation technique-MUSIC and one of the road blocks is the snapshot concept. Nearly every paper retrieved from google assumes an prior understanding of snapshots and am yet to find one that tells me what exactly is it.

So, what are snapshots in DoA estimation?


texts:
Ok, so I'm kind of moving forward with implementing the STFT. I realise that I was taking the entire FFT and not the first half.

The problem that I'm having now is, understanding how to get the Freq bins, what calculation should I use? For example, I have been using this formula:

(10 * log10(sqrt(v[i][j].re * v[i][j].re + v[i][j].im * v[i][j].im)));

Which, I realise that this is for decibels. Currently, my graph looks like the following:

I know this makes sense to me now, this is finally starting to come together!

Would this algorithm work?:

I take the magnitude of each of the elements inside the bin `mags = (re * re + im + im) I then find the maximum index of the bin, multiply it by the Fs and divide by the bin size?

If I did this, however, would this only give me 1 result for each block?


texts: Why MATLAB fft cos makes imaginary parts?

texts:
The cos fourier transform has no imaginary parts, but in this code it has imaginary parts that little big.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7);

fs=128;

t=-4*pi:1/fs:4*pi-1/fs;

x=1*cos(0.5*t);

X=fft(x);

N=length(x);

n=0:N-1;

f=n*fs/N;

phs=angle(fftshift(X));

subplot(2,1,1); plot(f,abs(fftshift(X))/N,'LineWidth',1.5);

title('FFT'); ylabel('Magnitude'); xlabel('Frequency(Hz)');

subplot(2,1,2); plot(f,phs*180/pi,'-o'); xlabel('Frequency(Hz)');

ylabel('Phase(Degree)');

figure(8);

Re=real(X); Im=imag(X);

subplot(2,1,1); stem(f,Re); ylabel('Re(G(f))'); xlabel('Frequency(Hz)');

subplot(2,1,2); stem(f,Im); ylabel('Im(G(f))'); xlabel('Frequency(Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

If I make time range a multiple of 0.5, the results of imaginary parts has very little value that is e-14~15.

Let me know why this problem happened...


texts:
How can I understand in which time position in STFT function output am I?

For example, if I have 3987 time frame in the output of STFT, my window length is 625 (hamming), my hop size is 125 and the length of signal is 2 second.

How can I estimate for example the first 5 milisecond or a window of 10 milisecond in the ouput of STFT? I know that the the number of time frames are calculated by
coln = 1+fix((xlen-wlen)/h);
But how is the equation in my situation?
Can I say that each time frame has a length of hop size?


texts: Difference between filters in serial and parallel

texts:
I am trying to understand the difference between filters in serial and filters in parallel. As I understand it, serial is the more natural approach: A is applied then B, and the result is the product. Parallel filters confuse me though. If I think about it, I'd assume their outputs are added. Is this like a resistor in an electric circuit?
The book I'm reading (DSP by Jonathan Stein) says capicitors in series are high pass, while in parallel are low pass. Inductors the other way around. Not knowing much about electrical engineering, this confuses me.
Is the understanding I have shared at least correct? Any help with the capacitor and inductor example? Thanks!


texts:
my name is niladri, I am new to image processing(actually this is my first code). I want to implement shock filter using structure tensor. I have rough idea of what structure tensor is and implemented in MATLAB. But to design a shock I need to calculate the sign using the dominating eigenvector. But as for my understanding dominating eigenvector will be number(real and complex) for any matrix and structure tensor is basically a 2X2 matrix for each point(x,y). please correct me if my understanding is wrong. My code so far is
Edit:11/08/16: I have completed the Code but still not getting the correct response. I am unable to understand my misconception

 clc;
 clear all;
 I = imread('C:\Users\Niladri\Desktop\miramarp2gs.bmp');

 [m,n] = size(I);
 struct_scale = 3;
 %inte_scale = 5;

 x = -2*struct_scale:2*struct_scale;
 g  = exp(-0.5*(x/struct_scale).^2);
 g  = g/sum(g);
 gd = -x.*g/(struct_scale^2);
 Ix = conv2(g', conv2(gd, double(I)));

 Iy = conv2(g, conv2(gd', double(I)));

 sx = Ix^2;
 sxy = Ix*Iy;
 sy = Iy^2;

 %fileID = fopen('D:\IP\test.txt','w');
 fileID1 = fopen('D:\IP\test2.txt','w');

 for i = 1:512
    for j = 1:512
    s = [sx(i,j) sxy(i,j);sxy(i,j) sy(i,j)];
        %[v ,d] = sort(eigs(s));
        [v1, d1] = sort(diag(real(eig(s))));
        %[u , s1, d] = svd(s);
        v_perp = [0 -1;1 0 ]*v1(:,1);
        t = abs(dot(v1(:,1),v_perp));

        if v1(dot(v_perp, u(:,1)))>1E-9
            disp(t);
        end

        c  = [v1(:,1),v_perp];
        temp = det(c);

        if temp > 0
            res = 1;
        elseif temp < 0
            res = 0;
        end

        fprintf(fileID1,'%1d',res);
   end
   fprintf(fileID1,'\n');
end
fclose(fileID1);


texts: What's wrong with this implementation of my IIR filter?

texts:
I posted this originally in the electronics stack exchange, as it's for an electronic purpose (digitally controlled boost) but it was pointed out to me that this would be the more appropriate place - so here goes!

As mentioned, I recently started developing a digital regulator for a DCM boost and I've implemented the following discrete transfer function in C:

$$ G_{Z} = \frac{1.39 z^{-1} - 1.332z^{-2} }{2 - 2.002{z^{-1}} + 0.002031z^{-2}} $$

and, in terms of a difference equation:

$$ y[n] = 0.695\cdot y[n-1] - 0.666 \cdot y[n-2] + 1.001\cdot x[n-1] - 0.0010155 \cdot x[n-2] $$

My code is as follows:

// Digi Reg Coeffs
const float a1 = 1.001;
const float a2 = 0.0010155;

const float b1 = 0.6975;
const float b2 = 0.666;

typedef struct {
  float y0; // new reg value
  float y1; // last reg value
  float y2; // last^2 reg value

  float x0; // latest sample
  float x1; // last sample
  float x2; // last^2 sample
} voltage_control;

uint16_t DigitalRegulator(float error) {
  float duty = 0;

  voltage_control.x2 = voltage_control.x1;
  voltage_control.x1 = voltage_control.x0;
  voltage_control.x0 = error; /// this is the new sample

  voltage_control.y2 = voltage_control.y1;
  voltage_control.y1 = voltage_control.y0;

  // calc new y0
  voltage_control.y0 = (a1 * voltage_control.y1) - (a2 * voltage_control.y2) +
                       (b1 * voltage_control.x1) - (b2 * voltage_control.x2);

  // new value for PWM
  duty = (voltage_control.y0 / PWM_GAIN) * PERIOD_VALUE;

  return duty;
}

Is this the correct way to implement a difference equation? I'm asking because the resulting controller seems unstable, and I'd like to know if it's my implementation or something else.


texts: OpenCV/C++ connect nearby contours based on distance between them

texts:
I have to connect nearby contours in an image based on the distance between them which specifies if the contours are to be connected.

Now there's a question already on the same problem here https://stackoverflow.com/questions/8973017/opencv-c-obj-c-connect-nearby-contours but here he merges all contours into a single one. This I don't want. I don't think that there is some function in opencv for this but you can suggest an algorithm for that. My application goes like this:

I'm detecting hands, so I used a skin detection algorithm to determine them but since my skin ain't white and maybe because of lightening conditions sometime the contour breaks of at elbow. So i want the nearby contours to be connected but not all of them (because both both my hands will be there in contours.) (By hands I mean from shoulder to palm.)

Moreover I think that by using some edge detection I will get my hands boundaries and detect if some of this patch inside this boundary is detected as skin then whole region within this boundary will be detected as skin but i'm not sure how to do this part.

Any Help will be appreciated. Thanks in advance

Sample image:

In this Image I wanna connect points (8 connectivity) which are less than say 40 pixels of distance so that I'll get my left hand as a single contour

My aim is to only get hand's contour (I don't care about any other region)


texts: Quadrature Amplitude Modulation: I and Q Same

texts:
I am trying to understand a previous design that uses Amplitude Modulation.  Instead of just using the real part of the signal, the designers put the voice signal on both I and Q.  Both I and Q have the same data.

How would this affect the spectrum of the resulting AM signal and how would this affect the performance of the demodulator (assuming we are demodulating using a simple envelope detector)?


texts: Meaning of Rect and Train of Rect Spectra

texts:
The Fourier transform of $x(t)=\operatorname{rect}(t)$ is $X(f)=\operatorname{sinc}(f)$

The Fourier transform of a periodic train of rectangular pulses $x(t)=\sum\limits_{n=-\infty}^{\infty}\operatorname{rect}(t-n)$ is a sampled sinc:

Picture taken from here.
From a math point of view it is fine because:

The periodic train of rects is a periodic signal that can be expressed as a discrete sum of sine waves according to Fourier Series Expansion.

The periodic train of rects is the convolution of a rect and a train of dirac deltas. The product of their Fourier Transforms is a sampled sinc since they are a sinc and again train of dirac deltas in frequency domain.

So, my question is a about the physical meaning of this. A single rect has a continuous spectrum, that means that it has a certain amout of power (and energy) at each possible frequency (from $-\infty$ to $+\infty$). Why, if I repeat the rect infinite times, do I remove some frequencies from its spectrum and leave only a discrete set of frequencies? Why does repeating periodically a signal means removing some frequency intervals from the single pulse spectrum?
In time domain, frequency may be seen as the speed of variation of the signal. Well, it is the same for a single rect and a train of rects (if we suppose the same fall and rising times and same rect duration).


texts: Are there any difference between the FFT magnitude of a signal and average FFT magnitudes of the signal broken into multiple frames?

texts:
Suppose I calculate the FFT of a 5 seconds sound wave. If I'm only interested in the FFT magnitudes across frequencies, I should obtain the same result if I segmented the sound wave into five 1-second frames, obtain the FFT magnitude for each, and taking the average of the magnitude across the frequencies?

Calculating the FFT across a longer signal would allows for smaller frequency resolution, assuming sampling rate is constant.

Does averaging across the frames offer advantages in reducing white noise in the signal? Or does it not matter in that one would get the same result with a long FFT? When might one uses the average magnitudes across frames?

Thanks


texts: Generate time domain signal from frequency domain filter

texts:
I am familiar with using the Fourier transform to take a signal from the time domain to the frequency domain. What I would like to do is the reverse: describe a signal in the frequency domain and then take the IFT to generate the time domain signal.

Explicitly, I want to take a gaussian in the frequency domain and transform this into a pulse in time domain.  Here is the code that I have so far:

# Description of time domain sampling.
sampling_frequency = 10e6
sample_length = 100
time_index = np.arange(0, sample_length) * 1 / sampling_frequency

# Frequency domain sampling. Half the length of time domain since we only describe the positive frequency components.

frequency_domain_sample_length = sample_length // 2
sampling_frequency_domain = sampling_frequency / frequency_domain_sample_length

# Frequency domain filter.
center = 5e6
sigma = 1e6
frequency_index = np.arange(0, 10e6, bin_size_frequency_domain)
frequency_signal = np.exp(-(index - center) ** 2. / (2 * sigma ** 2.))

# Inverse fourier transform.
time_signal = np.fft.irfft(frequency_signal, sample_length)

I expect this to generate a pulse with a gaussian envelope in the time domain.

However, when I generate my signal it appears to have the "flipped" positive and negative sections that I usually correct using fftshift.  However, I was under the impression that this should not be necessary when doing an inverse real fft. Here is what I get when I plot my time and frequency signals:

My questions are:

1) Is this an appropriate method for generating time domain signals? Since I am only describing a real component of the frequency content, this is more like a "filter" or a "power spectral density" than the result one would actually get when taking the FT of a time domain signal.

2) Why is the signal split and halves reversed? How do I correct this -- fftshift or is something else required.


texts: How to calculate the noise power for a non-stationary noise?

texts:
With stationary noise we have constant mean and variance (let's assume it is Gaussian noise). My first question is, how is the noise power calculated and how it is related to the variance?

Now, I have a non-stationary noise which its variance is not constant and changes periodically over time. My second question is, how can I calculate the noise power for this noise and how it is related to the time-variant variance?

Any help and remark on the topic would be appreciated.

Thank you all.


texts: Inverse Fourier transform of complex exponential with frequency dependent shift

texts:
In the case of a constant delay $\tau$, we have the following equality:
$$\begin{align}\mathcal{F^{-1}}\left\{e^{-j\omega \tau}\right\}=\delta(t-\tau)\end{align}$$
If the delay is frequency dependent $\tau(\omega)$, can ${F^{-1}}\left\{e^{-j\omega \tau(f)}\right\}$ be expressed as a sum of diracs?


texts:
I am reading this paper Perceptual Based Adversarial Audio Attacks I have a doubt what is meant by weighting in the frequency domain and what does Section 3.4-equation-4 wants to convey .Please refer the document for detailed explanation https://arxiv.org/pdf/1906.06355.pdf

Thanks in advance


texts: Why does Fourier zero-fill interpolation pad the middle frequencies?

texts:
In this: http://dspguru.com/dsp/howtos/how-to-interpolate-in-time-domain-by-zero-padding-in-frequency-domain
the author says that we can interpolate a signal by zero-padding the middle frequencies in the DFT. But then the high frequency components are pushed even higher, so shouldn't the result look more wiggly? Why is that not so?


texts: What effect does rotation in the spatial domain has on phase in Fourier transforms?

texts:
More precisely, let's say I apply a 45 degrees rotation to an image (in the spatial domain) say, in Matlab :

Ir=imrotate(myImage,45,'crop');
FT_I=fft2(I);

In the magnitude, i.e. abs(FT_I) is it clear that the spectrum has be rotated too.

For example, trying to reconstruct/sort of "unrotate", wrt center (with shifted DC in center), in the frequency domaine like:

FT_Is=fftshift(FT_I);
I_rec_FT=abs(ifftshift(imrotate(FT_Is,45,'crop'))).*exp(-i*angle(ifftshift(imrotate(FT_Is,-45,'crop'))));

...does not make sense (ifft2()) does not show an image that seems at all sort of unrotated.

Therefore, what exactly happens with the phase part?

"unrotated" unsuccessfully, image.

From these spectra, one can clearly see that the magnitude spectrum has been rotated by the same amount. However, the phase has a random character, and one cannot visually see anything (from the phase image). I would like to know what mathematically and intuitively happens with the phase?

Moreover, it is mentioned in this lecture that rotation in spatial (or temporal) domain results in rotation in frequency domain, by the same angle, but unfortunately no precision are given w.r.t. the phase.

Moreover its stated in this lecture of Verona University and I heard it from image processing professors...

http://www.di.univr.it/documenti/OccorrenzaIns/matdid/matdid916567.pdf

https://www.slideshare.net/chinnannanperiasamy/fourier-transform-44374579


texts:
I want to use waverec to evaluate a linear combination of the scaled and shifted wavelet  of the form

$$ \sum_{i=1}^n\sum_{j=0}^{2^n-1}d_{nj} 2^{-n}\psi_{ij} (\frac{t-j}{2^n}) +c_{00}\phi(t)  $$
at $t=k2^{-n} $ for $k=0\dots2^{n}-1$. I know how do it for the Haar wavelet, namely

[c,l]=wavedec(zeros(1,8),3,'haar')
waverec(2^(3/2)*[c00 d10 d20 d21 d30 d31 d32 d33],l,'haar')

But how can you do the same thing for Daubechies wavelets? The problem being that the array c does not have 8 entries but 28:

[c,l]=wavedec(zeros(1,8),3,'db4')

c =

  Columns 1 through 14

     0     0     0     0     0     0     0     0     0     0     0     0     0     0

  Columns 15 through 28

     0     0     0     0     0     0     0     0     0     0     0     0     0     0

l =

     7     7     7     7     8

The Question is related to Number of Daubechies coefficients, but I do not understand the answer.


texts:
I have done some research on the Internet and I have found that a given 2D mask is separable if it exists only a singular value of that matrix. For example, given the following matrix:

$$A = \begin{pmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{pmatrix}$$

If I do for instance in matlab svd(A), it gives me the vector $[0\; 0\; 6]$, thus, the filter is separable. Nevertheless, I would like to know if it is possible to calculate the vectors that multiplied (or convolved) give raise to $A$. As a matter of fact, I know that:

$$\begin{pmatrix} 1  \\ 2  \\ 1  \end{pmatrix} \begin{pmatrix} 1  & 2  & 1  \end{pmatrix}= \begin{pmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{pmatrix} $$

However, I would like to calculate those vector for any separable matrix.

Thank you so much for your responses.


texts: Basic “filtering” in digital communications

texts:
I am completely stuck on solving some problems for my digital communications class. Junior level undergrad class, so fairly basic. I have been searching around all day for related articles etc but can't seem to find any that are helpful. Any help would be much appreciated. Here are the three questions:

For the coefficients my initial hunch is just to make c2 and c3 = 0, which would leave $\hat{y}[1] = y[1]$. My other attempt is to create a system of equations so the $x[2]$ terms are cancelled and the $x[1]$ terms are the only ones left over, but this system proves difficult to me and I am left with a mess of noise. For some reason I can only get two equations with three unknowns.

$$c_3= -3 c_2 \\
c_1 + \frac13  c_2 + \frac14 c_3 = 1$$
  where $c_1,c_2,c_3$ are constants.

Sorry for my terrible formatting. I believe I completed the first question correctly, here was my process.

Moving on to question 2, for part A I'm fairly sure it is what I was attempting to do before so I got c3=-3 * c2

for part b, I am stuck again because my brain keeps going back to the idea that I should just make c2 and c3 both 0, because it appears to me our goal is to make x1 recoverable.


texts: Signal bandwidth of an 802.15.4 signal at 2.4 GHz

texts: This signal uses Offset QPSK with a half sine pulse shaping filter for the purpose of  getting a constant signal envelop. Does anyone know the exact bandwidth of this signal after being pulse shaped, i.e. the nul to nul bandwidth of the main lobe?


texts: Novice question about mid-rise quantization

texts:
Accoarding to wikipage, the function for quantization has the ceiling of x/Ds.
Is this the ceiling of something multiplied by Ds?


texts: Determine filter type (HPF or LPF) by its coefficients

texts:
I want to determine optimal initial states for a FIR/IIR filter to get rid of oscillations and shifts at the initial interval.

I've found that these initial states depend on filter type in the following way (up to a constant):

for k=1:order % (Direct Form II Transposed)
   z_hp(k) = sum(1*b(1+k:end)- 0*a(1+k:end)); % for HPF:  x=const, y=0
   z_lp(k) = sum(1*b(1+k:end)- 1*a(1+k:end)); % for LPF:  x=const, y=const
end

Then I should choose z_hp or z_lp - how can I do this given only filter coefficients?

By now the solution is empirical - passing 10 points of constant value to the filter and compare the shift in output:

x = ones(1,10);
y_hp = filter(b,a,x,z_hp);
y_lp = filter(b,a,x,z_lp);

if (1-y_lp) > y_hp
    z = z_hp;
else
    z = z_lp;
end

I wonder if I there is a simple formula for this.


texts:
Is this piece of text (written by me) true about Discrete Cosine Transform ?
Discrete Cosine Transform (DCT) is a lossy data compression algorithm that is used in many compressed image and video formats,  including JPEG,  MJPEG, DV and MPEG. In this algorithm, special DCT coefficients  are calculated for each 8x8 image block.  Then the coefficients are quantized and the image block is represented as a matrix of these quantized coefficients.  The algorithm utilizes the fact that human visual system does not distinguish small changes in color or intensity.


texts: What can I do to improve quality of sound recording?

texts:
I have an issue with a sound recording I made with my Olympus LS-12.

The recorder is recording a constant unknown sound at about 10kHz. The sound is visible in the sonagram below as a faint dark band.

I uploaded a .wav file made with my Olympus LS-12 to Dropbox

Can anyone suggest what this noise is and how I can get rid of it? I've tried adjusting the gain of the recorder and recording away from my phone and computer, but this hasn't solved the problem.

I covered the microphones of the recorder with hands, to try to determine if the unknown sound is internal. This was the result:


texts: How to shift input for testing shift invariance in a system?

texts:
As I know, when we want to test if a system is shift invariant, we define system as $g(x) = H(f(x))$ which $f(x)$ is the input function, $H$ is the system and $g(x)$ is the output function. Then we use $f(x+x_0)$ as input and compare $H(f(x+x_0))$ with $g(x+x_0)$, but my problem is that why $x_0$ is not considered as input, to be clear, I give the following example from this link:
Consider the system $H(x) = f(2x)$, in the solution part, it is said that $x_0$ is just a constant and should not be considered in input and $H(f(x+x_0)) = f(2x+x_0)$, but I can't understand why? if we change the variable ($y = x+x_0$) then $H(y) = f(2y) = f(2x+2x_0)$. Can anybody please explain this for me?


texts: Can edge detection be done in the frequency domain?

texts:
Can we take advantage of the fact that high frequency components in the FFT of an image generally correspond to edges, to implement an edge detection algorithm in the fourier domain?
I did try multiplying a high pass filter with the FFT of an image. Although the resultant image sort of corresponded to edges, it wasn't exactly the edge detection established using convolution matrices.
So is there any way you could do edge detection in the fourier domain, or it's not possible at all?


texts: Is there a generically accepted method of removing motion blur from an image. For a simple case, we can assume that motion happens in a straight line. I assume that it's a two part process composed of motion estimation and then deconvolution, but how is it actually done?


texts: How to detect the presence one- and two-dimensional patterns in distributions of spatial points?

texts:
I have many point patterns, and I want to detect those patterns that have points arranged in lines that may be curved, intersected or broken. For example, on the following picture, the upper left pattern is not interesting but other three patterns are those that I want to detect.

Not all points in the interesting patterns are aligned. Some points in those patterns are still random. The lines do not have any parametric form, so I cannot use the Hough transform. I searched some literature on spatial statistics, for example this one and only found metrics for assessing
spatial uniformity and detecting clusters such as the index of dispersion. Which metrics exist for detecting the presence of not only clusters but curved lines?
I also need to do the same analysis for three-dimensional points. How to detect those 3D patterns that have points roughly arranged in lines and surfaces that can be wiggly?


texts:
I am not from DSP background and also not sure if I am posting this in right community. I have basic understanding of resistors, capacitors and filters.
My requirement is as follows-

The supply voltage is fixed and I would like to simulate the output of the filter using Python. With a fixed input voltage, I would like to know what is the output voltage at R3. The current is sampled at a frequency of 100 milliseconds. The input is a DC voltage, which can change.

I am confused where should I start from. Should I use scipy library?
Since I want a time domain output, what is the best approach?
Is there some specific reading material for such problem?
Which kind of filter would be this?

Please excuse for a basic question.


texts:
The algorithm i use to apply a discrete wavelet transform only accept raw-data with a certain number of samples (2,4,8,..128..).

I have to either interpolate the raw signal or padding it with zeros (or discard samples if size is e.g 129).

Interpolation would change the samle-rate of the entire data. Therefore I assume that zero-padding would have less effect on the "nature" of the raw signal.

How would zero-padding influence the "nature" of a raw-signal in respect to the DWT?


texts: Graphical fourier series of a square wave

texts: This is probably off-topic since it isn't really a question, but I thought that this GIF of the fourier series of a square wave was too cool not to share.


texts:
To get a panning effect I know that there is a need for a Circular Matrix.
The angle is can be set as a default or I can let the user to select the angle.

So if the angle is $\theta$, the circular matrix will be:

$$
\begin{pmatrix}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{pmatrix}
$$
and the gain is actually set in which channel (Right or Left) the audio will be.

So if I'm sampling the signal from a guitar and get the sequence $X$,
all I need to do is to do is to multiply each sample value with the circular matrix?

This operation (multiply each sampled value with the circular matrix) calculates the gain right? So I will know in which channel - right or left to play the signal.

I don't understand when the signal is converted again to analog. After the calculation of the gain? Because if you play with a panning effect, you will hear the sound in the left or the right channel very fast.

I don't understand if it's a multiplication the whole signal with the circular matrix and then just playing it back on the different channels. Or is it actually a multiplication of each sampled value of the signal by the matrix, conversion to analog and then processing of  other sampled values.


texts: Can I compute convolution through the overlap and add method without the IFFT?

texts:
So, I'm (I think) aware of how the overlap and add algorithm for linear convolution works, but my question is that, suppose I have a FFT-ed set of sequences that belong to a large sequence.

Can I add the overlapping parts before taking the IFFT of each sequence?
Can I add the overlapping parts and then perform another convolution without taking another set of IFFTs and FFTs?
Let's say that I have a function which multiplies by 3 if the number is 0 in the time domain, but multiplies by 6 if it isn't. I can perform this with only one IFFT as opposed to an IFFT+ a FFT by taking the IFFT and then multiplying in the frequency domain by the required number. Can I do a convolution through overlap and add, apply the aforementioned function through the "one IFFT only technique" and then do another convolution through the overlap and add method without having to do another set of FFTs for the next convolution through the overlap and add?


texts: Is HSV Color Space Sufficient for Rudimentary Color ID and Edge Detection

texts:
I want to, as part of a C++ school project (no OpenCV), work with an RGB image taken from a camera filming smarties (round colored pellets) from above and extract the number of candies in each color.

Is converting the (reasonably lit and clear) image to HSV sufficient to perform edge detection (with Sobel on the V channel) and color identification (with H and S), following noise removal and possibly histogram (or other kinds of) equalization? Or should I stick with RGB color space and do more thorough filtering in order to properly isolate the same color of differently lit smarties.

Hope this makes sense, thank you for your time


texts: Evaluation of phase difference between two acquired signals (expected to be sine waves) with the same frequency

texts:
I would like to calculate the phase difference between two signals $F(t)$ and $R(t)$ which are expected to be:
\begin{align}
F(t)&=A_F\cos\left(\omega t+\phi_F\right)\\
R(t)&=A_R\cos\left(\omega t+\phi_R\right)
\end{align} What I do at the moment is to acquire these two signals, multiply them and then by averaging I can get the phase difference as follow:

$$\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$$

Now this works fine, but I don't like the fact the I can represent the phase difference only between $0^{\circ}$ and $180^{\circ}$.

Is there an alternative method that I can use to evalutae the phase difference on the full circle $0^{\circ}$ to $360^{\circ}$?

I already tried using FFT and calculate a Cross Power Spectrum, but unfortunately this method is too slow compared to the one mentioned above. I use LabVIEW for the whole acquisition and analysis. The calculation time is crucial for me and I would like to keep it minimal.

EDIT: Here some additional information about the application. I acquire the two $F(t)$ and $R(t)$ with a sampling frequency of 10 MHz. $F(t)$ is a signal (sinusoidal) generated by a signal generator and it drives a piezoelectric element. $R(t)$ is the response of the device that I drive with the piezo. Both signals are acquired simultaneously and I'm interested in the phase difference between the two signals. After acquiring them I follow the following algorithm (in LabVIEW):

Estimation of $A_F$ and $A_R$ from their RMS values.
Multiplication  $X(t)=F(t)R(t)$
Calculation of $2\textrm{Mean}[F(t)R(t)]$ (twice the DC level of $X(t)$)
Calculation of $\cos({\phi})=\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R}$
Calculation of inverse cosine $\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$

I then repeat the algorithm for a list of different driving frequencies to build a phase spectrum. The phase spectrum looks as expected (Simple Harmonic Oscillator) but it wraps around $180^{\circ}$. The problem is that it doesn't wraps by restarting from $0^{\circ}$ but once it reaches $180^{\circ}$ it goes backward towards $0^{\circ}$. This means that the phase spectrum is continuous but it presents something like a cusp around $180^{\circ}$. I tried some unwrapping, like flipping the portion of the phase spectrum that should be after $180^{\circ}$ but the signals have some noise and they are not perfect sinusoidal signals. This means that I cannot flip by $180^{\circ}$ because the spectrum then doesn't look continuous but with a small step. Here some screenshots of the unwrapping techniques I tried so far:

a. Calculated phase spectrum:

b. Unwrapping around $180^{\circ}$:

c. ZOOM around $180^{\circ}$:

As you can see if I unwrap around $180^{\circ}$ I get the small step. This is a big problem for me because I then fit a simple harmonic oscillator model to my phase spectrum and with that step the fit is not optimal.

EDIT 2: Some additional parameters:

a. Typical frequency ranges of $F(t)$ are $280-500\ \textrm{kHz}$ in particular the frequency range of the images above is $370-420\ \textrm{kHz}$.

b. The frequency around $180^\circ$ is $403\ \textrm{kHz}$.

c. The two signals are acquired at $10\ \textrm{MHz}$ for $1\ \textrm{ms}$ which means there are $10^{3}$ data points for each signal.

d. The amplitude of the reference signal $F(t)$ is $\pm 2\ \textrm{V}$.

EDIT 3: $F(t)$ is generated by a function generator and acquired by a digitizer which acquires simultaneously $F(t)$ and $R(t)$. The letters were chosen this way because for me $F(t)$ is the driving Force of my mechanical system and $R(t)$ is the Response of the mechanical system. The mechanical system is an oscillating cantilever, which is a beam fixed at one end. The beam has a rectangular cross section with length, width and thickness of $500\ \mu \textrm{m}$, $100\ \mu \textrm{m}$ and $1\ \mu \textrm{m}$, respectively.


texts: DFT normalization for amplitude estimation

texts:
After an interesting recent answer, I'm doing some research on proper DFT normalization for sinusoidal peak estimation. It's clear that to get the correct amplitude from DFT peaks we need to normalize by the window gain
\begin{equation}
S_w \triangleq \sum_{n=0}^{N-1} w[k]
\end{equation}
Which becomes $S_w = N$ for the rect window.

Now I'd like to find some rough demonstration, Harris has one (see equations 12 and 13 in the article below), but he chooses the easiest path by making exponentials cancel out. Well he actually does a pretty good job, it's me that thought it could be inferred to some more general relation.

For the rectangular window I've come up with this.
Let's say we have a sampled analytic sinusoidal signal $x[n]$, with frequency of $m/N$ cycles per sample (or $m$ cycles in $N$ samples), and compute the DFT $X[k]$

\begin{equation}
x[n] = A \ e^{j 2\pi m n/N}
\end{equation}

\begin{align}
X[k] & = A\sum_{n=0}^{N-1}e^{j \frac{2\pi n}{N} m} e^{-j \frac{2\pi
  n}{N} k} \\
\\
& =  A\sum_{n=0}^{N-1}e^{-j \frac{2\pi n}{N} (k-m)}\notag \\
\\
& = A \ \dfrac{1-e^{-j2\pi(k-m)}}{1-e^{-j\frac{2\pi}{N}(k-m)}} \\
\\
& = AN\delta \left[ k-m\right]
\end{align}

where $\delta[k]$ is the Kronecker delta.
To get the correct amplitude in freq space we need to divide by $N$.

Now, let's say we have a windowed signal, let me call it $x[n]$ again

\begin{equation}
x[n] = w[n] \ A \ e^{j 2\pi m n/N}
\end{equation}

and its transform, again overloading notation

\begin{equation}
X[k] = A\sum_{n=0}^{N-1} w[n] e^{-j \frac{2\pi n}{N} (k-m)}
\end{equation}

From here I'm not sure how to proceed, guess the result should be
$$
\require{cancel}\cancel{X[k] \overset{?}{=} A \ S_w \ \delta[k-m]}
$$

Edit: after some discussion with Robert I believe the previous equation is just wrong. It works for the rectangular window only, while for the other windows the only thing we can say is about the peak $k=m$, which can be easily demonstrated just from the sum.

  G. Heinzel et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform


texts: How do I deconstruct a signal with fast fourier transformation in Python?

texts:
I want to use the fast Fourier transformation in python to deconstruct a periodic signal into its harmonic periodicities.

Since the FFT of the original signal is the same as adding the FFT of the periodicities, I was calculating the power spectrum of the FFT of the original signal and then tried to extract the several signals. You can see my code for the extraction below:

def harmonic_period_calc(old_fft, peak_index, peak_value):
"""
    description: this function recalculates the fft so that only one peak in
                the psd should be visible

                 furthermore it uses the ifft to calculate a harmonic signal,
                 from which you are able to calculate the time indexes via
                 its peaks

    old_fft: the old fourier tansformation of our signal
    peak_index: the index, of peaks in the power spectrum
    peak_value: scalar-factor, with which we are going to shorten the other
            values with

"""

    new_fft_one = old_fft[:int(len(old_fft)/2)]
    new_fft_two = old_fft[int(len(old_fft)/2):]

    new_fft_one[np.delete(np.arange(new_fft_one.shape[0]), peak_index)] *= 1/peak_value
    new_fft_two[np.delete(np.arange(new_fft_two.shape[0]), -peak_index)] *= 1/peak_value

    new_fft = np.concatenate((new_fft_one, new_fft_two))

    period = scipy.fftpack.ifft(new_fft)

    return period

Unfortunately this is only working okay. Because the peaks of the inverse FFT is slightly off from the original period.

So my question is, does anybody know, how I will be able to do this?

If necessary, this is how I calculate the FFT and the power spectrum:

temp_fft = scipy.fftpack.fft(signal)
temp_psd = np.abs(temp_fft) ** 2
fftfreq = scipy.fftpack.fftfreq(len(temp_psd), 1/len(signal))
i = fftfreq > 0

Also the technical information from me:

Python 3.7
scipy 1.3.1
pandas 0.25.3
numpy 1.17.2


texts: feature selection for classification for rare classes

texts:
I have a spectrogram of a signal and want to classify it. I used different features such as MFFC, the density of power spectrum and other common features and classify using the random forest, but I did not get good results.

I think that my features are not well selected and I guess that there is no relation between my aimed classes and the selected features.
my classes are so rare.

could anyone suggest a way for obtaining features?(for very rare data)
also, I used some methods of feature extraction from biological images, but I think that the data are different. my data are double and those are some 3 dimensions data.


texts: Determine reflections from received signal

texts:
I have a reference signal $r(t)$ and the correlation between that reference signal and the received signal : $C_{XR}(\tau)$. The signal I receive contains reflections on walls. I have to build a synthetic signal $s(t)$ with the addition of attenuated and time-shifted versions of the reference signal. The goal is to build the synthetic signal in a way that it is as close to the received signal as possible. If I do that, I get informations about the walls on which the signal has been reflected.

My idea is to write the synthetic signal as follow : $s(t) = \sum_{i=0}^{N} A_i r(t-T_i)$ where $N$ is the number of reflections we consider.

Then, the correlation between $s(t)$ and $r(t)$ is :

$$C_{RS}(\tau) = \sum_i A_i \ IFFT(R(\omega) R^*(\omega) e^{j\omega T_i}) = \sum_i A_i C_{R}(\tau+T_i) \simeq C_{XR}(\tau)$$

After that, I don't know how to determine $A_i$ and $T_i$ in order to get the last equality.


texts: Detecting a frequency swept sinusoid and its parameters?

texts:
Given a (FFT-sized) frame of data, and detection of a spectral component statistically above the noise floor in the FFT of this window, what characteristics or signal analysis could be used to determine that this spectral component is more likely to be a linearly swept sinusoid, rather than one that is stationary across the frame?

And, assuming the dF/dt sweep across the data window is small (from a fraction of an FFT bin to a couple bins), how can one estimate the sweep parameter (but, beyond answer offered to this question, assuming this is an estimation of a detected signal in noise.)

One offered solution seems to be to segment the FFT frame into several shorter subframes, do shorter STFTs, and look for a linear best fit among the resulting set of subframe FFT peak magnitude frequency estimates (which all have poorer frequency resolution due to the shorter subframes).  Are there any other or better options for detection and estimation?


texts: Inverse Fourier transform Of a triangular impulse

texts: I have to find the expression of this graphic and after find the inverse Fourier transform of it. First of all I found that the expression of the graphic is $$ X(f) = \frac{1}{2} tri (\frac{f+f_0}{B}) - \frac{1}{2} tri(\frac{f-f_0}{B})$$. Now to find inverse Fourier transform , my book give me the advice to multiply numerator and denominator for i. But $$ i=e^{i \frac{\pi}{2}} $$ and $$ -i=e^{-i \frac{\pi}{2}} $$. Now I know that the Fourier transform of a triangular impulse is $$ (sinc(f)^{2}) $$ and that $$ \frac{d}{dt} tri(t) = rect ( t + \frac{1}{2}) - rect ( t - \frac{1}{2}) $$ but I don’t know how to apply correctly integration property of my x(t). Thank you so much


texts: Convolution of $\sin(\omega t)$ and $\cos(\omega t)$?

texts:
If $x(t)=\sin(\frac{\pi t}{4})$ and $y(t)=\cos(\frac{\pi t}{4})$ then i need to find the Convolution $$z(t)=x(t) \circledast y(t)$$

So convolution will be

$$\begin{align}
z(t) &=\int_{-\infty}^{\infty}x(\tau)y(t-\tau)d\tau \\
\\
 &=\int_{-\infty}^{\infty}\sin\left(\frac{\pi \tau}{4}\right)\cos\left(\frac{\pi (t-\tau)}{4}\right)d\tau \\
\\
 &= \tfrac12 \cos\left(\frac{\pi t}{4}\right)\int_{-\infty}^{\infty}\sin\left(\frac{\pi \tau}{2}\right)d\tau+\sin\left(\frac{\pi \tau}{4}\right)\int_{-\infty}^{\infty}\sin^2\left(\frac{\pi t}{4}\right)d\tau \\
\end{align}$$

now first integral will come $0$ because odd signal and 2nd integral value will come $\infty$ that is undefined.That means convolution is not possible for these signals.

But if i do it by fourier transform then i will get
$$\begin{align}
Z(j\omega) &= X(j\omega)Y(j\omega) \\
 &= \frac{\pi^2}{j}( \delta(\omega- \tfrac{\pi}{4})-\delta(\omega+\tfrac{\pi}{4}) ) \\
 &= 2\pi^2\sin \left(\frac{\pi t}{4} \right) \\
\end{align}$$

How its possible that i can get it by fourier transform but not by direct convolution?

However the answer given is $4\sin(\frac{\pi t}{4})$.

What is the mistake i am doing please help me?


texts: Periodic signals in Continuous and discrete time

texts:
Is there any signal which is periodic in Continuous Time but not in Discrete Time?
I have this doubt prevailing in me for a long time.
Are all CT periodic signals periodic in DT?
If so, how is sin wave a periodic signal?
It is periodic only at intervals of pi radians, which is an irrational number and not an integer.
But for a DT signal to be periodic shouldn't the signal be periodic at discrete intervals of time?


texts: How can I extract information about timbre from a .wav file. In which form can I study those features?

texts: To study the frequencies of a piece of audio, we can generate a spectrogram from that audio file. I'm wondering how the features of timbre are saved in such tiny wave files, and how we can extract those features and study the timbre? Thus we know what instruments are played in the audio?


texts:
I'm new to audio processing and I hope here be a good place to ask about the definition of terms in this context.

My question is what is difference between sound, audio, voice, speech and tone, from both technical and mathematical point of view?


texts: How to calculate particles sizes for spheric particles with overlapping and superposition? (Example image included)

texts:
I'm trying to automatically calculate particles sizes histogram for images like this one:

Common threshold and segmentation doesn't seem to work. What filters and algorithms do you suggest to be used?


texts:
DSP-Newbie here... Why are we using the ratio of the center frequency of the signal to the sampling frequency of the signal spectrum when we downconvert the signal?

I have an Spectrum with Carriers which reaches from $-f_{s}/2 ... f_{s}/2$ ($f_{s}$ is the Sampling-frequency) (see picture)

The downconversion/shifting of the spectrum is

$$e^{-j2\pi(LO/f_s)}$$


texts: Cutoff frequency of a first order recursive filter

texts:
I've implemented a very simple first order recursive low pass filter in c using the algorithm:

double a, b, prevOutput;

double lowPassFilter(double input)
{
    double output = a * input + b * prevOutput;
    prevOutput = output;
    return output;
}

where:

a = 1.0 - b

According to The Scientist and Engineer's Guide to
Digital Signal Processing, the cutoff frequency of the filter is defined by the relation:

  $$b = e^{-2πf_c}$$

$f_c$ is the -3db cutoff frequency as a fraction of the sampling frequency. Using $f_c = 0.25$, I calculated that $b = 0.20788$ (approximately) and plugged it into my filter. I wanted to verify my filter, so I passed an impulse function through it and recorded the output:

double result[128];
for (int i = 0; i < 128; ++i)
{
    if (0 == i) result[i] = lowPassFilter(1.0);
    else result[i] = lowPassFilter(0.0);
}

I then performed FFT on the result to find the frequency response:

complex_t fftResult* = fft(result);
for (int i = 0; i < 65; ++i)
{
    double amplitude = fftResult[i].real * fftResult[i].real;
    amplitude += fftResult[i].img * fftResult[i].img;
    amplitude = sqrt(magnitude);
    printf("f = %.3f: %f\n", i / 128.0, amplitude );
}

The result I saw was somewhat expected but also not exactly what I expected. In the DC bin ($f = 0$) I saw a amplitude very close to one and I saw the amplitude attenuate as it approached 0.5. However, I also expected that at $f = 0.25$ the amplitude should have been around $0.70795$ (i.e. -3db). but instead the amplitude was around $0.775$:

 ...
f = 0.242: 0.782501
f = 0.250: 0.774971
f = 0.258: 0.767380
 ...

Aside from the inefficient code, is there something wrong with the method I am using to obtain the frequency response or am I interpreting the cut-off frequency relation incorrectly or something else?


texts: Combining two filters (audio processing)

texts:
I'm trying to combine two filters into a single filter command for a homework problem and was looking for advice. I don't necessarily need an explicit answer as long as someone can point me in the right direction.

In the following code, assume that 'y' is a sound file that has been pre-loaded (10kHz). Also, both the hlpf and hnotch filters work together separately to produce an audio signal with noise and a beeping tone removed. However, when I try to combine the two into a single filter, I do not achieve this result. Here is my code:

Fs = 10000;
f0 = 1000;                %#notch frequency
fn = Fs/2;              %#Nyquist frequency
freqRatio = f0/fn;      %#ratio of notch freq. to Nyquist freq.
notchWidth = 0.1;       %#width of the notch

%Compute zeros
zeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];

%Compute poles
poles = (1-notchWidth) * zeros;
b = poly( zeros ); %# Get moving average filter coefficients
a = poly( poles ); %# Get autoregressive filter coefficients

hnotch = [poly(poles)];
alpha = 0.42;
hlpf = fir1(100, alpha);
hcombo = conv(hlpf,hnotch);
s2 = filter(hcombo,1,y);
soundsc(s2,Fs)

Any help would be greatly appreciated.


texts: how to compute a discrete fourier transform on fragmented data

texts:
I have a list of collected data points that I need to take a DFT of, however with the problem that about one quarter of the data points in the middle are missing, and so even though the existing data points are evenly spaced their is still a large time gap that needs to be accounted for. I looked up some things on a variable sampling rate DFT, but it was first of all over my head and secondly I am not sure if it is applicable. My question then is how would I account for that time gap in my calculation. It would be best to have everything explained in as close to laymen's terms as possible as I am out of my depth in this.

Thank you for any help.


texts: Why is the noise on an eye diagram not in the centre i.e. at 0?

texts: I don't yet completely understand eye diagrams. Why isn't the noise in the centre of this polar BPSK diagram and instead are 2 margins at -4 and 4.. surely that's not where the noise would be in relation to the carrier. Also, how do you work out the SNR looking at this, I'm assuming it's the proportion of the full height of the eye to the noise height?


texts: What is the algorithm to generate sine waves of arbitrary frequency in the STFT domain?

texts:
I'd like to write a DSP algorithm to do additive synthesis using arbitrary sine waves with inverse rectangular FFTs.
This requires two things:

The ability to generate phase/amplitude lists that cause the IFFT to create sine waves that are not integer multiples of 1/(IFFT length)

Example:  How do I generate 1000 Hz with 16 sample IFFTs at 44100 Hz?

The ability combine two frames pre-IFFT so they are added post-fft.

What are the necessary steps here?


texts: How do you calculate RMSE in degrees for a DOA estimator system with multiple incoming signals?

texts:
Let's say you have a system with 1 transmitter and 5 receivers. We'll use the MUSIC algorithm to determine DOA of incoming signals.

If we have 1 signal (at 0° for example), it's easy to calculate RMSE vs SNR. Simply take the peak of the MUSIC spectrum (an angle) and find the difference between the actual target angle.

However, what if we have 2 signals? One at 0° and a second at 10° for example. You can no longer look at the MUSIC spectrum and simply pick the 2 greatest points, since the absolute greatest points may be at 0.0° and 0.1°, even if there's a clear (lower) peak at 10°.

What would be the/an appropriate way to measure the error when multiple signals are present? Or perhaps another metric (besides RMSE) is needed, such as probability of observing 2 peaks, etc.


texts: Is there a good metric to compare character images by similarity?

texts:
Good evening,
I'm thinking about a good way to compare images of characters with a metric.
Infos:
Image colors - black/white; one may think of it as a matrix with a value from 0 to 255 for each pixel
Image structure - contains in most cases 1, but it may contain some more characters (up to 5)
Image sizes - unluckily very different (i. e. matrices have different sizes; between around 3x3 up to 90x25)
I've looked into these two articles, but they don't give me any idea.
Are there any good "conceptual distances" for comparing two images?
What distance metric can I use for comparing images?


texts: How to Find the Probability of Kalman Filter States?

texts:
I am working on a video object tracking problem. I am using Kalman filter to predict and correct the object position return by an algorithm such as CamShift. I want to adjust the likelihood probability of the states so as to customize the weight of the states before feeding to the Kalman filter. Currently this the code I have:

#Define kalman filter
kalman = cv2.KalmanFilter(4, 2, 2)
# initial position: (x0,y0)
state = np.array([col+w/2,row+h/2,0,0],dtype='float64')
kalman.transitionMatrix = np.array([[1.,0., .1, 0.],
                                [0., 1., 0., .1],
                                [0., 0., 1., 0.],
                                [0., 0., 0., 1.]])
kalman.measurementMatrix = 1. * np.eye(2, 4)
kalman.processNoiseCov = 1e-5 * np.eye(4, 4)
kalman.measurementNoiseCov = 1e-3 * np.eye(2, 2)
kalman.errorCovPost = 1e-1 * np.eye(4, 4)
kalman.statePost = state

Here, 'row' and 'col' are the top-left row and column indices, respectively of the current rectangular box containing the object position, and 'w' and 'h' are the width and height, respectively of the rectangular box containing the object.

Now, suppose (x1,y1) is the position returned by my algorithm (e.g. CamShift,Optical flow etc.).

Now I pass (x1,y1) and the explicitly calculated velocity (v) and acceleration(a) to Kalman filter as:

#Predict using kalman filter
prediction = kalman.predict()

# Update kalman state using optical flow results
pts=np.array([np.float32(x1),np.float32(y1),np.float32(v),np.float32(a)],
dtype=np.float64)

state=pts
measurement = (np.dot(kalman.measurementNoiseCov, np.random.randn(2,
1))).reshape(-1)
measurement = np.dot(kalman.measurementMatrix, state) + measurement
# use optical flow to correct kalman filter
prediction_new=kalman.correct(measurement)

It is all fine till here. But what if I want to give only less than 0.5 probability weightage to the previous state (=[x0,y0,0,0]) compared to the present state (=[x1,y1,v,a]) while computing 'prediction_new'. I have extensively searched the Kalman filter OpenCV library as well as online, but couldn't find an answer. Appreciate help on this. Thanks in advance.


texts: Why does spectral leakage arise in an FFT?

texts:
I've tried googling and wikipedia-ing it, but I haven't gotten any answers beyond 'it's because the frequency of the input signal is sitting between two bins'.

I understand that this is the reason, but what I can't understand is why the leakage seems to extend to several adjacent bins rather than only one adjacent bin.

To illustrate what I'm talking about, here's some simulated data (code at the end of the post) :

Above is the FFT spectrum (plotted on a log scale) of a sine wave of frequency 10. The sampling rate is one, and the number of samples is 100. The graph has been FFT-shifted. There's clearly only a peak at bin 10, and the rest is on the order of numerical error, or there about.

This is the frequency spectrum at a generated frequency of 10.1.
Clearly there is 'leakage' into more bins than just the immediately adjacent bin.

This is the plot for a frequency of 10.5.

Question: Why is there this leakage, and why does it extend to all the other bins, rather than the immediate adjacent bin?

Code, for anyone who's interested (Python code)

import numpy as np
import matplotlib.pyplot as plt

xFreq = 10.5
xSize = 100.0
xPeriod = xSize/xFreq
x = np.linspace(1,xSize,xSize)

data = np.sin(2*np.pi*x/xPeriod)
fft = np.fft.fft(data)
fft = np.fft.fftshift(fft)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(abs(fft), "o")
ax.set_yscale('log')
plt.show()

I changed the xFreq value from 10.0 to 10.5, etc.


texts: Binary classifier for a small range of the audio spectrum

texts:
I am trying to detect a pitch in a narrow range of the audio spectrum with minimal samples. This corresponds to the rattling of certain mechanical systems. For example, the hum of an engine.

To do this, I am hoping to identify if there is a frequency at at a certain range such as 550hz-555hz.

I spoke with a learned fellow about this problem and he mentioned that there is a variation of the FFT that only targets the desired frequency range. What is this called?
I was hoping to avoid an amplitude (volume) training period by distinguishing peaks from white-noise in the frequency domain. Can anybody point me to a tunable parameter that will enable me to identify if a certain range has a peak? I was thinking of comparing the max height to the average.
I was wondering if anybody knew of a good method to do this.


texts: Understanding time delay estimation for GPS?

texts:
I really need help to understand these questions which is highlight, for the GPS situation, I hope anyone who has been working in this situation before, please help me to have a clear understanding about "time delay estimation":

Here is the scenario:

Suppose I know precisely I have transmit the signal (data), I know the sampling, I know the time correspond to the very 1st point, I can go through the sign the time every single measurement, that I just transmitted. Now, I received the signal (data), I have section of data and in there somewhere is the signal that transmitted. Unlike, as the transmit, I know the time correspond to the very 1st point in my receive data:

what would I do to get the time delay?

what steps do we follow?

In there, somewhere, I don’t know where exactly, but in there somewhere is transmitted the signal and I want to estimate time delay between the transmit and receive. So, I want to look at:

what you need to do that?

What math. operation do I need to use to compute that time delay?

If I say:"I want to compute xcorr between 2 vectors". So, this xcorr is twice sum of the length – 1:

what do I do when I have xcorr output?

I take 2 vectors, and I compute even longer vector use xcorr, so now I compute the longer vector but I want a number, I don’t want a vector, I want a time delay estimation:

what do I do with the xcorr solved like a MATLAB?

If I say xcorr(x,y), x has certain length and y has certain length, the result is the length of x plus length of y and minus 1. So, now I have a long vector:

what do I do then to determine the time delay?

what the signal look like and what xcorr look like?

what do I do when I have xcorr to determine the time delay?

6.I know that have xcorr with bunch of complex number but I need only one real number. The real number that I want is the time delay. I take input data, I take received signal, I compute xcorr, then I do something to help myself to get the time delay. Result to that is single real number. But these are general ideas, I need to answer those question.

I have read this link: Time delay estimation of oscilloscope signals using cross correlation . It's really help my understanding, but they way he answered very high level, I think I need to understand the basic first, then continue with this link.

Please help me for these questions.

Thank you very much, everyone.


texts: How to process video with DFT/FFT - getting the frequency of green color changes?

texts:
Goal: I have a video stream from a webcam. I need get frequency of the green color changes on the video.

Video color space is RGB
The signal looks sinusoidal with a frequency range 0.5 - 2.4 Hz

I performed the following steps:

Capture Video frames into collection.
Processing frames by Gaussian pyramid down level 4. (frame size now is ~ 37x32)
use DFT/FFT with sample rate 30
get spectrum
calculate frequency

on DFT I put data from each frame:
video frames quantity = 60

Method 1:
I collect green channel from frame and join with the green channel from a next frame, in result I have a big data array.

37*32=1184

1184*60=71040 the size of the sample buffer

one long array with 71040 data point of all pixels in the green channel of each frame.
(all green pixels from one frame + from next frame)

Method 2:
I collect the green channel from an each pixel on the frame and join with the same green channel pixel from a next frame, as a result I have a big data array.

1184 arrays with 60 data points
(each green pixel from one frame + from next frame)

Which method is the correct one?


texts: Band limited signal: Infinite time duration, time limited signal : infinite bandwidth

texts:
In signal processing, there is a time-frequency duality. Band-limited signals have infinite time duration and time-limited signals have infinite bandwidth.

Where can I find a proof of this and is there some intuitive explanation for this?

I am not sure if it has to do something with fouriers uncertainty principle?


texts: Relationship between Wavelet transform and Fourier Power Spectral Density

texts:
Is there anyway to obtain the Fourier Power Spectral Density from a wavelet transform of a time series?

I am particularly interested in this problem because I was wondering if there is any possibility to obtain the local Power Spectral Density from the wavelet transform.

If I am not wrong, according to Torrence and Compo, the average of all the local wavelet spectra tends to approach the Fourier Spectrum of the time series.

However when I compute the Wavelet spectra the results is much widther than the one given by the Fourier Transform.

Here there is an example:

I made an oscillatory time series with some noise (first panel)
I compute the Wavelet Continious transform and then I compute the wavelet power spectra (second panel) using a Morlet Transform.
According to Torrence and Compo, we can compute the Fourier periods with the formulae than they provide, and also we can compute the Fourier frequencies. They said that a vertical cut in any time is the local Fourier Spectrum, and if we average all the local wavelet spectra we will obtain the Fourier Spectra.
However when I do this the result obtained from the wavelet is much wider than the Power Spectra obtained from Fourier:

The main frequencies seems to be OK, but I don't understand why the spectrum from the wavelets is so wide.

To rephrase the question again, I would like to know if I can approximate the Fourier Spectra with the wavelet transform without this wideness.
And again: I am particularly interested in this problem because I was wondering if there is any possibility to obtain the local Power Spectral Density from the wavelet transform.

Thank you very much.


texts: Problems discretizing a notch filter - Chebyshev Type I

texts:
I'm trying to implement and discretize a notch filter (Chebyshev type I) these specifications:
Ws = 2pi[2000 2500];
Wp = 2pi[1000 7000];
Rp = 0.5;
Rs = 40;
But when I try to use discretization methods, Tustin has a weird phase graph, and ZOH has also a weird behaviour. I couldn't find anything online that would explain those behaviours. Could anyone help me?

Discretization matlab code:
% continuos to discrete
Ws = 2*pi*[2000 2500];
Wp = 2*pi*[1000 7000];
Rp = 0.5;
Rs = 40;
[n,Wp] = cheb1ord(Wp,Ws,Rp,Rs,'s');
[b,a] = cheby1(n,Rp,Wp,'stop','s');
H = tf(b,a);

Hbi = c2d(H,1/20000,'tustin');
Hzoh = c2d(H,1/20000,'zoh');
Hm = c2d(H,1/20000,'matched');
Hdp = c2d(H,1/20000,'prewarp',2*pi*2250);
%Hii = c2d(H,1/20000,'impulse');

options = bodeoptions;
options.FreqUnits = 'Hz'; % or 'rad/second', 'rpm', etc.

bode(H,'b',Hzoh,'r',Hbi,'m',Hm,'g', Hdp, 'y',{2*pi*100 2*pi*10000}, options), grid
legend('Continuous','ZOH','Tustin','Matched','Tustin w/ pre-warping')

Thanks in advance!


texts: How to insert \ell in text box on figure in matlab?

texts: I have a plot in matlab and want to insert a text box write inside ell symbol. Till now i couldn't do it. I wrote '\ell','interpreter','latex' but it doesn't accept. I appreciate If anyone could help. Thanks.


texts:
I've succeeded to compute complex FFT of a data array of interleaved values using the arm_cfft_f32() API of the CMSIS DSP library for Cortex-M4, and got in the same array the FFT results as frequency domain values with the same interleaving.

I chose a result value of the Frequency Domain (let say a point) based on a specific criteria, and want to go back to the initial index (in the Time Domain) before calculating the FFT and passing through arm_cfft_radix8by4_f32() API.

Is there anyway to go in the opposite sens, and identify a point based on its FFT result value ?

Thanks and regards,
tinymanta.


texts: Given an image how should I detect this is a day picture or night picture?

texts:
I have a 2D image which contains R G and B values. I think this should be purely based on averaging of pixels and selecting the threshold to decided whether this is day picture or night one.

Can you provide any kind of different ideas to solve this problem?

I am shifting this question from the https://stackoverflow.com/questions/28157035/given-an-image-how-should-i-detect-this-is-a-day-picture-or-night-picture?noredirect=1#comment44685520_28157035


texts: Dimensional analysis of integrated white noise process

texts:
This question is somewhat related to this post.

Let us consider we have a white noise current source $i_n(t)$, with a variance $\sigma_i^2$, and mean, $\mu_n=0$.

Assume that this current is passed through a system, which integrates for a finite time duration ($T_{int}$) on a capacitor $C$. One can then write an equation for the power of the noise voltage (integration of current produces voltage $V=\frac{1}{C}\int_{0}^{T}i(t)dt $  ) produced from this integration process ($\sigma_v^2$):

$$\begin{align}\sigma_{v}^{2} &=E\bigg(\frac{1}{C}\int_{0}^{T_{int}}i_{n}(t_{1})dt_{1}\frac{1}{C}\int_{0}^{T_{int}}i_{n}(t_{2})dt_{2}\bigg)
\\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\int_{0}^{T_{int}}E(i_{n}(t_{1})i_{n}(t_{2}))dt_{1}dt_{2}
 \\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\int_{0}^{T_{int}}\sigma_i^2\delta(t_{2}-t_{1})dt_{1}dt_{2}
 \\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\sigma_i^2dt_{2}
 =\frac{1}{C^{2}}\sigma_i^2\cdot T_{int}\end{align}$$

Here $\delta$ denotes the dirac delta impulse function.

This is a similar result obtained to what was obtained in the post linked above. The problem with the above equation is that it is dimensionally inconsistent (there should be a $T_{int}^2$ factor in the numerator to make it consistent.)

Can someone kindly point out the error in the above argument?


texts: Contradiction between complex baseband and real-valued baseband

texts:
Suppose a direct-conversion transmitter+receiver (ideal transmitter) along with its filters. In complex baseband, the filters in the signal chain can be modeled as a complex-valued FIR filter ($y$, $c$, $x$ are complex-valued!):

$$
y = \sum_{k=0}^K c_k x[n-k]
$$

Rewriting this in cartesian coordinates gives:

$$
y_i + j y_q = \underbrace{\sum_{k=0}^K \left(a_k x_i[n-k] - b_k x_q[n-k]\right)}_{y_i} + j \sum_{k=0}^K \left( b_k x_i[n-k] + a_k x_q[n-k] \right)
$$

Now only considering the real part, it can be seen that it consists of two different filters with which the real and imaginary input signals are filtered.

Now consider the same system from a practical perspective: After the DAC and upconversion, the RF signal looks as follows:

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

The output sequence $y_i$ is obtained by multiplying $x_{rf}$ with $\cos(\omega_c t +\phi)$ and filtering it with a filter called $H$:

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

The filter $H$ (or $G$) is modeled as a FIR filter:

$$
y_i = \sum_{k=0}^K \left( g_k x_i[n-k] + g_k \tan\phi x_q [n-k] \right)
$$

From the equation above it can be seen that the real and imaginary part are filtered through only one filter (they differ just by the constant $\tan\phi$ !).

Where does this contradiction come from?

PS: I know that the first approach is the correct one because it gives the correct results. I do not understand why it is not consistent with my second approach


texts: Kalman Filter: continuous state space, discrete observations

texts: This is just an idea. How can we model the kalman filter to get the state representation in continuous space when the observations to the system are actually from the discrete space. The discrete observation, for example, can be sequence of digits (0-9) where the digits themselves do not have numerical meaning but the sequence of observations are correlated. I know we can try Hidden Markov Model (HMM) for this, but I would like to get the continuous (real-valued vector) space state representation rather than the discrete space as in HMM.


texts:
I'm doing an assignment for the course Signal analysis where I have to analyse a signal. I've tried quite some things now but it's still bothering me that the FFT is looking weird, and not looks like the 'normal look' FFT's we learned in class.

FFT (absolute values of complex values):

FFT (no absolute values):

The FFT seen in the image above is zoomed in on the frequency range 0-30Hz. The rest of the frequency range does not show a lot of (high) peaks, which probably are caused by noise.

The signal is created during a method of welding, using an oscilloscope with a sampling frequency of 1000Hz. I've filtered the signal to remove noise, and after that the signal is converted to the frequency spectrum using the fft function of MATLAB.

Signal before and after filtering:

My general question is, can the shown FFT be valid or did I make a mistake? I estimated the ground frequency to be around 5.5Hz, can I say this when I take one period of the big sinusoidal wave? I also noticed there are about 64 little sinusoidal waves inside one (ground??) period, is this an high harmonic wave form?.

If my theory is right, what causes the fft to be a damped sinusoidal form?

The code I use is basically the following. I leave the part of the noise filtering out because I don't think it's necessary for this question. The dataset is an matrix of 40100 rows.

fs = 1000;
cleanSignaal = data(:,4);
fftSignal = fft(cleanSignaal)/lenght(cleanSignaal);
f = fs/(2*length(fftSignal)):fs/length(fftSignal):fs;
plot(f,abs(fftSignal));
xlim([0 fs(m)/2]);
title('Fast Fourier Transform')
xlabel('Frequentie (Hz)')
ylabel('Magnitude')


texts:
I had a question about normalizing OFDM transmit power after IFFT across each subcarrier. I have simulated OFDM with BPSK and according to different website I have used NFFT/sqrt(Nsubcarrier_used)*ifft(x) where x is the modulated signal and NFFT is the size of ifft and Nsubcarrier_used is total number of subcarrier that is used. This gives me the transmit power of $1$ across each subcarrier.

My question is: How to know the scaling factor, such as NFFT/sqrt(Nsubcarrier_used), of different modulation schemes?

Because I have tried the same scaling factor for 16-QAM here and it did not give me $1$ as a transmit power. Therefore I would really appreciate if anyone can explain about the method to find the scaling factor for FFT and IFFT in OFDM.


texts: Which (Fourier transform?) processing method to use for time variant audio processing? (newbie)

texts:
Edit:
I've since been able to understand the how the relvant parts of Fourier transform work and how they relate to my problem / what I want to do.
And have come to the conclusion that Fourier transform is not the best suited to my problem and that I have an alternative way to do what I want that is much simpler to control / behave in the way I want in practice.

Since I'm new here I don't know what's the right thing to do with my question on this site. Just in case anybody was following this / was preparing an answer, I'm leaving this question up for one more day with this edit and after that I'll delete the question.

I've devloped a basic algorithm for processing audio. I have some programming experience but none with Fourier transforms (but will dedicate all of 2020 to learn this full time if need be).
My question relates to how best to program my algorithm in practice.

What I need to do is the following:

Analyze the frequency content at exact locations of a 44.1kHz PCM
audio stream / file.
Detect and modify the amplitude of specific frequencies at specific
time locations (I do not need to modify phase or have knowledge of
phase to do this).
Then turn this modified version back into audio without audible artifacts from processing errors/limitations.

As I said, I currenty do not know much about Fourier transform. So far my understanding is that an FFT (with windowing) will allow me to analyze the frequency spectrum of a block, detect and modify the amplitude of specific frequencies. But not modify the amplitude of specific frequencies at a specific time in that block and not the same frequency at a different time in that same block? For instance if for 44.1kHz sample rate and block size 8192 as far as I understand it normal FFT will not allow me to detect for instance two seperate 10kHz transients in that block as seperate and modify only one of them?

If my above understanding is correct. Can anybody suggest proper (Fourier transform like?) processing method(s) to use / investigate to achieve my goal?

I've throught of a rough way to do what I want with FFT/IFFT by dividing the frequency spectrum (5Hz to Nyquist for instance) into seperate octave bands all processed seperately, for instance first band 11.025kHz to 22.05kHz under 44.1kHz sample rate and block size 4 for instance not sure of the right block size just an estimate (so I get the transients in one block and still sufficient frequency resolution for that frequency band), then 5.5125kHz to 11.025kHz band under 22.05kHz sample rate (downsampled) also with block size 4 as an estimate, etc with a total of 12 bands down to about 5Hz for instance.
But this will be a lot of work and will take me a lot of time to develop and surely someone must have done something similar and smarter a long time ago to do the type of processing I wish to do :) If anybody can point me in the right direction (and perhaps if I got some assumptions wrong correct me) I will be very thankful!


texts: Software Defined Radio samples analysis using gnuradio and python

texts:
I'm trying to replicate what gnuradio is doing in "pure" python with a hackrf. Having mixed result, I decided to use the same data to compare the outputs I obtain. My gnuradio workflow is the following:
It aims to create two files: one i'll plot directly with matplotlib and python (the "out" file sink), and one that has the format of the hackrf_transfer output (the "out.iq" file sink) that i'll

modify to do the multiplication and complex conjugate thing that was done in gnuradio;
plot using matplotlib and python.

To plot the first file, I use:

#!/bin/env python3
import scipy
import matplotlib as mp
import matplotlib.pyplot as plt
import numpy as np

f = scipy.fromfile(open("out"), dtype=scipy.float32)

window_len = 5
w=np.ones(window_len, 'd')
y=np.convolve(w/w.sum(), f, mode='valid')

plt.plot(y)
plt.show()

To plot the second file, I use:

#!/bin/env python3
import scipy
import numpy as np
import matplotlib as mp
import matplotlib.pyplot as plt

raw_data = scipy.fromfile(open("out.iq"), dtype=scipy.uint8)
raw_data = (raw_data.astype(np.float32)/255)
q = raw_data[::2]
i = raw_data[1::2]
N = len(i)

i = 2*(i-0.5)
q = 2*(q-0.5)

signal = np.zeros(N)
signal = i + 1j*q

energy = signal * np.conj(signal)

window_len = 5
w=np.ones(window_len, 'd')
y=np.convolve(w/w.sum(), energy, mode='valid')

plt.plot(2-y)
plt.show()

It seems to me I'm doing the same operations on the data in the two cases, yet I obtain really different data visualization:
First case:

Second case:

Do you know where I messed up?


texts: Delay of non-polyphase vs polyphase FIR filters

texts:
I implemented a polyphase filter function for interpolating and decimating. It works great (I compared it to what I'd get by interp/decim by inserting L-1 zeros and I get the same results) except for when I compensate delay.

I mean, what I'm doing basically is pulse-shaping with a b=rcosdesign(0.35,4,3,'sqrt'); matlab built-in function (13 taps, so delay is 6). Then I proceed as follows:

Non-polyphase

Upsample by M=3 then x=filter(b.*M,1,upsampled_data). After that, I filter again r=filter(b,1,x) and then I downsample downsample(r(mean(grpdelay(b))*2+1:end),M). As you can see, I'm removing the delay generetad by both filters.

Polypase

I have my custom function so in two lines I run:

x_p = polyfilter( data, b.*M, M, 'interpolation' );
rd_p= polyfilter( x_p(mean(grpdelay(b))*2+1:end), b, M, 'decimation' );

But when I plot the results, only the non-polyphase process 100% matches the original data, while the polyphase one doesn't. You can see such plot here:

When I'm not taking the delay into account, both non-polyphase and polyphase filters output the exact same thing. I kind of understand that I'm not removing the delay at the exact same part of the process in both cases, because I'm removing it right before removing L-1 zeros but after filtering in the non-polyphase case, and in polyphase case I'm removing the delay right before filtering, and I think the correct one is the non-polyphase one but I'm not too sure when I should remove the delay in the polyphase case.

Thanks.


texts: Estimating the deformation of a printed label in an image

texts:
For a university honours project, I am developing an algorithm to read the text on nutrition labels.

One difficulty I have encountered is how to 'flatten' the text in images of labels that are very deformed, such as the one on this ball of mozzarella:

A lot of research has been done on shape estimation of printed documents using generalised cylindrical models, but I don't really think that model is applicable here. It would also not work for images of nutrition labels on chip packets, for example.

However, by virtue of the rectangular nature of the label (when originally printed), I had the idea that it might be possible to approximately estimate the deformation by looking at the borders and lines in the table. I am willing to assume the presence of such lines in any image I need to flatten.

Can anyone point me to some relevant research or give any other advice?


texts: Should zero element be considered to calculate mean when smoothing spectrogram by window?

texts:
My (mass spectrometry) spectrogram does not measure the intensity constantly. I obtain the spectrogram, for example,

time(seconds) = [0,1,2,3,4,5,6,7,8,9]
intensity_profile = [0,10,0,0,50,75,0,60,20]

The 0 in intensity is when the machine does not do measurement (but in fact, it has some signal). If the window is 3 seconds and I choose to use mean, should I average only non-zero element?

smooth_intensity_profile = [10,62.5,40]

Or I should count zero as elements for denominator, and the result would be

smooth_intensity_profile = [3.33,41.67,26.67]


texts:
So, I was trying to plot BER (Bit Error Rate) vs. Eb/No (in dB) curve using Python. Here is what I get,

Eb/N0 is in dB on the x-axis

But I know that this is not the correct graph.

Here is my code.

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erfc

def Q(x):
    return 0.5*erfc(x*(2**(-0.5)))

ebnodb = np.arange(0.1, 12, 0.01)
ebno = [10**(i*0.1) for i in ebnodb]
ebno = np.array(ebno)
y = Q(np.sqrt(2*ebno))

plt.plot(ebnodb, y)
plt.xlabel("Eb/N0")
plt.ylabel("BER")
plt.show()

Can someone figure out, what am I doing wrong?


texts: Can a software defined radio with a frequency range of 300 MHz – 3.8 GHz detect AM (535 to 1705 KHz) and FM signals (88 to 108 MHz)?

texts: I am interested in knowing whether BladeRF SDR with a frequency range of (300 MHz - 3.8 GHz) detect AM and FM frequencies? Also, how important is the sampling rate of the SDR. For example, BladeRF has a sampling rate of 80 KSPS - 40 MSPS. What does this number signify? will every signal going out/in the SDR be sampled with this sample rate? Similarly, what does the maximum bandwidth signify?


texts: Avoiding edge repeats in a FFT convolution

texts:
I am trying to use FFT as a fast way to convolve a image, however one of the downsides of this is that the "kernel" image will repeat when reaching the edges of the sample.

Is there a way to avoid this to make it more like a classic convolution, that cuts off at the edges?


texts: What does this notation $g(x) = S\{\delta(t) \}(x)$ mean?

texts:
What does this notation mean $g(x) = S\{\delta(t) \}(x)$ mean, where $S$ is a "system", $\delta(t)$ is the impulse function and $g(x)$ the output function of the system? I am really not familiar with signals and signal processing. I know the definition of the impulse function

$$
\delta(t) =
\begin{cases}
1, \; t = 0,\\
0, \; t \neq 0
\end{cases}
$$

I understood that a system is a function that gets as input a function and returns another function (so it should be a functional). However, I do not get this notation $g(x) = S\{\delta(t) \}(x)$? For example, why the parentheses $\{ \}$? Why do we apply the $\delta(t)$ inside $S\{ \cdot \}$? What is the precedence of those function applications? How should I read that formula? Why do we even use an impulse here, what is the intuition behind this?

In my case, also have that

$$g(x) = S\{\delta(t) \}(x) = (A\delta(t))(x) = A_{x, 0}$$

Why are the notations $S\{\delta(t) \}(x)$ and $(A\delta(t))(x)$ equivalent? What is this $A_{x, 0}$? $A$ should be a Toeplitz matrix, but I don't understand why. How should this matrix $A$ look like?

$A_{x, 0}$ should represent coefficient associated to instant $0$ for computing the output at instant $x$. I also don't get this, even though I think we are trying to convert a continuous signal to a discrete one.

You can assume that $S$ is a linear shift-invariant system.


texts: When integral control is added to proportional controller (i.e. PI control), why steady state error becomes 0 and  why oscillations appear , how actually the integrator works ?


texts: Controlling a signals amplitude in the frequency domain

texts:
I want to control the amplitude of a signal I'm creating from a user drawn spectrum by scaling the magnitude values in the frequency domain.  Here is my scenario.

Sample rate $F_s= 44100\textrm{ Hz}$
FFT size $NFFT = 512$
Desired waveform frequency: $86.1328125\textrm{ Hz}$ ($44100/512$ so a single cycle)

My user input screen has the magnitudes presented like a bar chart and they are all stored as values with a range of $0.0$ to $1.0$.

How do I scale these in the frequency domain so that the output signal is a $0\textrm{ dB}$ signal?


texts: Why use both a high pass and a low pass filter in a Butterworth implementation for noise reduction?

texts:
I am a new to audio processing .  I want to develop an application that reduces audio noise using a Butterworth filter. I found some existing code doing this, but I still do not understand the use of 2 filters (low and high pass filter) in a reversed order.  What I know from my experience, the high pass filter will eliminates the low frequency samples, so there is no need to apply the low pass filter.
The code begins with :
b,a = scipy.signal.butter(5, 1100/(Frequency/2), btype='highpass') # ButterWorth filter highpass
filteredSignal = scipy.signal.lfilter(b,a,NewSound) # Applying the filter to the signal

and ends with:
c,d = scipy.signal.butter(5, 100/(Frequency/2), btype='lowpass') # ButterWorth low-filter
newFilteredSignal = scipy.signal.lfilter(c,d,filteredSignal) # Applying the filter to the signal


texts: Reducing number of taps in FIR for echo cancellation

texts: I'm trying to implement the NLMS algorithm on an FPGA. The problem is that I'm only allowed to use a maximum number of 2048 coefficients for my FIR filter, my sampling rate is 16KHz, however, the room impulse response sampled at 16KHz can be 10 times longer than 2048, is there any way to decrease the number of taps to remove long echos with only 2048 taps? I'm new to FIR subject, I know that my FIR filter IP allows decimation, interpolation, and can have different sampling rate than the input, can I exploit this to decrease the number of taps? I can also compress room impulse response and round many samples that are close to zero and make them zero but can this help to reduce the number of coefficients?


texts: Random (Over) Sampling signal and perfect Reconstruction in Nyquist form?

texts:
Imagine we have band limited signal with bandwidth of $B$, so the required Nyquist rate would be $f_{nyq}>2B$ that is oversampled with rates $f_s$ where $f_s = M*f_{nyq}$ and $M$ is random and $M>=1$. By that, I mean the sampling rate is at least at Nyquist rate but most of the time its way more than Nyquist and it is random. I have heard and read about  iterative  My questions are:

Is there a method to construct the signal as if it was taken in Nyquist more, without degradation of signal quality? After all, we have more samples than needed, so I do not expect any degradation, right?
Is there are some loss, how much is it, what factors it depends on?
Is there a toolbox or something that I can feed my data in and have some experiments of my own? ,

4.Is the technique used in MATLAB lossy?

I appreciate any kind of contribution here :)


texts: Series connection of digital PT1 filters

texts:
I was trying design $N^{th}$ order filter, I was suggested a series connection of low-pass PT1 filter. Basically, the filter coefficients are like this:

$$a_i = \dbinom{n}{i} T^{i}$$

This works for $N^{th}$ order, while $a_1, a_2, a_3....a_n$ being the filter coefficient.

I don't know how this works, in theory. I tried searching, but could not come up with anything. if anyone has any idea about this. I want to know how this works.


texts: How to Use the Walsh Hadamard Transform to Filter 24 bit RGB Color Values in a Localized Spatial Region?

texts:
After reading this Signal processing Stack Exchange article tonight, What is the Walsh-Hadamard Transform and what is it good for?, I became interested in applying the Walsh-Hadamard transform to filter one of three classes of 24 bit RGB color values, True Blue or True Green or Tree Red, in a localized spatial region of the color image seen from a commercial airplane cockpit.

In the book Color Image Processing and Applications, the authors Konstantinos Plataniotis, Anastasios N. Venetsanopoulos state that the RGB space is not sn efficient representation for compression because there is a significant correlation between the three color components since the image energy is distributed equally among them both spatially and spectrally.

A solution is to apply an orthogonal decomposition of the RGB color signals in order to compact the color image data into fewer channels. The commonly used YIQ, YUV and YCBCR color spaces are examples of color spaces based on these principles. The resulting luminance chrominance decomposition exhibits unequal  energy distribution from the luminance component in which the vast majority of fine detail high frequencies can be found

I would like to apply color image YUV or YUV (or YCbCrYCbCr which is like HSV) compression to theoretically attain almost real time filtering of RGB color coordinates in a localized spatial region of the color image seen from a commercial airplane cockpit. The reason for the spatial localization is that we do not want to subtract green from everwhere uniformly in the color image projected to the commercial airliner flight team.

Please consider the paper, NTSC component separation via Hadamard transform,
Published in: Vision, Image and Signal Processing, IEE Proceedings -  (Volume:141 ,  Issue: 1 ) Feb. 1994 page 27-32 where the authors, M. Coley and J.Barba of the New York Institute of Technology present a technique for separating the NTSC composite signal into its components (Y, I, Q) by applying an N×N Hadamard transform directly to the composite signal sampled at four times the colour subcarrier.  They show that the Y, I, Q component signals are mapped into specific areas in the Hadamard domain. Component separation is achieved by assigning particular transform coefficients to each component signal. Simulation results of component separation applied to five colour images and one black and white test pattern
are presented. Both signal to noise ratios and subjective results are presented. The results indicate that no noticeable degradation is expected in typical colour images.

  I am interested in restricting RGB color filtering to small localized
  spatial regions of the color image acquired by per pixel processing of
  an array of sensors developed by either by M.I.T Lincoln Labs [27]
  Kelly, M. W. and Blackwell, M. H., “Digital-pixel FPAs enhance
  infrared imaging capabilities,” Laser Focus World 49(1), 90 (2013).
      or Stanford University's DEpartment of Electrical Engineering and Computer Science, [28] Wan, G., Li, X., Agranov, G., Levoy, M., and
  Horowitz, M., “CMOS Image Sensors With Multi-Bucket Pixels for
  Computational Photography,” IEEE J. Solid State Circuits 47(4),
  1031–1042 (2012).

Thank you.


texts: Relation between probability of a signal and its magnitude spectrum

texts:
I am a newbie. In these days I am searching if there are reference about properties that can relate the probability of the signal and the magnitude spectrum. I know what are both. But I was curious wheter there are connection between them. Can you help me
Given a signal y(t), the magnitude spectrum  is
$$Y(f)=\left\lvert\int y(t)e^{-2\pi ft}\,\mathrm dt\right\rvert,$$
in other words the absolute value of the FT of the signal.


texts: How to Apply Window Functions in Frequency Domain?

texts:
I am trying to apply a window function to smooth my signal and finally obtain the peaks of the signal. Yet, for spectral leakage I apply window function in time domain. After taking the Fourier amplitude spectrum ı want to further smooth the signal. Yet I do not know how to apply smoothing windows for frequency domain. Should I use convolution for the Fourier amplitude spectrum of the window and Fourier amplitude spectrum of the data?

What is the meaning of smoothing window length with 0.1 Hz or 0.2 Hz?


texts: Why is image blurring called convolution and not cross correlation?

texts:
I've always thought that image convolution is a weighted sum of the nearby pixels using the convolution kernel as weights. But from the description on Wikipedia, this is what cross correlation does. Convolution seems to be equivalent to 1) reversing the convolution kernel 2) doing cross correlation with this reversed kernel.

If I understand it correctly – why is image blurring called convolution? What's the point of reversing the kernel?


texts:
I got this image:

I need to keep only the part in the center (I mean, the circle with texture inside). How would you do it in Python? I'm trying to use Hough Transform, but some times it detect wrong circles when using other similar images, and I need the right circle all the time.
Also, the solucion would have to work in different images, each one with different size and location of the circle.
(For those who are interested, is a blastocyst image. I have the original too, but I guess this processed one is easier to handle).


texts:
I have the following code, which I am told is a low-pass Butterworth filter. How do I derive the frequency response and phase for this? I only have experience in analog electronics and so cannot relate this to a real filter response. I'm not sure why everything is scaled up by 2^ 16 - it must be to do with maintaining adequate resolution during the calculations?

#define FILTER_TC(FC, FS) ((411775 * FC) / FS)'

where TC = (2^16) * 2pi * Fc * Fs, with Fc the cut-off frequency in Hz, and Fs the sample frequency in Hz.

#define FILTER_LP(FILTER_OUT, FILTER_IN, FC, FS) \
    (FILTER_OUT += ((signed long)FILTER_TC(FC, FS) * \
    ((signed long)FILTER_IN - (signed long)(FILTER_OUT)) >> 16)'

Thank you for any help with this.

jmv


texts: How to sketch/plot this discrete time signal?

texts:
$(n+5)u[n+5]+2nu[n]+(n-10)u[n-10]$ is the question.

The $u[n]$ is messing me up. If I substitute $n$ as $0$, what would be $u[5]$ or $u[-10]$?


texts: How to approximate gaussian kernel for image blur

texts:
From wiki, a $3 \times3$ gaussian kernel is approximated as:
$$\frac{1}{16}\begin{bmatrix}1&2&1\\2&4&2\\1&2&1 \end{bmatrix}.$$
Applying this kernel in an image equals to applying an one-dimensional kernel in x-direction then again in y-direction, so one-dimensional kernel is
$$\frac{1}{4}\begin{bmatrix}1&2&1\end{bmatrix}.\qquad \text{for 3*3 kernel}$$
$$\frac{1}{16}\begin{bmatrix}1&4&6&4&1\end{bmatrix}.\qquad \text{for 5*5 kernel}$$
My question is how to derive this approximation?


texts: What are the mathematical consequences of a discrete signal living in Hilbert space?

texts:
I am doing the Coursera course on DSP by Prandoni and Vetterli, and was excited to learn that we can use the tools of vector spaces to analyze discrete signals. At this early point in the course the professors did mention that since the signals live in Hilbert space, we can naturally talk of orthogonality, for one.

I am curious: what other linear algebra operations/concepts do we use with discrete signals, and what is the rationale behind them?


texts: How do decide which order my butterworth filter should have?

texts:
I have a signal that was recorded with around 300 samples in a second. Now I want to run a low pass and a high pass butterworth filter over it. The frequency for the high pass filter is around 0.05 Hz and the one for the low pass around 5 Hz.

But how can I decide which order my butterworth filters should have? Is there any formula?


texts:
I am trying to design a FIR bandpass filter for a STM32F407 microcontroller with a passband of 1kHz between 59500 Hz and 60500 Hz or something similar (since I don't really know how tight I can squeeze it). I am trying to design the filter in MATLAB so I can get the co-efficents for later implementing them in C.

While simulating it in MATLAB the problem seems to be with the high sampling frequency. The microcontroller that I'm using is sampling at 5.25 MSPS so the FIR filter will have a pretty high filter order and I am wondering if there is something that I can change or can a filter really have this high order? Because it kind of looks like it would be working if the order is around 3000-4000.

Am I doing this correct and are these parameters even realistic?

% MATLAB Code
% Generated by MATLAB(R) 8.3 and the DSP System Toolbox 8.6.
% Generated on: 28-Jun-2016 09:43:38

% FIR Window Bandpass filter designed using the FIR1 function.

% All frequency values are in Hz.
Fs = 5250000;  % Sampling Frequency

N    = 3500;     % Order
Fc1  = 59500;    % First Cutoff Frequency
Fc2  = 60500;    % Second Cutoff Frequency
flag = 'scale';  % Sampling Flag
% Create the window vector for the design algorithm.
win = blackman(N+1);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);
Hd = dfilt.dffir(b);
%----------------------------------------------------------

t = 0:1/Fs:0.01;    % time vector

signal1  = sin(2*pi*58000*t);
signal2  = sin(2*pi*60000*t);
signal3  = sin(2*pi*62000*t);
signal4  = sin(2*pi*80000*t);

hold on;

length = size(t);

subplot(4,1,1);
plot(filter(Hd,signal1));
axis([0 length(2) -1 1]);
title('58kHz');

subplot(4,1,2);
plot(filter(Hd, signal2));
axis([0 length(2) -1 1]);
title('60kHz');

subplot(4,1,3);
plot(filter(Hd, signal3));
axis([0 length(2) -1 1]);
title('62kHz');

subplot(4,1,4);
plot(filter(Hd, signal4));
axis([0 length(2) -1 1]);
title('80kHz');

hold off;


texts: Discrete wavelet transform; how to interpret approximation and detail coefficients?

texts:
I am trying to understand Wavelet transform. So far I have understood the basic theory of it. But I am not able to get my head around how to interpret both coefficients.

I am using PyWavelets package of python, I have a time-series data for 1 year equally divided into 15-min time interval. I am trying to get high-frequency components for each day. Fourier transform gives the information about high-frequency. But I would like to get time-information also i.e. at what part during the day there occurs a spike(high-frequency). If someone could please help me figure out this and how to plot time-frequency. I have searched a lot but couldn't find a perfect answer. I am not looking for a free solution without trying; Here is everything I have tried

tdayde=d2012_2['Date'].iloc[0]
endde=d2012_2['Date'].iloc[-1]
skiphr=datetime.timedelta(hours=24)
blockde=tdayde+skiphr

I have timestamps attached to data, so to plot for each day I am using

healtharrayde=d2012_2
healtharrayde=healtharrayde[(healtharrayde['Date'] >= tdayde) & (healtharrayde['Date'] < blockde)]
sig=list(healtharrayde['MCP'])
  cA, cD = pywt.dwt(sig, 'db2')
    plt.plot(cD)
    plt.plot(cA)

So if there are 96 values in array, both cD and cA are 48 values each. I would like to understand and plot cD(detail components) with respect to time, So how should I do? because there are 96 timestamps(MCP is column name for values)

Sorry if it's a silly question but I really need help. Thank you