texts: Frequency analysis to determine low-pass cut off frequency

texts:
I collected some data for a practical application, where the signal represents force data obtained from an impact of a punch against a force plate attached to a quasi-rigid rig (it moves once the impact occurs). I have a few questions to understand how to deal with my data.

Features of the signal:

the signal (sampling 1000Hz) has a length of 2 seconds, but the impact last about 60ms. The recording is 2s for practical reasons though. Therefore, before the impact a baseline is recorded (values close to zero), and after the impact vibrations occur until data set back to baseline levels. I need to figure out a reasonable cut-off frequency to use combined with a 4th order butterworth filter to properly filter the signal.

When exploring (FFT) the frequency of that signal should I use the whole 5s long signal or only the 30ms of interest? Due to sampling frequency, data points over 30ms maybe not enough?
Vibrations occurring due to the impact are not of interest, but the frequency is likely to be in part similar to the frequency of the portion of the signal of interest. How to deal with this? I could select from before the impact (arbitrary) to when the signal becomes zero after the impact. Then use from that point onwards to explore frequency of the unwanted signal only?
Impacts are at least 10 per subject, when exploring the frequency of that signal, should I average all trials and all subjects out to obtain an average frequency content or should I perform the analysis for each trial and subject indepedently?

Attached a figure of a typical signal (x = time; y = Newtons)


texts: what is the difference between QPSK and 4QAM?

texts: i need to know the main difference between QPSK and 4QAM, i have searched over the internet and still i have not found a good explenation.


texts: Twisted convolution - simplifying multiple integral via substitution

texts:
I need help with understanding a proof from a paper. Following is the description:

Let
$$r_{1}(t) = \int \int h_{1}\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} s\left(t - \tau\right) d\tau d\nu$$
$$r(t) = \int \int h_{2}\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} r_{1}\left(t - \tau\right) d\tau d\nu$$
Substituting $r_{1}(t)$ in $r_{2}(t)$  we obtain after some algebraic
manipulations,
$$r(t) = \int \int h\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} s\left(t - \tau\right) d\tau d\nu$$
where, $f(\tau, \nu)$ is given by
$$h(\tau, \nu) = \int \int h_{2}\left(\tau', \nu' \right) h_{1}\left(\tau - \tau', \nu-\nu' \right) e^{j 2\pi \nu' \left(\tau - \tau' \right)} d\tau' d\nu'$$

I am unable to use appropriate algebraic manipulations to arrive at the final expression for $r(t)$. Can anyone help me? The first step is arrived at by substituting $r_{1}(t-\tau)$ in $r_{2}(t)$. But I am unable to proceed further. (Incorporated corrections from Tanya Choudhary's answer)
$$r(t) = \int \int h_2(\tau,\nu)e^{j 2\pi\nu(t-\tau)} \int \int h_1(\tau_1,\nu_1)e^{j 2\pi\nu(t-\tau-\tau_1)}s(t-\tau-\tau_1)d\tau_1 d\nu_1 d\tau d\nu$$
Substituting $\tau = \tau', \tau_1 = \tau - \tau', \nu = \nu', \nu_1 = \nu - \nu'$,
$$r(t) = \int \int h_2(\tau',\nu')e^{j 2\pi\nu'(t-\tau')} \int \int h_1(\tau - \tau',\nu - \nu')e^{j 2\pi\nu(t-\tau)}s(t-\tau)d\tau' d\nu' d\tau d\nu$$


texts: Do Causal Discrete-time systems have proper transfer functions?

texts:
In the case of continuous-time systems, if the system is causal, its Laplace transfer function is strictly proper (the degree of the numerator is less than the degree of the denominator).

Is this true also in case of discrete-time systems? That is, if the system is causal, is its Z transfer function strictly proper?


texts:
I was reading a paper which uses a filter bank of gabor filters with different frequencies and orientations. However, I often see the construction of a filter bank of gabor filters with different frequencies and scales.

I am not sure on the construction of it. How do I construct it, any generic example will be ok. Thanks!


texts: Sensor fusion under unknown correlations: can covariance intersection account for delays?

texts:
Of late, there has been some interest in cooperative estimation algorithms in robotics, where the information sources are usually sensors such as cameras. When multiple robots observe surrounding environments through cameras, they can share information about salient features and correct the poses of each other. At the same time, these measurements are not independent because the source of information may have been the same, or the robots may have exchanged some data in the past and have not kept track of cross correlation terms.

In this context, there have been many papers discussing covariance intersection as a way to fuse relative measurements with individual measurements. In my example for this question, there are two robots $i$ and $j$ which are capable of performing their own position estimation at every time step. Occasionally, robot $i$ can compute the relative position to robot $j$. Say at time step $k$, $i$ generates a relative measurement to $j$ (say, $z_k^{ij}$) and using this, generates a propagated state-covariance pair for the robot $j$ as follows:

\begin{gather}
    \mathbf{x}_k^{j'} = \mathbf{x}_{k|k-1}^{i} + \mathbf{M}_{k}^{i,j} \mathbf{z}_{k}^{i,j} \\
    \mathbf{P}_{k}^{j'} = \mathbf{H}_{k}^{i,j} \mathbf{P}_{k|k-1}^i \mathbf{H}_{k}^{{i,j}^\top} + \mathbf{M}_{k}^{i,j} \mathbf{R}_{k}^{i,j} \mathbf{M}_k^{{i,j}^\top}
\end{gather}

On the other side, the robot $j$ already has a locally propagated state and covariance estimate, and once the above pair is received, it fuses this relative data with its own estimate using covariance intersection.

\begin{gather}
    \mathbf{P}_{k|k}^j = {\Big[\pmb{\omega}(\mathbf{P}_{k|k-1}^j)^{-1} + (1-\pmb{\omega})(\mathbf{P}_{k}^{j'})^{-1}\Big]}^{-1} \\
    \mathbf{x}_{k|k}^j = \mathbf{P}_{k|k}^j \Big[\pmb{\omega}(\mathbf{P}_{k|k-1}^j)^{-1} \mathbf{x}_{k|k-1}^j + (1-\pmb{\omega})(\mathbf{P}_{k}^{j'})^{-1} \mathbf{x}_k^{j'} \Big]
\end{gather}

My questions:

In real systems, measurements and corrections, especially when involving communication, has a significant delay. While there has been some discussion on how to incorporate delayed measurements in the conventional KF update steps, how would a CI based Kalman filter deal with delays? Say robot $i$s state covariance pair for robot $j$ reaches robot $j$ not at time step $k$, as would have been ideal, but at time step $k+\delta$. Can a modification to the CI algorithm take that into account? Assuming that robot $j$ has access to its own measurements meanwhile.
What happens to the usual update step of robot $j$ when it receives a relative measurement? Does it happen after the CI fusion step? In that case, should this also be treated as a 'delayed' measurement?


texts: Does the determination of the ergodicity of a signal force any changes in methodology?

texts:
In mathematics, the term ergodic is used to describe a dynamical system which, broadly speaking, has the same behavior averaged over time as averaged over space.  -from wikipedia

From the perspective of engineering a signal processing system, does knowing whether a signal is ergodic change the plan of attack for analyzing the signals?  I've always found this a fascinating concept, but I don't really know what do do with it once the determination is made.  What analysis methods would be more pertinent knowing this information?


texts: Are discrete time functions expected to be this inaccurate?

texts:
I am attempting to test some discrete time functions against their continuous counterparts and I am surprised by how inaccurate the discrete versions are coming out.

I am modeling the displacement, velocity, and acceleration of an ideal string struck by a mass as per: https://www.dsprelated.com/freebooks/pasp/Ideal_String_Struck_Mass.html

This gives displacement as:

$y(t) = \frac{v*m}{2R}*(1-e^{\frac{-2Rt}{m}})$

Velocity as a function of time is:

$v(t) = v*e^{\frac{-2Rt}{m}}$

Acceleration can be derived from velocity:

$a(t) = \frac{-2Rve^{\frac{-2Rt}{m}}}{m}$

A discrete time formula for velocity using the displacement output from the first continuous function is:

$v(s) = sy(s)$

$v[n] = \frac{y[n] - y[n-1]}{T}$

A discrete formula for acceleration from displacement is:

$a(s) = s^{2}y(s)$

$a[n] = \frac{y[n-2] - 2y[n-1] + y[n]}{T^{2}}$

Mostly the discrete formulas seem to give normal output that matches the continuous equations. It seems close at the start and finish. However, at some points it's shocking how far off it is.

Here's some basic code to test this:

displacement_2 = displacement_1;
displacement_1 = displacement;
displacement = ((velocity * mass)/ (2 * impedance)) * (1 - exp(-2 * impedance * timer / mass));

velocityContinuous = velocity * exp((-2 * impedance * timer) / mass);
velocityDiscrete = (displacement - displacement_1) / deltaTime;

accelContinuous =  (-2 * impedance * velocity * exp((-2 * impedance * timer) / mass)) / mass;
accelDiscrete = (displacement_2 - (2 * displacement_1) + displacement)/(deltaTime*deltaTime);

timer = timer + deltaTime;

...
deltaTime = 1/sampleRate;
velocity = 8.f;
impedance = 0.5f;
mass = 0.004f;

At some points I'll get output like:

 velocityContinuous: 2.81216 velocityDiscrete: 3.47256
 velocityContinuous: 2.78302 velocityDiscrete: 3.45452
...

 accelContinuous: -840.176 accelDiscrete: -1089.63
 accelContinuous: -831.469 accelDiscrete: -1083.97
...

That's a pretty massive error. Is this normal? Or have I done something wrong? Sample rate is 96 kHz.

On another note, accelDiscrete puts out 0 first sample, then an insanely large number of the wrong sign second sample, then "correct" output after that. velocityDiscrete puts out 0 first sample then normal output. I assume this is because displacement_2 and displacement_1 are 0 at initialization and it can't calculate properly from them. Is there any typical fix for this type of behavior, particularly accelDiscrete which is putting out a crazy number on the second sample?


texts:
I think this question is so trivial that it's nowhere explained (I would not find it) and the answer is probably trivial too. But I still wonder how codes in a typical CDMA system are distributed? They act like symmetric keys in a crypto system: both parties need to know them and there is a complexity of O(N) to share them (N being all possible users).

For example, in a WCDMA cellphone system, there is an "infinite" number N of possible nodes (e.g. a roaming user coming to the US). The codes can't be statically assigned. And the wireless base station can't just try all possible codes.

So I assume that the first communication of a handset and base station does not go over CDMA but the base station somehow detects that there is a new handset, creates a new code and sends this code to the handset. Afterwards the communication goes over CDMA.
For handover between different stations, the base station could distribute the code among its neighbor stations.

How does it work in practice?
And am I correct that the seed to an LSFR generating the chipping sequence is distributed as opposed to the chipping sequence itself?


texts: audio analog input + usb data cable = sound noises

texts:
I connect my android smartphone from audio jack output to my pc audio input and get a good analog sound signal. But when I connect this smartphone to  the pc with usb data cable the analog sound becomes hissing e.g. music from smartphone overlaps with noises on pc.
Why does this happen?
How can I eliminate these noises? maybe I have to use some shielded cables? If yes - then audio or usb? or I need some sound direct box or something like that?


texts: How to find the group delay of IIR filters in Python/SciPy

texts: I sent a square wave signal (red) through a IIR filter (butterworth in this case), to make it more "realistic". Then there's a slight group delay in my output signal (blue). Is there a way to find out what the group delay is? I only need the total group delay, not spectrum of group delay. I need the number so I can shift my decision point accordingly, so I am still sampling at the middle of the eye diagram to achieve optimum result.


texts: Feature extraction for sound classification

texts:
I'm trying to extract features from a sound file and classify the sound as belonging to a particular category (eg : dog bark, vehicle engine e.t.c). I'd like some clarity on the following things :

1) Is this doable at all? There are programs that can recognize speech, and differentiate between different types of dog bark. But is it possible to have a program that can receive a sound sample and just say what kind of a sound it is? (Assume there's a database containing a lot of sound samples to refer to). The input sound samples can be a bit noisy (microphone input).

2) I assume that the first step is audio feature extraction. This article suggests extracting MFCCs and feeding them to a machine learning algorithm. Is MFCC enough? Are there any other features that are generally used for sound classification?

Thank you for your time.


texts: How to find correlation between 3 signals

texts: I have 3 signals, S1 S2 S3, with a same time domain and same time steps. I'm wondering if there is any correlation between these 3 signals. more specifically I want to know if it is possible to make Signal S3 by any combination of signals S1 and S2.


texts: How can I obtain the transfer function from these two transfer functions in this case?

texts:
Using Scipy signal in Python I want to obtain a transfer function for an input which passes first from sys1 and then sys2. So explicitly
$$\text{input}{\longrightarrow}\boxed{\tt sys1}{\longrightarrow}\boxed{\tt sys2}{\longrightarrow}\text{output}$$
sys1 = signal.lti([a], [b, c])

sys1 = signal.lti([x], [y, z])

How can I cascade sys1 and sys2 programatically to obtain the overall transfer function (call it sys_total) between the input and output?


texts: Reconstructing a partially deleted image through wavelets

texts:
I am trying to form an approximation of the wavelet transform from a partially sampled image.

Reconstruction in the 1D case is easy. We have $w = h x$, with $w$ as the wavelet coefficients, $h$ as the inverse of the Haar matrix and $x$ our signal. $w$ and $x$ are $n \times 1$ and $h$ is $n \times n$. Where we don't sample the signal, we delete those columns from $h$ and those rows/elements from the full signal, $x$. Since we have to have as many equations as unknowns, we sample at $m$ locations and approximate the first $m$ wavelet coefficients. All we do is delete the corresponding columns of where we do not sample.

I have $w = c x r$, where $c$ is the matrix that performs the wavelet transform on the columns, and $r$ the rows. I am having a hard time forming in approximation in this 2D case. I know this is a simple solution and am having a hard time getting the math behind it. I know it's just solving a linear system of equations, but am having a hard time coming up with the matrix to represent that linear system.


texts: How to match a thumbnail to an area in the original image?

texts:
Suppose that you have a thumbnail and the original image from which it was generated. You need to match the thumbnail to an area in the original image that represents the original selection from which that thumbnail was generated. Humans are really good at this sort of matching, but I have no idea how to go about writing an algorithm to do this automatically.

Here are my thoughts on the matter so far:

A thumbnail is created by first cropping then scaling the original selection. Logically, some scaling is required to reverse the process.
A thumbnail's dimensions could be greater than that of the original selection. Therefore, this algorithm should account not only for the common downscaled thumbnails, but also for cases where the cropped original selection has been upscaled to the thumbnail's intended dimensions.
Theoretically, it is possible that a thumbnail was generated such that the entire original image was encompassed inside the original selection, akin to background-size: contain in CSS. All four corners of the original selection thus fall outside of the original image. The extra space could have been filled either with a solid color (possibly transparent) or with Photoshop's Content Aware Fill or a similar algorithm (inpainting). This is an edge case; however, an ideal reverse-thumbnail algorithm should be able to handle such cases.
This algorithm should return the x,y pixel coordinates of the top-left corner of the original selection as it falls on the original image, as well as the width and length of the selection.

I suspect that feature detection is required to make this work, but I don't know where to start. Partial answers are acceptable, but if this is an easier problem than I think it is, please do give it a shot. Any language is fine, but the implementation should not rely on proprietary software, and Python or JavaScript are tentatively preferred.

I've put together a sample set of 30 images from my current project for testing  potential algorithms. I tried to include a variety of images to challenge any feature-recognition approach. All selections in this collection have been downscaled to the intended thumbnail size of 300x300 after cropping. I included an info.csv, which contains id,x,y,w,h representing the coordinates and dimensions of the original selection relative to the full-sized image. These values were obtained by applying Deep Zoom coordinates saved from the source tiff files to the full images. Due to floating point inaccuracy, w and h differ by one pixel in three cases (13, 17, 21). These images were included intentionally to test how your algorithm might handle accuracy concerns. A second info_ur.csv includes unrounded values for reference.

Download test images: https://app.box.com/s/65240dlgvfvevw2grhyaqqjjc29jzv65

Finally, here's an illustration of the four cases such an algorithm might encounter:


texts:
I build a flowgraph in gnuradio, generated python code, change a little and run it using python.
My question is, does building the flowgraph in cpp will provide better results (time)? (in other words does GIL and other stuff that make python slower also infect here?)(and if not i would like to read why)


texts: I'm sure there must be an easy way to do this, but given the Fourier transform of an isotropic filter kernel, $\hat{f}(\mathbf{u}) = \mathcal{F}f(\mathbf{z})$, can one calculate the value of the kernel at $\mathbf{z} = 0$?


texts:
So I have an incoming signal of 137.9 MHz and a sampling frequency of 1.4 MHz. Once sampled how do I calculate the aliased signal's new frequency? Bandwidth is 38 kHz with a subcarrier of 2400 Hz, if that information is relevant. Just looking for an equation to calculate the aliased frequency.

The hardware used is a RTL-SDR into a computer sampling at 1.4 MSPS.

This is the signal's information page:

https://www.sigidwiki.com/wiki/Automatic_Picture_Transmission_(APT)


texts:
I have multiple time series coming from pings of phone applications in a specific area.
I am suspecting that the signal is sometimes bigger and sometimes lower not because the actual signal is bigger but due to technical things (like more people downloaded the application, or the application optimized something in their process so I get more data).
I want to identify these trends in the signal and remove it in order to clean the signal from these inconsistencies.
On one side, I want to get rid of the trend created by the technical issues, on the other side I want keep the signal due to the actual number of phones in the area.
I guess I could use some kind of comparison between the signals and understand which trend is due to the real number and each is due to the provider. But I am not sure how to proceed ?
How would you do it ?
Thanks !


texts: What Is the Purpose of Bundle Adjustment?

texts:
So, my understanding of stereovision is that you essentially have two cameras. Then, offline, you perform camera calibration with checkerboards to get intrinsic camera properties as well as the geometric transformation between cameras. Then, you can take a stereo pair of images, find the disparity map, and using some geometric constraints, can find the 3D positions of points in the scene.

I've read recently about "bundle adjustment" and that it should be performed at the very end? What's the purpose?


texts: Difference of Gaussian filter but using Gaussian CDFs not PDFs

texts:
A very common filter in signal / image processing is the Laplacian of Gaussian filter, and its approximation, the difference of gaussians filter.

The DoG filter is the difference of 2 Gaussian PDFs that have slightly different standard deviations.

My question is: what would a filter that uses the difference of 2 Gaussian CDFs (instead of PDFs) look like and what properties would it have?


texts: RRC filter or rectangular pulse shaping filter in LFM chirp signal

texts: I am working with a LFM up chirp signal, and I wanted to embed some data on it using a BPSK constellation. My doubt was, ¿should I use a RRC filter, like in usual communications schemes, or should I use a rectangular pulse shaping filter? ¿Is there any advantage to use a rectangular pulse shaping filter in this case?


texts: Compute dynamic range in a linear quantization system

texts:
I know the dynamic range is the ratio between the maximum amplitude and the minimum amplitude of a signal, calculating the dynamic range in dB I get a difference of logarithms:

But why the dynamic range of a N-bit digital system with linear quantization is 6N dB?


texts: Why Does the Median Filter Minimize the Absolute Value Error $L_1$ Cost Function?

texts:
I can easily prove that the mean filter minimizes the square error $L_2$ cost function using simple calculus.
However, how do you prove that the median filter is optimal with respect the absolute error $L_1$ norm?


texts: Validating cross corrolation between sine and cos, shouldn't pahse lag be pi/2?

texts:
I am trying to use cross correlation instead of FFT to find lags between 2 signals. I am trying to understand cross correlation by using it on sin and cos and expecting phase lag to be pi/2 but it is not:

clc
clear

frequency =1000;
t= linspace(-1000,0.1,1000);

sine_wave = sin(2*pi*frequency*t);
cos_wave  = cos(2*pi*frequency*t);

[xc,lags] = xcorr(cos_wave,sine_wave);
[~,I] = max(abs(xc));

figure
stem(lags,xc,'filled')
legend(sprintf('Maximum at lag %d',lags(I)))
title('Sample Cross-Correlation Sequence')

I am getting a maximum at lag 22 (if you copy and paste the code provided you ll see it in the plot) and I understand that phase lag = 2*pif (time delay) so in my case:

phase lag = 2*pi*1000*22 which is by no means close to pi/2...
What am I missing here?
Thank you


texts:
I am currently working with a set of wavelet transforms that I came across in a paper. The paper suggests that I first apply the FFT on the source image and then perform the downsampling. So it says something like:

Consider a source image of size 128*128.

P = fft(source);
P(1:64,1:64) = P(1:64,1:64) + P(65:128,65:128);
P(65:128,1:64) = P(65:128,1:64) + P(1:64,65:128);
P = P(1:128,1:64);

Well this solution definitely works and I get the results, but could anyone give a possible explanation as to why I am supposed to add in this manner?


texts:
I have 2D pictures of an airplane of a specific model, taken with a simple camera, and a 3D point cloud model of this airplane.
I would like to project the pictures onto the 3D model (cloud of points) based on control points that I manually locate on both ends. This should allow me to calculate distances on the 2D pictures.
Is there a programming library that does this? I program in Python mostly.


texts: Echo cancelling using autocorrelation function

texts:
I was given a problem, but I couldn't solve. I did some researches but I still didn't figure it out. Here is the problem:

An audio signal $ s(t) $ is generated by a speaker reflects in a wall with reflection coefficient $ \mu (\mu < 1) $. The signal $ x(t) $ is recorded by a microphone close to the speaker and far from the wall, after sampled, is given by:

$ x(t) = s(t) + \mu s(t - k) $

Where $ k $ is the delay given in samples due to the echo.

So, my goal here is to estimate $ \mu $ and $ k $, observing the autocorrelation function of $ x(t) $, $ r_{xx}(l) $, that can be writen in function of the autocorrelation function of $ s(t) $, $ r_{ss}(l) $. I was able to compute this autocorrelation function and I got:

$ r_{xx}(l) = (1 + \mu^{2})r_{ss}(l) + \mu r_{ss}(l-k) + \mu r_{ss}(l+k) $

So, it obvious has three peaks: when $ l = 0 $, when $ l = k $ and when $ l = -k $.

Assuming that $ r_{ss}(l) = 0 $ for $ |l| > k/10 $, I think that I can obtain the value of k observing the side peaks, and the value of $ \mu $ observing the central peak, but i wasn't able to formalize it mathematically. Could anyone help me to do so?

Thanks in advance.


texts: What is benefit when convert signed 16 bit to unsigned 8 bit in medical image?

texts: I have an image with its type is 16 bit signed. I see some paper, the author often convert it to unsigned 8 bit. But they did not explain why they do it. Could you explain what is benefit of this work? And How to implement it by matlab code?


texts: Segment highly packed nuclei with dotted high intensity signal inside each object

texts:
I collected high magnification images on cell nuclei and I need to segment them (red circle in the image). The edges of the nuclei are visible but inside each nuclear object there are bright spots of positive signal (Not an imaging artefact, if I will retake the images it will be the same).

I tested different way to segment the image (using python scripts based on scikit-image, scipy.ndimage and opencv):
- Threshold->morphology (close/opening)-> watershed (after gaussian filtering or median filtering to try to remove the dots)
- Graph-cuts (Al-Kofahi algorithm implemented in the Farsight package: https://github.com/gbkedar/farsight)

All the methods lead to oversegmentation because of the internal bright spots.

To solve that I thought was a good idea to isolate the dots with a laplacian of gaussian and used the opencv function cv.inpaint to fill the dots with the intensity of the surrounded pixel values. The result was not good either.

The questions are:
- Do you think that the best strategy is to isolate the bright dots and then try to replace them with the surrounding signal?
- Any strategy that will lead to a better result?

Any help will be greatly appreciated


texts: Why is the output of the LPF half that of the real/imaginary part of the complex envelope?

texts:
The passband signal $\ u_p(t) =
I_{[−1,1]}(t) \cos 400πt $ is passed through an LPF as $\ u_p(t) \cos(401\pi t).$ The question is to find the output. The example does the following:

If we let $\ u_1 = u_{c1} +ju_{s1}  $ denote the complex envelope with respect to the reference $\ e^{j401t} $,  the output of the LPF is $\ u_{c1}/2.$ Why is it half?

The passband signal can be rewritten as,
$\ u_p(t) = \Re(I_{[−1,1]}(t)e^{−jt}e^{j401t})$ to find $\ u_{c1}$
but how is this justified if we want to find the output of $\ u_p(t) \cos(401\pi t) $ and not just $\ u_p $ that is just rearranged.  Shouldn't the translation induced by the cosine be taken into account. I'd appreciate any help.


texts: Quick question about difference between two implementations of MFCC

texts: I'm trying to learn more about voice recognition, and I was looking this library here and noticed there were two implementations of MFCC (one called mfcc and another called mfcc2). I ran both of these on the same data and got different results, and was wondering why there are two different version of MFCC that give two different results. Thanks.


texts: How to extract information about each constituent 'wave' of a signal?

texts:
I have a vector containing wind speeds sampled at 1Hz over ten minutes.

I have applied a lowpass filter of 0.25Hz to get only the wind speeds I am interested in.

I now want to use MATLAB to identify information on all of the 'waves' which combine to form the complete ten minute wind speed profile. I'd like to know the waves' minimum wind speed, mean wind speed and frequency/time period.

I have tried using cross-correlation, however I couldn't get this to show me whether any of the test waves I'd applied were found in the original signal.

I will attach a figure containing filtered wind speed (600s), the 'test wave' of duration 9 seconds, and the result of the cross-correlation. Any test wave used results in same output of cross-correlation whereas id like to input various test waves and find out if they are present in the original signal. Maybe resampling?


texts:
So, suppose you have a superposition of cosines of the same magnitude:
$\sum_{i}\cos[\omega_i t]$

When you take the Fourier transform, theoretically you'd expect to find $\delta$-like peaks of the same magnitude. But in practice, what I find is this:

How would you explain this variation in the magnitude? Given such variation, would it be possible to recover the amplitudes of each frequency component in a superposition $\sum_{i}A_i\cos[\omega_i t]$

EDIT: I suspect it has to do with the relationship between the sampling frequency and the frequencies of the components. The fact that the theoretical $\delta$ response exists is a consequence of the orthogonality of the FT kernel. So, I suppose that in the discrete case, the frequencies of the DFT lines don't exactly "match" the sampled frequencies of the signal. I don't know how to articulate this properly or if it's even correct though.

For reference, here's the MATLAB code I used to generate this

Fs = 100;
Sp = 1/Fs;
L = 2000;
t = (0:L-1)*Sp;

y2 = cos(5.*t)+cos(25.*t)+cos(10.*t)+cos(1.*t)+cos(100.*t)+cos(200.*t)+cos(300.*t);

Y = fft(y2);
P2 = abs(Y/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

figure
plot(f,P1,'-k','Linewidth',1);


texts: DFT sample point k < N has negative frequency

texts:
(From: Schaum's DSP outline, 2nd edition, page 254, problem 6.35)

A signal $x_a(t)$ that is bandlimited to 10 kHz is sampled with a sampling frequency of $f_s = 20$ kHz.  The DFT of N=1000 samples of x[n] is then computed, that is:

$$
X[k]=\sum_{n=0}^{N-1} x[n] e^{-j(2\pi / N)nk}
$$

with N=1000.

(a) to what analog frequency does the index … k=150, and k=800 correspond?

A few equations:

$$
\omega_k = \frac{2 \pi k}{N}
$$

$$
f_k = \frac{\omega_k f_s}{2 \pi}
$$

For k=150:

$$
\omega_k = \frac{2 \pi}{1000} 150 = \frac{3}{10}\pi
$$

$$
f_k= \frac{3\pi 20000}{2 \pi 10} = 3000\ \  hz
$$

Ok no problem so far:

For k=300:

$$
\omega_k = \frac{2 \pi}{1000} 800 = \frac{4}{5}\pi
$$

Here's where i hit a problem:

book says,  for k=800, we need to be careful. Because $X(e^{j\omega})$ is periodic:

$$X(e^{j\omega}) = X(e^{j\omega + 2\pi})
$$

Therefore, k=800 corresponds to the frequency:

$$
\omega_k = \frac{2\pi}{N}k = \frac{2\pi}{N}(k - N) = -200\frac{2\pi}{N}
$$

My question is this:

what's the reason we need to use a negative frequency here?  Since k=800 < N=1000?  Wouldn't it just be:

$$
\omega_k = \frac{2 \pi 800}{1000} = \frac{8\pi}{5} = 5.02...
$$

Is this just the book being strange or am i missing some understanding of DFT?  I would have thought you wouldn't need to worry about $2\pi$ unless k > 1000?


texts:
I have two signals with nearly the same frequency, and i would like to find the phase shift between them, since i am using micro-controller, i would like to do it in the less expensive way.

I know how to calculate Goertzel for a specific frequency bin (for each signal), but how do i calculate the angle for this bin by using Goertzel?

If you know a lower computational method, please advice.

Thanks


texts: How to find smoothed estimates of the derivative and second derivative of a signal?

texts:
I have a signal sampled at $\Delta t$: $f_i(t_i=i\Delta t)$ where $i = 0,\ldots,n-1$.
I want to find the first and second derivative of the signal: $f'(t)$ and $f''(t)$.

My first thought was to estimate the derivatives by central differences:

\begin{align}
f'(t_{i})&=\frac{f(t_{i+1})-f(t_{i-1})}{2\Delta t}\\
f''(t_{i})&=\frac{f(t_{i+1})-2f(t_{i})+f(t_{i-1})}{(\Delta t)^2}
\end{align}

However the signal may have a lot of high frequency noise that may cause quick fluctuations in $f'$ and $f''$.

What would be the best way to find "smoothed" estimates of $f'$ and $f''$?


texts:
I'm writing a gunshot detector software. I'm using this article as reference:https://arxiv.org/ftp/arxiv/papers/1706/1706.08759.pdf.

For now, my software is able to recognize an impulsive sound. Now, I need to take this impulsive sound and check if it's and gunshot or not.

For gunshot detection, the author of this article uses MFCC features, but, I'm having a lot of difficulty to extract these MFCC features.

I'm using the Sphinx 4 Java library for this. Here is the base code:

// Windown size is 99 (16 KHz sample rate)
private void calculateMFCC(double[] windowAfterFFT) {
    final ArrayList<DataProcessor> pipeline = new ArrayList<>();
    pipeline.add(new SingleValueAudioSource(windowAfterFFT));
    pipeline.add(new MelFrequencyFilterBank(130, 6800, 1)); // params: (minFreq, maxFreq, numberFilters)
    pipeline.add(new DiscreteCosineTransform2(4, 20)); // params: (numberMelFilters, cepstrumSize)
    FrontEnd f = new FrontEnd(pipeline);
    log.info("Result: " + f.getData());
}

My mainly issue is that I can't figure out these items:

How do I define the minimum frequency, maximum frequency and the number of filters that I need for a gunshot detection on the mel frequency filter bank?
How do I define the cepstrum size that I need for gunshot detection?

Most of the articles that I read are not very detailed about these topics.


texts: Quantized error in scatter plots for a hamming windowed LTE signal

texts:
I have a time domain signal based on LTE standards.
Carrier spacing = 15KHz, Symbols per slot = 6, extended Cyclic Prefix = 512 samples
FFT size = 2048, No. of sub carriers = 1200, channel BW = 20MHz, Samp Freq = 2048*15KHz = 30.72 MHz, Total Bits = 14400 (Symbols * slots * bits used to MOD * total subcarriers).

It is QPSK modulated (using qammod ("IntegerSymbols",2^bits,binary) in matlab);
This QAM modulated data is IFFT'ed after zero padding on the sides and DC carrier insertion is done and then CP is added (512 samples to a 2048 sample bin).

This signal now is windowed (Point wise multiplied by a window function(only in the TRANSMITTER side) (FOR DIFFERENT REASONS) in time domain) "BY A HAMMING WINDOW FUNTION" Keeping all the overlap add constraint in mind.

On the receiver side, no noise is added, no specific channel is modeled. Simple reception. On the receiver side, CP is removed and then the stream is FFT'ed and only the central info carrying subcarriers are taken. This stream at this point is still modulated. However when I try to plot a "scatter plot figure" of the QAM modulated data from the txer side and the FFT'ed data in the receiver, I get a pattern as shown in the figure above which seems to be some sort of quantized thing. I am not sure what it is or what is it caused by. Could you please help me on it?
On decoding Error/BER is 0 but EVM is -18dB.

P.S: I've asked this ques before but I reckon the details weren't clear. So fingers crossed this time.


texts: Comparing spectrogram similarity over time

texts:
I have a data set of spectrograms from moving audio recorders. As the recorders get closer, I expect the spectrograms to become more "similar". I expect this "similarity" to be a function of distance and frequency.
I would like to create a 3d plot or image that has distance and frequency on the x and y, and "similarity" on z. I am not sure how to define this similarity. Is there a standardized way of comparing similar signals in this way?


texts: Principal Component Analysis (PCA) on Convolutional Neural Network (CNN) Features

texts:
Please, I have a question regarding PCA and features which are extracted from a convolutional layer. link

if we have a test dataset , and we extract all conv features of all images at test dataset called feat. Then we do the following normalization and PCA to all conv features

feats = normalize(feats)
pca = PCA(512 ,whiten=True)
pca.fit(feats)

And then use pca model for test dataset

# PCA MODEL - use paris for oxford data and vice versa
self.pca = pickle.load(open('_oxford.pkl', 'rb'))

And then applying PCA on the test dataset using the previous pca

print "Applying PCA"
self.db_feats = normalize(self.db_feats)
if self.pooling is 'sum':
self.db_feats = self.pca.transform(self.db_feats)
self.db_feats = normalize(self.db_feats)

What is the idea of using pca model of another dataset to transform the features of test dataset?


texts: What is the difference between signal processing and data processing?

texts:
I often find the terms 'signal processing' and 'data processing', and sometimes 'information processing' in articles. But what are the differences between these notations. In my opinion, signal processing is more related to a physical background but data processing or information processing is only related to raw time series, irrelevant of the data sources. For example, in array signal processing, directional of arrival (doa) is signal processing since it is closely connected to wave propagation, but smoothing the doa results is data processing.
Is it correct? Hope you can help clarify these notations.


texts:
I'm working on a sound source localization project and am interested in using GCC-PHAT as a feature. I'm quite new to this transformation and I've been reading a few papers about it and I'm still struggling with understanding what it actually is. One of the papers I'm interested in contains the following:

The maximum possible time delay between any 2 microphones is
$τ = 0.2/340 = 0.5882ms$ where the sound speed is assume to be
$340m/s$. Suppose we are using a sampling rate of 16kHz, then the
maximum delay in samples is $n = 16000τ ≈ 21$. Hence, for each
microphone pair, only the 21 correlation coefficients near the center
contains useful information for DOA estimation.

First of all I'm not sure where $n = 21$ is coming from, seing as $0.2/340 * 16000 = 9.411$, but perhaps its simply an error on their part, but I was curious if this has something to do with GCC-PHAT. I've found several implementations and when I take a 10ms audio sample with a sample rate of 44.1kHz (441 samples) I get the following plot:

So this is where I start to have questions as to what the actual feature is. But more to the point about the quote above, what exactly are the coefficients that are being referenced? Are these simply the correlation values around the center? And by transforming my x-axis to milliseconds, it seems that the peak is around 6ms. So does this mean that the signal arrived at microphone-1 6ms after microphone 0? Also any literature on the topic would be greatly appreciated, as I've so far only really been able to find code implementations but I'm trying to get a more intuitive understanding of the thing itself. Any help is greatly appreciated.


texts:
I have to find a convolution of two signals

$h[n] = 0.5^nu[n]$

$x[n] = u[n]-u[n-3]$

the final sum, which is correct is:

$$\sum_{m=n-2}^n 0.5^mu[m] $$

note that i replaced n-k with m, that is
$ m = n-k $

So, in regards to parameter n, i have to decypher the formula.
The result is:

$$0, n<0$$
$$2(1-0.5^n*0.5), 0\le n \lt 2$$
$$1.75*0.5^n/0.5^2, n \ge 2$$

But, i don't know how to get those results.
Can you please explain.

The results are correct (maybe i just typed something wrong, if that's the case i'm sorry in advance)


texts:
Can someone here explain the meaning of,  64 point FFT Hamming window interpolated to 720-point sequences and its inverse?

720 point is from 0 to $2 \pi$ or (0 to 720 degree). I get the 64 point hamming window. What I don't get is the interpolation to 720-point sequences. In what order should it be calculated?

Link to get more info the math and technology its not a problem just the FFT and its inverse.

What I am looking for is the resultant plot or graph is fig 18 on page 43.

I add this to answer A_A and make it more clear.
And the fact he made my mind more clear about the (hamming window) typo!.
Skuld have been moving Hanning window.
As A_A pointed out i may have been confused why it was implemented.
yes and no is my answer to that clearly it shows the variations in the,
Intake pressure when valve is both open and closed that is one part of the "problem" i have.
The other "problem" witch is also related to the same issue a wave and pressure.
When the valve is closed the intake air will hit a wall and punchback and usual we say.
4 waves. One cycle of the intake stroke is total 720 degree.

In between the opening of the valve. Lets say engine rpm is 18000, 60 sec / rpm = 0.0033 sec.

i forgot to mention before ivo (intake valve open) ivc (intake valve closed) in between we have
180 crankshaft degree tdc(piston strait up) degree = 0 bdc(piston at bottom) degree = 180

and we can set camshaft lobe center to 108 after tdc

crankshaft ivo = 20 degree before tdc + 180 and ivc 50 degree after bdc.
here is also one important thing camshaft rotate 1/2 the speed of crankshaft.

camshaft opening is 250 degree time closed is 720-250 = 470 / 360 1.30 revolution = 0.0043 sec
Speed of sound in air at 30 degrees Celsius is 349.08 m/sec 0.0043 sec time * speed of sound wave length 1.05 meter
impossible to take care of in a car. So shorten to 1/4 the final length is then 0.37598825 meter
Normal referred to as quarter of a wave, witch is the last of the 4 waves.
When that arrives at the valve i want that to be a sinusoidal wave and know its frequency.
Why? it tells me howe fast the intake valve needs to move,
along with what is pointed out in 3.3 The acoustics of the intake system.
The pressure is not constant during the intake stroke, also important to know.
First time i have seen fft used in context of intake flow.
So if any other suggestion i will be happy to listen and also code it.

Hope you can help


texts: Why does my calculated SNR have so much variance?

texts:
I am running a SNR test for a real system as part of a final system test. There is real noise. The noise is calculated as the ratio between two measurements.
I noticed however that my SNR value is actually quite noisy and all over the place. For example, I have a minimum target of $8000$ for my SNR, but after running the SNR test several hundred times I get a messy distribution of SNR values, like $7500$ with a standard deviation of $600$.


texts:
This signal represents a drawn line. For some reason, the device showed interference as shown in the figure below:

The question is: how can I remove the interference, making the signal as continuous as possible?

Any suggestions?


texts: Is sift a good way to extract features from an image?

texts:
I am trying to extract features from an image, but I have failed to get the points that I want to extract, and hence my image fails to match with the template.

Are there any checkpoints that I have to go through before applying SIFT to my image so that I get good results? For example the template image is;

targeted image;


texts: What determines the accuracy of the phase result in a DFT bin?

texts:
What are the factors that affect the accuracy and precision for the phase that is given by the DFT?

Just thinking medium-hard about this, it occurs to me that it must have something to do with the window length $L$ and the sampling frequency $F_s$, and maybe the fundamental frequency that I'm interested in extracting?

I'm going to keep working on some of the math, but I'd like to start a discussion and get some insight for which way to look.


texts:
I'm pretty new to signal processing and have encountered strange behaviour.

I have a vibration signal from a non-stationary gearbox component that I need to resample to do order analysis. The original sampling frequency is 25.6kHz. It seems to be working fine for most signals but one output looks weird.

The image below shows what the FFT transformation of a signal should look like.

This image shows the weird transformed signal.

Although the FFT looks similar, the axis are wrong. Does anybody know what could cause this? All help appreciated!


texts: What is used to initialise difference equations that require data from the output?

texts:
Given some difference equation for a filter like

$y[k] = ax[k] + bx[k-1]+ cy[k-1]$,

how would you initialize it? Since it needs an old output value (feedback) to calculate the new output value, it seems like the equation will keep on chasing its own tail. I would be tempted to just let

$y[k-1] = x[k-1]$

initially to calculate a value for $y[k]$ after which we could just use that for $y[k-1]$. There might be a better technique for this problem or even a really obvious solution to it, but I am at a loss on how. Most of the digital filter books/papers that I have read have done a great job on explaining how the filter gets derived, but the practical implementation requirements seem to be absent.

This is also my first time posting on Stack Exchange, so any comments on etiquette or style would be appreciated. Just be gentle ;)


texts: Understanding Pitch Detection with Autocorrelation Methods

texts:
I've been reading through A Smarter Way to Find Pitch, describing its pitch detection algorithm using autocorrelation.

I've been having trouble understanding the accuracy claims. It says:

  MPM runs in real time with a standard 44.1 kHz sampling rate. It operates without using low-pass filtering so
  it can work on sound with high harmonic frequencies such
  as a violin and it can display pitch changes of one cent reliably.

I cannot fathom how detecting changes of one cent is possible. Considering a note played at 441Hz, sampled at 44.1kHz, one would measure a first peak at tao=100 samples (delay space).

The next closest note which could be measured, by my understanding, would have a peak at the 101st sample. This would be a frequency of 437Hz, or about 16 cents lower than the first note. What am I missing here?


texts: I have an image (that is a matrix), let's say of dimensions NxN. I then want to expand this matrix into M basis matrices (for the moment I'm still unsure how many M of these basis matrices I should assume) of the same dimensions NxN, that is analog to basis expansion of a vector into basis vectors. But as for the base matrices, I want them to be either only the translation (fixed scaling) or only the scaling (fixed translation) of a given mother basis matrix. Is such an expansion justifiable, or probably even feasible in the first place? I'm interested in searching for a set of basis matrices whose corresponding expansion coefficients form a sparse vector.


texts: Converting floating point number to Q-notation fixed point number in C/Java?

texts: I am trying to convert from a floating point number to a Q-notation signed 16-bit fixed point number (e.g. Q15). How can this be done in C or Java? I searched online but could not find any available code.


texts: What Is the Transfer Function of a Moving Average (FIR Filter)?

texts:
To make post-processing easier, I export scope measurements as CSV files, which are then post-processed (mostly in Microsoft Excel, which is not the best tool for the job, but it is all I have at my disposal).

One of the post processing used is a simple moving average, and I was wondering how it does transform the source measurement. I expect it to behave similarly to a low-pass filter, but I guess it is not the same. Hence, my question.

I am specifically looking for a Bode plot. Thanks in advance.


texts: How to distinguish between received quadratures

texts:
How can you distinguish between two quadratures in QAM when the phase is not known at the receiver? For example, if the quadratures are sin and cos, how can the receiver tell which one is which, since both look the same except for the phases?

I've found the idea of random-phase carrier but I'm not sure if this is related to this question, or how it can be used if it is.


texts: Sub nyquist sampling, required number of samples for time sparse grouped signals

texts:
Question: Does it make sense to perform compressed sampling if the non zero samples are grouped in time? If so, what is the minimal length of the vector x that should be acquired to allow full signal reconstruction (of at least the non zero pulse segments)?

$y=\phi x$, where y is the compressed vector and x is the signal (sampled at nyquist)

Situation: time sparse signal (+- 97 % of the samples at noise level), but the 3% non zero samples are grouped (periodic). I'm only interested in the pulse samples (so not in the noise floor samples).

Thanks a lot!


texts: Implications of $X( -j\omega ) = X^*(j\omega)$

texts:
What are the implications of:

  If $x(t)$ is real and $x(-t) = x^*(t)$, then $X(-j\omega) = X^*(j\omega)$ and $X(j\omega)$ is real.

I am trying to understand it and I would like to research it further (what's the theorem that guarantees this?).

I'd like to understand why it's important aside from the fact that it guarantees a symmetric Fourier transform or frequency response if we're talking about that of the impulse response.

Thanks,


texts:
While there are plenty of methods to detect R peaks it appears that detecting S peaks is less discussed. Is there a recommended way of detecting S peaks?

Additionally, could you recommend an implementation in python?


texts: Root MUSIC and solution of conjugate palindromic polynomials

texts:
The root-MUSIC algorithm gives direction of arrival angles directly and avoids the computational overhead of computing the MUSIC estimator curve (see Direction of Arrival Estimation using a Root-MUSIC Algorithm in the Proceedings of the International Multi Conference of Engineers and Computer Scientists 2008 Vol II).

However, for a $P$ element antenna, the root-MUSIC algorithm constructs a conjugate palindromic polynomial  of degree $2(P-1)$ from the diagonals of a Hermitian matrix constructed from the noise space eigenvectors. The $2(P-1)$ roots occur in conjugate reciprocal pairs $z$ and $\frac{1}{z^*}$ by nature of the conjugate palindromic symmetry of the polynomial and their arguments give candidate DOA angles.

I'm solving this polynomial (of degree $2(P-1)$) by construction of the companion matrix and then using a general complex LAPACK eigensolver. This is obviously inefficient for two reasons:

the companion matrix has size $2(P-1) \times 2(P-1)$ whereas the problem is fundamentally of order $(P-1)$ and
I'm not able to use the constraint that the polynomial has conjugate palindromic symmetric leading to reciprocal pair roots.

Does anyone know of

a modification to the root MUSIC algorithm that constructs a polynomial of only degree $(P-1)$ containing just the roots inside or outside the unit circle or
having constructed the conventional root MUSIC polynomial of degree $2(P-1)$ is there an elegant way of factorizing it into two polynomials each of degree $(P-1)$: $$p(z)=z^{P-1}.q(z).q^*(\frac{1}{z^*})$$?


texts: Discrete-time Fourier transform of $a^{|n|} u[n]$

texts:
I have a problem calculating the DTFT of this pair:
Could anyone tell me why the DTFT for $a^{|n|} u[n]$ is different from $a^{n} u[n]$'s?


texts: FFT: do I need to include gaps in data in my FFT domain?

texts:
I am currently doing an FFT on a very large amount of vibration data, which I am going to sample with no pattern which I want to PSD eventually

My question is in essence can I go from my data

[1, 3, 5, 10, 2, 4, 3 , 23, 24]

then sample it to

A = [0, 0, 5, 10, 0, 0, 0, 23, 24]

If I fft(A),  will it give me the same answer as if I fft(B), where

B = [5, 10, 23, 24]

I think this is fine because if A[n] = 0 it provides no contribution to the sum in the DFT but this feels very wrong and I don't know anywhere near enough about this stuff to feel confident.

I can't access data set A because the way my data is sampled can only give me B.

More detail:

I have very large time series data for vibrations caused by air flow over an object, I am trying to parse that data to identify what happens to the vibrations at specific wind speeds and wind directions (this is the sampling I am talking about, and why its sporadic) .

However the data comes in a strange file format (.gcf) which I can only deal with using obspy, and the only way I can find to sample the data (because there is so much and my computer is rubbish) is using obspy's trim() function, which essentially deletes the rest of the not useful data. (I know this is stupid but I cant find a better way) And so this is what is reducing the dimension of the data .

And the end goal is to create a PSD plot but the problem is the fft really.


texts: question about the frequency of Gabor filters

texts:
The Fourier transform of 2D Gabor filters is shown below:

F is the radial frequency of the Gabor functions.
Another parameter Fm is denoted as highest central frequency.
Fr is the constant frequency ratio.

My question is : 1. central frequency is radial frequency?
                 2. so , Fm is the largest value among values of F?
                 3. Is F a series with common ratio Fr ?

Thank you very much for any help.


texts: How to calculate autocorrelation function of an image noise

texts:
I am trying to calculate the autocorrelation function of an image noise. Specifically, I want to calculate the NCORR parameter in the function deconvwnr(). I have knowledge about the noise, which is gaussian with mean 0 and variance 1. I found the examples in this website calculate the NCORR as follows:

NP = abs(fftn(noise)).^2;
NPOW = sum(NP(:))/prod(size(noise)); % noise power
NCORR = fftshift(real(ifftn(NP))); % noise ACF, centered

I tried following this method, but to do that, I first need to create an image of certain size and then add the gaussian noise into it. But the problem is that if I initialize the pixel values to be all 0, then after the application of gaussian noise, the image elements will only have intensity values above 0, which does not reflect the shape of the gaussian noise. On the other hand, if I initialize the pixel values to be all 1, then the image elements will only have intensity values below 1, which also does not reflect the shape of gaussian. Another way is to initialize the values to be 0.5, but then I was afraid that 0.5 will affect the NCORR values because it increases the noise power. Can someone show me how to do this or is there another way to calculate the autocorrelation function?

I = zeros([m, n]);
I = imnoise(I, 'gaussian', 0, 1);
# then follow the steps above


texts: Using Python with scipy.signal firwin, how do I know how many taps/coefficients should be used?

texts:
I need to alter an FIR filter which I know very little about. I have a signal which need to be altered, and looking at a Python script used for this task, I see that it uses Scipy's signal library, more specifically the firwin function.

If I input

cutoff_freq = 5400
num_taps = 51
coeffs = signal.firwin(numtaps = num_taps, cutoff = cutoff_freq, window = "hamming", nyq = 72000)

I get

  scipy\signal\filter_design.py:338: BadCoefficients: Badly conditioned
  filter coefficients (numerator): the results may be meaningless

The same occurs when num_taps is 52, 53, 54 and 55, but 56, 57, 58 and 59 does not give this warning.
What does this mean and how should I be sure to avoid it?


texts: Hamming window for filtering unwanted terms

texts:
I am trying to filter off two unwanted terms in this frequency spectrum (the spectrum on top).
First, I used a convolution between a rectangular window to make the unwanted terms equal to zero.

the first problem is when I tried to centralize the remaining part and multiplying it by a cos ( 2 * pi * f *t) I didn't get any result.

And my second question about hamming window is that in an article about a similar issue, a hamming and a blackman window were used. So, is it possible to use hamming window to filter off the unwanted terms? if yes, then how do I choose the starting and the final point of the hamming window to leave only the desired part! I really need help to understand how does it work and how to use it.


texts: How can I transform a time signal into a third octave spectrum?

texts:
I try to make a third-octave spectrum from a time signal and I have some troubles with the algorithm. I'm not an expert in signal processing so I don't know what is wrong.

Here what my algorithm looks like :

Time signal -> Z transform of the signal -> First octave band filtering -> Quadratic sum of the result's values => First octave band value.
Under factor 2 sampling ->Z transform of the signal -> Second octave band filtering -> Quadratic sum of the result's values =>Second octave band.

Now I have this, I try to test it on Matlab with a time signal about 2 millions points.

I have this code to calculate all the octave's bands :

    fc = fMin*2.^((0:bmax)*bw);             %Center frequencies.
    fl = fc*2^(-bw/2);                      %Lower cutoffs.
    fu = fc*2^(+bw/2);                      %Upper cutoffs.

And this one to calculate the Z transform of the time signal :

    x = [1 2 3];
    syms z;
    d = z.^(-1*(0:numel(x)-1));
    zTransform = sum(x.*d)

My main problem is, I'm not completely sure that my algorithm is right and if it is, how can I filter the z transform equation that the code above gives me with each third-octave band filters?

Thanks in advance, it's been a while I'm on this problem.


texts: I'm trying to create and process an OBPSK signal, which, as I understand it, transmits even bits on one axis and the odd bits on the opposite axis, resulting in a constellation with 4 unique points. OBPSK seems to be somewhat rare, and a lot of software doesn't have functions to handle it. I created it (if I understand it correctly) by modulating the even and odd bits as BPSK separately, and rotating one of them by 90 degrees. However, with 4 constellation points, it seems similar to QPSK, or more specifically OQPSK, which is common with software having built in functions for modulation, etc. Is there some way to process OPSK, perhaps with some sort of transformation and/or change in data rate, as OQPSK?


texts: Comparing two sound files in embedded-systems (C language)

texts:
what would be the most efficient ("cheap" in terms of system resources etc.) way to compare 2 sound files?

Here is the scenario. Sound A comes as an input. There is sound B stored in memory. We need to see if sound A is identical to sound B. We know that both sounds can be copies of the same sound file, but there is a chance that sound B was processed to reduce the size.

What would be the most efficient ways to compare these sounds and see if they are similar (or rather the same).

ADDED

"Cheap", architecture is Cortex ARM M0. 16kb ram, 256kb flash.
"Identical" means that these 2 sounds come from the same sound file
(same exact sound effect). THe difference may be that sound B (stored
on the microcontroller side) may be processed to reduced size.
Basically imagine we have 10 sound files, we copy N7 to MCU (maybe
process it to reduce size) and then play them and expect MCU to
recognize the sound file when it's played.

The kind of sound - technically it can be any sound (within the frequencies recognizable by humans). Sound effects (explosions, weapons firing, tires screeching), people saying words, melodies played by musical instruments. But the common thing is that they are all identical (instruments playing melodies would be same instruments same room, same time, recorded with the same microphone etc.).

Worth mentioning (to avoid confusion), that the sound is being sent to MCU from a sound output, not as a file.


texts:
I need FIR filter lengths in the order of 1e4 and above to obtain reasonable accuracy in desired frequency response. The problem is that optimisation in MATLAB (e.g. fircls or the Optimization Toolbox) is very slow (because the problem size becomes a bit silly). Does any software exist out there that solve long FIR filter synthesis problems quickly-ish (or quicker than MATLAB)?

Or are there better ways of generating super-narrow pass band filters (suitable for FPGAs)?

Sampling will be at about 1 GHz, output should be about 10 kHz, so 1/100000 pass to stop ratio. The FIR filter is to perform weighted averaging.


texts: Am I supposed to normalize FFT in Python?

texts:
I'm sorry if this question is posed often, I can't seem to understand the answers already out there. I'm working with several forms of the Fourier transform, including the FFT, PSD, and spectrograms. I'm not sure whether, upon computing the Fourier transform of my signal, I'm supposed to normalize the result by some factor. For context I don't have any interest in comparing signals, I just need my outputs to be mathematically accurate.

I've seen some answers that say the DFT has to be normalized by a 1/N factor - I'm interpreting this as the 1/N factor in the IFFT, is that correct? Or is there some other reason you would have to normalize by 1/N?

What's really confusing me though, is answers where I'm told to divide by the sampling frequency. This answer implies that it is equivalent to dividing by N- does this have something to do with the transform being computed over 1-second segments or something along those lines?

A tutorial I'm working from normalizes this way, and I'm not sure where it comes from. They construct a whitening procedure like so:

def whiten(strain, interp_psd, dt):
    Nt = len(strain)
    freqs = np.fft.rfftfreq(Nt, dt)

    # whitening: transform to freq domain, divide by asd, then transform back,
    # taking care to get normalization right.
    hf = np.fft.rfft(strain)
    norm = 1./np.sqrt(1./(dt*2))
    white_hf = hf / np.sqrt(interp_psd(freqs)) * norm
    white_ht = np.fft.irfft(white_hf, n=Nt)
    return white_ht

where interp_psd is a power spectral density interpolant of the noise, and is computed in the following way:

NFFT = 4*fs
Pxx_H1, freqs = mlab.psd(strain_H1, Fs = fs, NFFT = NFFT)
Pxx_L1, freqs = mlab.psd(strain_L1, Fs = fs, NFFT = NFFT)

# We will use interpolations of the ASDs computed above for whitening:
psd_H1 = interp1d(freqs, Pxx_H1)
psd_L1 = interp1d(freqs, Pxx_L1)

Why do they multiply the whitened strain in the frequency domain by dt*2? I understand that 2 could come from discarding negative frequencies, but wouldn't the RFFT take that into account?

The tutorial later on computes the spectrogram of the filtered data:

# Plot the H1 spectrogram:
plt.figure(figsize=(10,6))
spec_H1, freqs, bins, im = plt.specgram(strain_H1[indxt], NFFT=NFFT, Fs=fs, window=window,
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
plt.xlabel('time (s) since '+str(tevent))
plt.ylabel('Frequency (Hz)')
plt.colorbar()
plt.axis([-deltat, deltat, 0, 2000])
plt.title('aLIGO H1 strain data near '+eventname)
plt.savefig(eventname+'_H1_spectrogram.'+plottype)

Here, they don't normalize, but I recognize that they might not need to since they're only interested in the visual of the spectrogram.

Later they normalize by the sampling frequency when performing a matched-filter exercise, but then reverse it.

# Take the Fourier Transform (FFT) of the data and the template (with dwindow)
    data_fft = np.fft.fft(data*dwindow) / fs

    # -- Interpolate to get the PSD values at the needed frequencies
    power_vec = np.interp(np.abs(datafreq), freqs, data_psd)

    # -- Calculate the matched filter output in the time domain:
    # Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin.
    # Taking the Inverse Fourier Transform (IFFT) of the filter output puts it back in the time domain,
    # so the result will be plotted as a function of time off-set between the template and the data:
    optimal = data_fft * template_fft.conjugate() / power_vec
    optimal_time = 2*np.fft.ifft(optimal)*fs

I apologize if this is too much information. All in all, my questions are these: do I have to normalize the output of a FFT in python (numpy, scipy, matplotlib) in order to be mathematically accurate, and by what factor? And does that normalization differ for transforms like the PSD or the spectrogram?


texts:
I'm trying to distinguish eletrical engines by the magnetic field they create. I have no idea whether this is a smart idea, it's possible at all, whatsoever. However it brought me to a basic question about the Fourier transformation to which I couldn't find an answer that I understand :)

I am sampling with a magnetometer at 160Hz. I'm taking the 160 samples I get each second and calculate the DFT. When I plot the DFT over time I see pretty much a flat surface along all frequencies. The amplitudes vary depending of the motor or distance of the magnetometer but I am not seeing what I was hoping to see -  some characteristic frequencies of the motors to reliably distinguish them.

This got me wondering whether it is necessary to get rid of frequencies that cannot be detected before applying the fourier transformation. As I understand it, sampling with 160Hz allows me to detect frequencies up to 80Hz. Would it be necessary to use a lowpass filter and get rid of all frequencies below 80Hz? A frequency of e.g. 160Hz is a multiple of 80Hz, 40Hz,..., so wouldn't it "leak" into the 80Hz, 40Hz,... buckets?


texts: In a signal there is a peak, I want to find the start of the peak and the end of the peak

texts:
I want to the start time and end time for the peak shown below.
The start and end are shown by the marked circles in the fig
As there are number of such plots and want to know start and end time for each plot. So if its an automatic procdure thats good.
Can any one help me with this
Thanks


texts:
I am trying to figure out how to use the matlab gausswin function which constructs a Gaussian window for $N$ samples with a given standard deviation $\sigma$. The function is defined by

w = gausswin(N, alpha)

where alpha is $1/\sigma$. I would like my Gaussian window to contain $\pm 3 \sigma$ in the Gaussian. So for $N$ samples it seems I should do: $6 \sigma = N$ or $\sigma = N/6$. However, here are the results for $N = 100$ and $\sigma = 100/6 \approx 16.666$:

Clearly this Gaussian does not extend to $\pm 3 \sigma$. Can someone tell me what I'm doing wrong?

My guess it is has something to do with some length scale which isn't being properly specified. It seems to me that the alpha input to the gausswin function is some kind of normalized (unitless) standard deviation. Perhaps I'm not accounting for that.


texts: Getting an specific frequency component from sampled audio

texts:
I have a half second audio data sampled at $44.1\textrm{ kHz}$. If I FFT it I get frequency components in $2\textrm{ Hz}$ increments. Now suppose I want the component of an arbitrary frequency, like $32.45\textrm{ Hz}$, how can I get it?

Looking at the DFT formula from Wikipedia:

$$
X_k = \sum_{n=0}^{N-1} x_n e^{-i 2 \pi k n / N}
$$

I wondered, what if I calculate the relevant $k$ and plug it into the equation? Like:

$$
{32.45 \times 22050 \over 44100} = 16.225\quad \textrm{then}\quad X_{16.225} = \sum_{n=0}^{N-1} x_n e^{-i 2\pi \times 16.225 n / N}
$$

would that get me what I want?


texts:
In BPSK for demodulation at the receiver end we do phase recovery between received wave and carrier wave. We can find the phase difference between the carrier and received waves using :

PLL
And  Hilbert transform.

Phase locked loop recovers the phase difference after some iterations and becomes constant (when it gets a lock on exact phase difference).

This is the code for finding phase difference using Hilbert, the function take two waves as parameters and returns the phase difference between them.

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

We can find the phase difference using Hilbert function using whole waves(carrier and received) or phase difference at every bit( bit by bit )

My question is that:

In practice phase recovery/synchronization is done just at starting of two waves(carrier and received) or at every bit?

If i do the bit by bit phase synchronization:My demodulation result gives all zeros or ones here is why.

say my carrier wave: $S(t) = \cos( 2\pi f_c t )$ for bit 1, and $S'(t) = \cos( 2\pi f_c t + \pi )$ for bit 0

Lets say at the receiver end there is a phase difference ($\varphi$) at every bit, between received and carrier waves.

So recovered carrier looks like this $\cos(2\pi f_c t + \varphi)$ for 1 and $\cos(2\pi f_c t + \varphi + pi)$ for 0

Hence demodulation will be  $A_m(t)\cos(2\pi f_c t + \varphi)*A\cos(2\pi f_c t + \varphi) = A^2$. for bit 1

and $A_m(t)\cos(2\pi f_c t + \varphi + \pi)*A\cos(2\pi f_c t + \varphi + \pi) = A^2$ for bit 0.

Hence both bit zero and one give a positive value so how do we decide the decision boundary.

Should phase recovery be done only at the beginning of the received
wave?
If so how do we incorporate if there is a phase error with time?


texts: Are all exponential functions eigensignals of LTI systems?

texts:
I know that complex exponential functions are eigensignals to LTI systems. Do these include real exponential functions? E.g. $2^t, e^t, ...$
Thanks for the help!


texts: Can bandwidth of the closed loop be bigger than bandwidth of individual elements in it?

texts:
This is a question from control theory but I hope you can help me with it. I was pointed to this community.

So in my control loop I have elements that have cut-off frequency of 3-5kHz and I designed a PID controller for the inner loop and those elements and suddenly I get a closed loop bandwidth of 30-50kHz cut-off frequency.

So the loop is in the figure bellow.

So the PID controller, or voltage controller I designed is a PID with a derivative filter. Current sources are modeled as 1st order filters with 3-5kHz cutoff frequency. How is it possible to get such a high bandwidth of the closed loop if the inner elements have a lower cut-off frequency?

Thanks


texts:
We are currently implementing the velocity synchronous discrete Fourier transform (vsdft) for order
tracking rotating machinery in Python.

The formula is similar to dft. Simplified, it is given by

$$ \text{VSDFT}(\Omega) = C\sum_{n=0}^{N-1} x[n] \omega[n] e^{j\Omega \theta[n]} $$

where $C$ is a constant and

$$ \omega = \frac{d\theta}{dt} $$

Since dft is given by

$$ F_k = \sum_{n=0}^{N-1} f_n\cdot e^{-i 2 \pi k n / N} $$

and we wanted to use the optimized fft implementing of numpy, and was thinking if there is a way to write the input f so that we could
use fft directly. However, my math skills are failing me.

Can any of you see a way or any pointer to implement the above VSDFT without having to reimplement the Cooley-Tukey FFT algorithm?

Refrence to the original paper can be found in here.


texts: Magnification of a colored image through two color space: RGB & YCbCr

texts:
I have tried to magnify a colored image using bicubic method. I did this in two ways; First, I resize R, G, and B parts of it. Second, I transformed RGB to YCbCr color space and then resized Y, Cb, and Cr parts of it. When I use YCbCr color space, final image changes in color but not in whole shape of the image. I could not realize the reason. May you help me?

code:

clear all
clc

A=imread('Sun.jpg');
R=A(:,:,1);
G=A(:,:,2);
B=A(:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y  = 0.257.*R + 0.504.*G + 0.098.*B + 16;
Cb = -0.148.*R - 0.291.*G + 0.439.*B + 128;
Cr = 0.439.*R - 0.368.*G - 0.071.*B + 128;

Y1  = imresize(Y, 2, 'bicubic');
Cb1 = imresize(Cb, 2, 'bicubic');
Cr1 = imresize(Cr, 2, 'bicubic');

R1 = 1.164.*(Y1-16) + 1.596.*(Cr1-128);
G1 = 1.164.*(Y1-16) - 0.813.*(Cr1-128) - 0.391.*(Cb1-128);
B1 = 1.164.*(Y1-16) + 2.018.*(Cb1-128);

A1 = cat(3,R1,G1,B1);

imwrite(A1,'Sun1.jpg');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R2=imresize(R, 2, 'bicubic');
G2=imresize(G, 2, 'bicubic');
B2=imresize(B, 2, 'bicubic');

A2 = cat(3,R2,G2,B2);

imwrite(A2,'Sun2.jpg');

Original image:

Resized image(YCbCr color space):


texts:
Given the absolute value of the Fourier transform of a signal $g(t)$:
$|G(f)|$
If I compute the inverse Fourier transform of $|G(f)|$,
$$\int_{-\infty}^{+\infty} |G(f)|\,  e^{j2\pi ft}df$$
do I obtain $|g(t)|$?


texts: Finding angles from an IMU (Quaternions to Eulers)

texts:
Greetings and apologies if this has already been answered, but I haven't been able to get a clear answer from looking on Stack Exchange so-far.

I have Quaternions coming from two different IMU's located on the trunk and the back of the head, in order to describe relative rotations between both.

Of course if I would make use of the Euler transform for each of the IMU´s (and subtract to find the relative angles), gimbal locks manifest. But we SHOULDN'T have that problem if we were to represent the relative angles (which is the case, because the aim of this is to show an angular representation that a physiotherapist can understand), because relatives are quite small angles, and probably won't go near +-90.

So what I have in mind is taking the quaternion Q1 that represents the trunk IMU and the quaternion Q2 that represents the head, find the quaternion that rotates Q1 into Q2 (let´s call it Qr) and transform THAT ONE into Eulers.

And it... works. But for some reason I have a "crosstalk" between the pitch and the roll, meaning that they are NOT independent between them.

Is this procedure flawed? What would be a better approach to this?


texts: Need to extract phase of incoming signal, how to find FFT bin of frequency of interest in GNU Radio?

texts:
When receiving a higher frequency signal, say 2 GHz, the freq as shown in frequency sink is often off by some amount (1-2 MHz) on the ettus B200. So I can't just assume the location of the signal is at the halfway point in the FFT bin.
Right now I'm using a hacky method using argmax and probe blocks to find the correct bin from which to extract the phase. But this doesn't feel fundamentally sound and doesn't work as well as I'd like. For example if I tune the frequency mid-measurement things start to go wonky. Anyone have any advice? My screen is ridiculously small so I can't capture the whole screen at once, but here's the relevant portion of the flow graph

Thank you in advance
EDIT:
Problem was a misunderstanding in GRC. I'm fairly new to it so I didn't realize that the freq conversions happen at the USRP, not the waveform generators. The signal I was seeing on Frequency sink was VERY biased. I reformatted flowgraph to appropriate settings for the baseband frequency I was using, and the frequency of interest is exactly where it should be now,
fft_bin = (baseband_freq * fft_size / samp_rate) + fft_size/2
For those curious, this was just a pass-through sort of test comparing phase at USRP sink of Tx/Rx port with that of USRP source at Rx2 port. Serving as proving I understand what is happening before moving to next step.
EDIT2:
Leveled up: screen capture skills acquired


texts: FFT peak estimation: Zero-padding vs signal repetitioon

texts:
I am working on obtaining the vortex shedding frequencies behind a 2D square using CFD. Time-varying data is being extracted from a specific cell behind the square. I've written a little python script which performs FFT analysis on the signal and spits out the dominant peak frequency.

Because the shedding frequencies are so low (0.1 - 0.2 Hz) I added a zero-padding routine as well as a quadratic interpolator in order to have accuracy to several digits. My background is in fluid mechanics, not signal analysis. So my question is the following:

Would it be more beneficial to simply take my input signal (of say, 200 seconds of data) and repeat it several times (therefore increasing the length of time and subsequently the frequency resolution) or to keep the zero-padding technique?

Would the former method be much more computationally intensive? FYI, the signal is nearly sinusoidal.

Thanks


texts:
Let's say we have the following matrix:

0 0 0 0 0 0
0 1 0 0 0 0
1 1 0 0 0 0
1 1 0 0 0 0
0 0 0 0 0 0
0 0 0 1 1 1
0 0 0 0 1 1
0 0 0 0 0 0

I want to calculate the size of every connected region of 1s and also get the corresponding matrix elements of each region. The size should be computed as (width + heigth) / 2.

In the example the bottom left region would have size (2 + 3) / 2 and the bottom right region would have size (3 + 2) / 2.

Additional to the size of the regions I also want the corresponding entries of the regions, i.e. region 1 has entries (2,2), (3,1), (3,2) and so on.

How can this be done in Matlab?

I'm looking forward for the answers and wish everybody all the best.


texts: Which method to choose to get linear phase response zero group delay low-pass filter?

texts:
Background: I do time-dependent simulations of physical systems. Like this: http://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.115431 During the simulation, observables are recorded and later Fourier (usually sin-transformed) transformed to analyse things like absortion of a nanoparticle. The problem is that the system is huge in and all time steps cannot be stored to hard disk. But maybe I could save every 10th or 20th iteration to hard disk. This quantity is the density matrix. But for this questions purpose, each element of this matrix are to be filtered independently.

To prevent folding, one needs to low-pass filter the signal. Our analysis methods are extremely sensitive to phase, if the signal is sin, it needs to remain that way. Phase deviation (mixing of cos) will ruin the analysis. Then the filter has to have a zero delay (but it can access future data to some extent). Then it has to filter sharply before Nyquist frequency to prevent folding.

I played with matlab a year ago a bit, but was unable to produce such a filter. I have very litle experience with DSP, so I am therefore asking for suggestions. Is such a filter possible, what I require? To what extent one has to give up on the requirements? IIR or FIR? How to generate the filter coefficients? How long does the filter have to be?

Lets say the data is 10Hz, can I get a linear phase response with no delay up to 0.4Hz, if I resample a low-pass filtered signal at 1Hz.

I was also thinking that if this is not possible, perhaps one could send the preliminary filtered and presampled data through a "correction" filter to reconstruct the data at 10Hz.

(In reality, the signal is $10^{15}$ Hz)

Update: I tried to create symmetric FIR low-pass filter.

%symmetric FIR low-pass filter
N=200
M = N/2;
n=0:N;
wc=0.02 * 2 * pi;
h = sin(wc*(n-M)) ./ (pi*(n-M));
h(M+1) = wc/pi;
new_signal = conv(signal, h,'full');
plot(signal,'b-');
hold on
plot(new_signal,'r-');
% Resample
time=(0:(prod(size(new_signal))-1))*dt - M*dt;
time=time(1:10:end)';
new_signal = new_signal(1:10:end);

In the figure below with bold red, is a quantity calculated from the original signal. It has a frequency content till about 80 (eV, but unit is irrelevant). It has a characteristic feature that it is always positive. If I just resample the signal 10:1 without any filtering, I get the green curve. It has terrible folding effects, but it is actually not that bad. With blue is my first attempt of a symmetric filter, but I used conv(...,'same') so I basically tossed away the beginning of the signal. This resulted into some unwanted phase effects (negative values close to zero frequency). In purple curve, I keep the acausal part of the signal, and this resolves that issue (the conv(...,'full', and -M*dt parts of my script). One has also clearly a low-pass filter effect near frequency 10. Cyan curve the same thing, but only with 20 coefficient FIR filter, and that does not seem to work at all.

In the second figure one sees that basically all the filters cause unwanted negative signals close to zero frequency. So I have not yet found the filter I am looking for.

Also, none of the filters produce well the amplitude of the 'plasmon' peak at 4eV.

Any ideas?


texts:
I would like to compute how many zeros and ones (e.i. 255) in a binary image. The following code to generate a binary image by using Matlab.

I = imread('images.jpg');
level = graythresh(I);
BW = im2bw(I, level);
imshow(BW)
imwrite(BW, 'img.jpg');

This is my C++ code.

#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv\cv.h>
#include <map>
#include <iomanip>

std::map<int, int> computeHistogram(const cv::Mat& image)
{
    std::map<int, int> histogram;

    for ( int row = 0; row < image.rows; ++row)
    {
        for ( int col = 0; col < image.cols; ++col)
        {
          ++histogram[(int)image.at<uchar>(row, col)];

        }
    }

    return histogram;
}

void printHistogram(const std::map<int, int>& histogram)
{
    std::map<int, int>::const_iterator histogram_iter;
    std::cout << "\n------------------\n";
    for( histogram_iter = histogram.begin(); histogram_iter != histogram.end(); ++histogram_iter)
    {
      std::cout << std::setw(5) << histogram_iter->first <<  " : " << histogram_iter->second << "\n";
    }
    std::cout << "------------------\n";
}

int main(int argc, char **argv)
{
    cv::Mat img = cv::imread("img.jpg", CV_BGR2GRAY);
    printHistogram(computeHistogram(img));
    return 0;
}

This is what I got

------------------
    0 : 11451
    1 : 590
    2 : 602
    3 : 498
    4 : 428
    5 : 376
    6 : 387
    7 : 349
    8 : 314
    9 : 278
   10 : 229
   11 : 206
   12 : 212
   13 : 159
   14 : 142
   15 : 112
   16 : 106
   17 : 96
   18 : 86
   19 : 54
   20 : 44
   21 : 37
   22 : 26
   23 : 26
   24 : 14
   25 : 14
   26 : 10
   27 : 4
   28 : 7
   29 : 8
   30 : 2
   32 : 7
   33 : 2
   34 : 3
   36 : 2
   37 : 4
   38 : 1
  211 : 1
  216 : 1
  217 : 1
  219 : 4
  220 : 2
  221 : 2
  223 : 3
  224 : 5
  225 : 2
  226 : 4
  227 : 4
  228 : 5
  229 : 8
  230 : 15
  231 : 19
  232 : 22
  233 : 25
  234 : 36
  235 : 21
  236 : 54
  237 : 54
  238 : 64
  239 : 93
  240 : 113
  241 : 111
  242 : 127
  243 : 163
  244 : 193
  245 : 182
  246 : 202
  247 : 218
  248 : 263
  249 : 289
  250 : 291
  251 : 301
  252 : 403
  253 : 458
  254 : 469
  255 : 29211
------------------

Is this right? What I'm expecting is the number of 0s and 255s.


texts:
Have many audio files 1 minute long. Some of them are normal. Some of them are noise.

Here is a normal file:

This is a noise:

And this is noise too:

How do auto detect noise using python?

BTW: if you suggest any external console program which can do that it will be cool too.


texts:
I've been trying to solve the following problem, but I noticed that I get two different answers by means of two different techniques:

Find the impulse response, $h(n)$, of the following system, described in frequency:
$Y(e^{j\omega}) = e^{-j\omega}X(e^{j\omega}) + \frac{dX(e^{j\omega})}{d\omega}$

Try #1: I found the difference equation for the system and iterated while making $x(n) = \delta(n)$:

$y(n) = x(n-1)+\frac{n}{j}x(n)$

$y(0) = 0 + 0$

$y(1) = 0 + 1$

$y(2) = 0 + 0$

And so on. Which yields $h(n) = \delta(n-1)$.

I found this same result by noticing that $X(e^{j\omega}) = 1$ for $x(n) = \delta(n)$, which gives: $Y(e^{j\omega}) = e^{-j\omega}$ and by using the time shift property we find that $h(n) = \delta(n-1)$.

Try #2: I used the inversion formula for $Y(e^{j\omega}) = e^{-j\omega}$:

$y(n) = \frac{1}{2\pi}\int_{-\pi}^{\pi}e^{-j\omega}e^{j\omega n}d\omega$ which yields $y(n) = \frac{sin(\pi(n-1))}{\pi(n-1)}$

Did I miss something? What is the correct result?


texts: Length of padding for multiple samples stored in a single MP3 file

texts:
I have 20 short audio samples of a piano playing a single note. Due to internal limitations of my application, I need to concatenate and combine these into a single lossy file (either MP3 or AAC). This concatenated buffer will later be read, decoded to PCM, and split back into 20 buffers.

I realize that both VBR MP3 and AAC take advantage of various pyschoacoustic phenomena, including possible temporal masking. Is there a specific amount of silence that I must insert between the various samples in order to ensure that one does not affect the encoding of the other? Is this impossible to know due to the way encoders work?

I know that MP3 has a frame size of 1152 samples and AAC has a frame size of 1024 samples. Would 1152 samples of silence be enough?


texts: Amount of delay in audio signals which are echos of the original audio signal

texts: Could you please inform me, in simulations and also in practical cases how much is delay in audio signals which are echos of the original audio signal (I mean the reference assumption in the literature)? Thanks.


texts: What Is the MATLAB `imreconstruct()` Useful For?

texts: Could some one tell reasons to use the MATLAB imreconstruct function when processing images? I already have studied the topic and what the function actually does to the image. However, I cannot figure what we would use it for. What output should we expect from calling doing this function?


texts:
I am using a DAC to transmit 2 different voltage levels. High and low. At 1Gsps. There is no carrier wave or anything. These are just raw 0 and 1 time samples going across a wire.

On my ADC I am also sampling at 1Gsps and I'm able to fully reconstruct the transmitted signal. In other words I can tell which sample is a 0 or 1 level.

Does this not violate the Nyquist sampling theorem? I assume if I am transmitting at 1Gsps my receiver would need to sample at 2Gsps.

Where is my understanding wrong here?


texts: Fuse two sources of linear acceleration with a Kalman filter

texts:
How would I fuse two different sources of linear acceleration with a Kalman filter (perhaps linear acceleration readings from an IMU and from a dedicated accelerometer)?

My state is defined by distance and velocity; it is very similar to the following post, except for my availability of two sources of linear acceleration:

Using the Kalman filter given acceleration to estimate position and velocity

I'm not sure which acceleration reading I would use for my control vector ($u_{k}$) in the predict function.


texts: Filtering Corners Detected on Multiple Scales

texts:
I made a multi-scale Harris corner detector (inspired by MOPS) which is simply a Harris corner detector performed on several scales (i.e. subsampled versions of the same image).

This effectively make the detector more robust to noisy/blurred images and able to find more features.

However, there are cases where two or more features are detected at the same location on different scales. See the image - each circle represents an oriented corner, radius of the circle corresponds with its scale:

Because of this, some features have similar descriptors and this would theoretically weaken later correspondence matching.

Detectors like SIFT or Harris-Laplace deal with the problem by either detecting scale-space maxima only or finding characteristic scale of each detected corner.

I think it would be reasonable to group corners that are close to each other spatially and then choose just the strongest one (the one with the strongest response is nearest to its characteristic scale). I am already performing Adaptive Non-Maximal Suppression, but this happens on each scale separately (to keep moderate amount of features on each scale) so the corners with similar locations can still remain when they are on different scales.

Maybe just operating on descriptor vectors would be sufficient, but this problem have much higher dimensionality.

Any suggestions?


texts: How to get a non-equally spaced FFT back into the time domain

texts:
I have a signal that I STFT and then filter using an ERB spaced filterbank. At some point after this I want to get the signal back into the time domain, how can I go about this? Using a standard iSTFT function won't work because it assumed linearly spaced frequency bins, AFAIK? I've put a code snippet below.
I'm also not sure what what to tag this question as apart from fourier-transform
Y = stft(sig) # Y.shape = (1025,4000)
fb = filterbank() # fb.shape = (20,1025)

Y_erb = matrix_multiply(fb,Y) # Y_erb.shape = (20,4000)
```


texts: Correlation between frequency and initial phase estimates using FFT

texts:
Why is there a high negative correlation between frequency estimation errors and initial phase estimation errors when measured using FFT?
I have a simple code showing this negative correlation.
is there a known explanation for that?
here is the code

clear all;clc;close all;
%%
T = 0.2;
fs = 500;
f = 200;
ts = 1/fs;
snr = 5;
n = fs*T;
t = linspace(0,T-ts,n);
theta = 0.1;
y = exp(1i*2*pi*(f.*t + theta));
for count = 1:2000
    count;
    yy  =  awgn(y,snr,'measured');
    NFFT = 2^16;
    fVals = fs/NFFT.*[0:NFFT-1];
    X = fft(yy,NFFT)./n;
    plot(abs(X))
    [~,M] = findpeaks(abs(X),'MinPeakheight',0.8);
    measuredAngle = angle(X(M))/2/pi;
    errfreq(count) = fVals(M) - f;
    errph(count) = measuredAngle - theta;
end
cov(errfreq,errph)
corrcoef(errfreq,errph)

ans =
    1.0000   -0.8506
   -0.8506    1.0000


texts: Running a sox spectrogram on a bandpass filtered wav file

texts:
Let's apply a 400-2000Hz bandpass filter on a respiratory .wav file:

sox audio1.wav audio1.bandpass.wav sinc -t 10 400-2000

Now I'd like to generate a spectrogram that takes advantage of the reduced bandwidth to still generate a full size image. But how to to do that? Notice what the result is from running spectrogram on the bandpass file:

sox audio1.bandpass.wav -n spectrogram -r -o audio1.bandpass.png -m

That's clearly a big waste of image real estate. What step(s) am I missing here to use the full image size to focus on that already-limited frequency bands?

Update From accepted answer: works great

Add a step after the initial bandpass:

  sox audio1.bandpass.wav -r 4000  audio1.bandpass1.wav

Then

sox audio1.bandpass1.wav -n spectrogram -o audio1.bandpass.png -m


texts: Creating synthetic frequency domain box filter vs 4 pixels in corners

texts:
Creating two 10x10 black images.
1st image, adding a 2x2 white box in the middle (Low pass filter).
2nd image, adding 1 pixel in each of the corners.

Then performing the following steps in that order:
1. ifft2 (Transformation to the time domain)
2. ifftshift (Without it the results would be shifted, not sure why, and shift before the ifft2 would not help).
3. Normalize the results to emphasize the changes.
4. Show the results as images.

Why am I getting the same results out of the two images we created above?

I tried to explain my question the best I could, but here are images of the two frequency domain and then their result:


texts: How to design very narrow, sharp lowpass filters - Only DC needed?

texts:
I asked exact the same here one year ago: One year ago question I asked this question again here, because after one year, I have learned more about DSP and am thinking ways to improve it.

Filter description

Filter input data rate is 64kHz
Filter output data rate is 0.5Hz - 2Hz
The filter is preferably to have bandwidth of 1mHz with very sharp transition area, since only DC is needed and the noises at 0-5Hz is large.

The answer I got

The averaging data is the best DC I can get. I bought this answer. However, it now becomes plausible to me, since a filter always have a bandwidth and frequency responses. Averaging is essentially a FIR with no superior characteristics.

Questions

Shall the filter be designed in this way: 9 stages of halfband filters(64kHz -> 125Hz). Followed by three stages of cascaded integrator–comb (CIC) filters with R=5. Then the output will be fed into a 1st order IIR filter?
Maybe 16 stages of halfband filters? (64kHz ->  0.9765625Hz) And then be fed into an IIR filter.
How shall the filter be designed?


texts: HPS Algorithm for detecting the fundamental frequency of a guitar string

texts:
When you plucked on a string, the string is vibrate and except the fundamental frequency,
there are an integer multiple of the fundamental frequency, which are the harmonics.

I know that the HPS algorithm is a tool to find the fundamental frequency, but
I don't understand the downsample part. How can I know how many times I need to downsample
the signal?

I read somewhere that the fundamental frequency can be given by measuring the harmonic frequencies and their GCD (greatest common division)

another source says that there's a need to use the Hanning method.
can someone please explain this algorithm?


texts:
can you tell me how to write video data to file..
what is the alternate of these 4 commands in octave

videoreader('bg.mp4');
readframe(v);
videowriter('bg1.avi');
writevideo(newv,k);


texts: How can I calculate the frequency using different color bands?

texts:
Let's say that I have a 1 dimensional signal that has data which comes from one of 2 different spectral bands. The measurements come evenly mixed between the two bands. I'm interested in a method of extending the frequency calculation by combining information from both signals.

To show an example, let's just say that I have a Red band and a Green Band. The detector arrays would look like this: RGRGRGRGRG. Let's say that I had 512 R, and 512 G, just to give a concrete example. As the detectors are neighboring each other, and most signals have some component of red and green, it seems likely that a common frequency calculation could be made by combining the bands, and possibly double the frequency range which can be detected. How could I best go about doing this?


texts:
In the previous question,

Smoothing in frequency domain and impulse reponse of the filter in time-domain

I guess I have made a mistake in-terms for applying windowing to FFT.

After reading "Understanding Digital Signal Processing - Richard G. Lyons" few things are rather very clear now. Nevertheless to be certain I would like to ask whether my understanding is correct or not.

Frequency Domain Windowing is simply given by,

$$X_{\mathrm{three-term}}(m)=\alpha X(m) - \frac{\beta}{2}X(m-1)-\frac{\beta}{2}X(m+1)$$
where, $\alpha$ and $\beta$ are the co-efficients of Hamming window or Hann window.

Further
$$X_{\mathrm{three-term}}(0)=\alpha X(0)-\frac{\beta}{2}X(N-1)-\frac{\beta}{2}X(1)$$
and
$$X_{\mathrm{three-term}}(N-1)=\alpha X(N-1)-\frac{\beta}{2}X(N-2)-\frac{\beta}{2}X(0)$$

So after taking fft(X(m)), the windowed N-point FFT is simply the above three equations?

All these formulas are from Understanding Digital Signal Processing - Richard G. Lyons and all the credits goes to the author of the book.

Some questions

Is Windowed FFT called smoothing?
Do I need to do anything else apart from ifft(X_3term(m)) to perform inverse of $X_{\mathrm{three-term}}(m)$?


texts: A non-lossy aka lossless codec without patent-restrictions for audio?

texts:
The term non-lossy codec means the coding thing having no compression lost, more here. I need such format to process the audio, without destroying the quality. Wikipedia mentions in the context of lossy compression:

Ogg Vorbis (noted for its lack of patent restrictions)
Speex (noted for its lack of patent restrictions)

but the other wikipedia site for lossless things is open-ended, it mentions some patents and some expirations but not whether some codec has no patent-restrictions. Does there exist a lossless (non-lossy) codec without patent-restrictions or some sort of lists/matrix summing up the patent -things?

Perhaps related

Xiph.Org, the organization behind the multimedia formats and software tools here


texts: What does "convergence" mean for Thinning?

texts:
I've been using Thinning recently on line-like structures in an attempt to generate skeleton-like chains of pixels connected through their Moore neighborhoods.

The description of Thinning says that the algorithm "finds the skeletons of foreground regions in image by applying morphological thinning until convergence".  What does "convergence" mean?  If I start with a cluster of pixels connected through their Moore neighborhoods, what am I guaranteed after "convergence"?  Will the product always be a chain of pixels also connected through their Moore neighborhoods?